Tools/Prep/DemRaw2ascii/rawdem.c

   1 /* rawdem.c -- library of routines for processing raw dem files (30 arcsec)
   2  *
   3  * Written by Curtis Olson, started February 1998.
   4  *
   5  * Copyright (C) 1998  Curtis L. Olson  - curt@me.umn.edu
   6  *
   7  * This program is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 2 of the License, or
  10  * (at your option) any later version.
  11  *
  12  * This program is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with this program; if not, write to the Free Software
  19  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  20  *
  21  * $Id$
  22  */
  23
  24
  25 #ifdef HAVE_CONFIG_H
  26 #  include <config.h>
  27 #endif
  28
  29 #include <math.h>      /* rint() */
  30 #include <stdio.h>
  31 #include <stdlib.h>    /* atoi() atof() */
  32 #include <string.h>    /* swab() */
  33
  34 #include <sys/types.h> /* open() */
  35 #include <sys/stat.h>
  36 #include <fcntl.h>
  37 #include <unistd.h>    /* close() */
  38
  39 #include "rawdem.h"
  40
  41
  42 /* Read the DEM header to determine various key parameters for this
  43  * DEM file */
  44 void rawReadDemHdr( fgRAWDEM *raw, char *hdr_file ) {
  45     FILE *hdr;
  46     char line[256], key[256], value[256];
  47     int i, len, offset;
  48     double tmp;
  49
  50     if ( (hdr = fopen(hdr_file, "r")) == NULL ) {
  51         printf("Error opening DEM header file: %s\n", hdr_file);
  52         exit(-1);
  53     }
  54
  55     /* process each line */
  56     while ( (fgets(line, 256, hdr) != NULL) ) {
  57         /* printf("%s", line); */
  58         len = strlen(line);
  59
  60         /* extract key */
  61         i = 0;
  62         while ( (line[i] != ' ') && (i < len) ) {
  63             key[i] = line[i];
  64             i++;
  65         }
  66         key[i] = '\0';
  67
  68         /* skip middle space */
  69         while ( (line[i] == ' ') && (i < len) ) {
  70             i++;
  71         }
  72         offset = i;
  73
  74         /* extract value */
  75         while ( (line[i] != '\n') && (i < len) ) {
  76             value[i-offset] = line[i];
  77             i++;
  78         }
  79         value[i-offset] = '\0';
  80         /* printf("key='%s'  value='%s'\n", key, value); */
  81
  82         if ( strcmp(key, "NROWS") == 0 ) {
  83             raw->nrows = atoi(value);
  84         } else if ( strcmp(key, "NCOLS") == 0 ) {
  85             raw->ncols = atoi(value);
  86         } else if ( strcmp(key, "ULXMAP") == 0 ) {
  87             tmp = atof(value);
  88 #ifdef HAVE_RINT
  89             raw->ulxmap = (int)rint(tmp * 3600.0); /* convert to arcsec */
  90 #else
  91 #  error Port me rint()
  92 #endif
  93         } else if ( strcmp(key, "ULYMAP") == 0 ) {
  94             tmp = atof(value);
  95 #ifdef HAVE_RINT
  96             raw->ulymap = (int)rint(tmp * 3600.0); /* convert to arcsec */
  97 #else
  98 #  error Port me rint()
  99 #endif
 100         } else if ( strcmp(key, "XDIM") == 0 ) {
 101             tmp = atof(value);
 102 #ifdef HAVE_RINT
 103             raw->xdim = (int)rint(tmp * 3600.0);   /* convert to arcsec */
 104 #else
 105 #  error Port me rint()
 106 #endif
 107         } else if ( strcmp(key, "YDIM") == 0 ) {
 108             tmp = atof(value);
 109 #ifdef HAVE_RINT
 110             raw->ydim = (int)rint(tmp * 3600.0);   /* convert to arcsec */
 111 #else
 112 #  error Port me rint()
 113 #endif
 114         } else {
 115             /* ignore for now */
 116         }
 117     }
 118
 119     raw->rootx = raw->ulxmap - (raw->xdim / 2);
 120     raw->rooty = raw->ulymap + (raw->ydim / 2);
 121
 122     printf("%d %d %d %d %d %d %d %d\n", raw->nrows, raw->ncols,
 123            raw->ulxmap, raw->ulymap, raw->rootx, raw->rooty, raw->xdim,
 124            raw->ydim);
 125 }
 126
 127
 128 /* Open a raw DEM file. */
 129 void rawOpenDemFile( fgRAWDEM *raw, char *raw_dem_file ) {
 130     printf("Opening Raw DEM file: %s\n", raw_dem_file);
 131     if ( (raw->fd = open(raw_dem_file ,O_RDONLY)) == -1 ) {
 132         printf("Error opening Raw DEM file: %s\n", raw_dem_file);
 133         exit(-1);
 134     }
 135 }
 136
 137
 138 /* Close a raw DEM file. */
 139 void rawCloseDemFile( fgRAWDEM *raw ) {
 140     close(raw->fd);
 141 }
 142
 143
 144 /* Advance file pointer position to correct latitude (row) */
 145 void rawAdvancePosition( fgRAWDEM *raw, int arcsec ) {
 146     long offset, result;
 147
 148     offset = 2 * raw->ncols * ( arcsec  / raw->ydim );
 149
 150     if ( (result = lseek(raw->fd, offset, SEEK_SET)) == -1 ) {
 151         printf("Error lseek filed trying to offset by %ld\n", offset);
 152         exit(-1);
 153     }
 154
 155     printf("Successful seek ahead of %ld bytes\n", result);
 156 }
 157
 158
 159 /* Read the next row of data */
 160 void rawReadNextRow( fgRAWDEM *raw, int index ) {
 161     char buf[MAX_COLS_X_2];
 162     int i, result;
 163
 164     if ( raw->ncols > MAX_ROWS ) {
 165         printf("Error, buf needs to be bigger in rawReadNextRow()\n");
 166         exit(-1);
 167     }
 168
 169     /* printf("Attempting to read %d bytes\n", 2 * raw->ncols); */
 170     result = read(raw->fd, buf, 2 * raw->ncols);
 171     /* printf("Read %d bytes\n", result); */
 172
 173     /* reverse byte order */
 174     /* it would be nice to test in advance some how if we need to do
 175      * this */
 176     /* swab(frombuf, tobuf, 2 * raw->ncols); */
 177
 178     for ( i = 0; i < raw->ncols; i++ ) {
 179         /* printf("hi = %d  lo = %d\n", buf[2*i], buf[2*i + 1]); */
 180         raw->strip[index][i] = ( (buf[2*i] + 1) << 8 ) + buf[2*i + 1];
 181     }
 182 }
 183
 184
 185 /* Convert from pixel centered values to pixel corner values.  This is
 186    accomplished by taking the average of the closes center nodes.  In
 187    the following diagram "x" marks the data point location:
 188
 189    +-----+        x-----x
 190    |     |        |     |
 191    |  x  |  ===>  |     |
 192    |     |        |     |
 193    +-----+        x-----x
 194
 195    */
 196 void rawConvertCenter2Edge( fgRAWDEM *raw ) {
 197     int i, j;
 198
 199     /* derive corner nodes */
 200     raw->edge[0][0]     = raw->center[0][0];
 201     raw->edge[120][0]   = raw->center[119][0];
 202     raw->edge[120][120] = raw->center[119][119];
 203     raw->edge[0][120]   = raw->center[0][119];
 204
 205     /* derive edge nodes */
 206     for ( i = 1; i < 120; i++ ) {
 207         raw->edge[i][0] = (raw->center[i-1][0] + raw->center[i][0]) / 2.0;
 208         raw->edge[i][120] = (raw->center[i-1][119] + raw->center[i][119]) / 2.0;
 209         raw->edge[0][i] = (raw->center[0][i-1] + raw->center[0][i]) / 2.0;
 210         raw->edge[120][i] = (raw->center[119][i-1] + raw->center[119][i]) / 2.0;
 211     }
 212
 213     /* derive internal nodes */
 214     for ( j = 1; j < 120; j++ ) {
 215         for ( i = 1; i < 120; i++ ) {
 216             raw->edge[i][j] = ( raw->center[i-1][j-1] +
 217                                 raw->center[i]  [j-1] +
 218                                 raw->center[i]  [j]   +
 219                                 raw->center[i-1][j] ) / 4;
 220         }
 221     }
 222 }
 223
 224
 225 /* Dump out the ascii format DEM file */
 226 void rawDumpAsciiDEM( fgRAWDEM *raw, char *path, int ilon, int ilat ) {
 227     char outfile[256];
 228     char tmp[256];
 229     int lon, lat;
 230     char lon_sign, lat_sign;
 231     int i, j;
 232     FILE *fd;
 233
 234     /* Generate output file name */
 235
 236     if ( ilon >= 0 ) {
 237         lon = ilon;
 238         lon_sign = 'e';
 239     } else {
 240         lon = -ilon;
 241         lon_sign = 'w';
 242     }
 243
 244     if ( ilat >= 0 ) {
 245         lat = ilat;
 246         lat_sign = 'n';
 247     } else {
 248         lat = -ilat;
 249         lat_sign = 's';
 250     }
 251
 252     sprintf(outfile, "%s/%c%03d%c%03d.dem", path, lon_sign, lon, lat_sign, lat);
 253
 254     printf("outfile = %s\n", outfile);
 255
 256     if ( (fd = fopen(outfile, "w")) == NULL ) {
 257         printf("Error opening output file = %s\n", outfile);
 258         exit(-1);
 259     }
 260
 261     /* Dump the "A" record */
 262
 263     /* print descriptive header (144 characters) */
 264     sprintf(tmp, "%s - Generated from a 30 arcsec binary DEM", outfile);
 265     fprintf(fd, "%-144s", tmp);
 266
 267     /* DEM level code, 3 reflects processing by DMA */
 268     fprintf(fd, "%6d", 1);
 269
 270     /* Pattern code, 1 indicates a regular elevation pattern */
 271     fprintf(fd, "%6d", 1);
 272
 273     /* Planimetric reference system code, 0 indicates geographic
 274      * coordinate system. */
 275     fprintf(fd, "%6d", 0);
 276
 277     /* Zone code */
 278     fprintf(fd, "%6d", 0);
 279
 280     /* Map projection parameters (ignored) */
 281     for ( i = 0; i < 15; i++ ) {
 282         fprintf(fd, "%6.1f%18s", 0.0, "");
 283     }
 284
 285    /* Units code, 3 represents arc-seconds as the unit of measure for
 286      * ground planimetric coordinates throughout the file. */
 287     fprintf(fd, "%6d", 3);
 288
 289     /* Units code; 2 represents meters as the unit of measure for
 290      * elevation coordinates throughout the file. */
 291     fprintf(fd, "%6d", 2);
 292
 293     /* Number (n) of sides in the polygon which defines the coverage of
 294      * the DEM file (usually equal to 4). */
 295     fprintf(fd, "%6d", 4);
 296
 297     /* Ground coordinates of bounding box in arc-seconds */
 298     fprintf(fd, "%20.15fD+06", ilon * 3600.0 / 1000000.0);
 299     fprintf(fd, "%20.15fD+06", ilat * 3600.0 / 1000000.0);
 300
 301     fprintf(fd, "%20.15fD+06", ilon * 3600.0 / 1000000.0);
 302     fprintf(fd, "%20.15fD+06", (ilat+1) * 3600.0 / 1000000.0);
 303
 304     fprintf(fd, "%20.15fD+06", (ilon+1) * 3600.0 / 1000000.0);
 305     fprintf(fd, "%20.15fD+06", (ilat+1) * 3600.0 / 1000000.0);
 306
 307     fprintf(fd, "%20.15fD+06", (ilon+1) * 3600.0 / 1000000.0);
 308     fprintf(fd, "%20.15fD+06", (ilat) * 3600.0 / 1000000.0);
 309
 310     /* Minimum/maximum elevations in meters */
 311     fprintf(fd, "   %20.15E", (double)raw->tmp_min);
 312     fprintf(fd, "   %20.15E", (double)raw->tmp_max);
 313
 314     /* Counterclockwise angle from the primary axis of ground
 315      * planimetric referenced to the primary axis of the DEM local
 316      * reference system. */
 317     fprintf(fd, "%6.1f", 0.0);
 318
 319     /* Accuracy code; 0 indicates that a record of accuracy does not
 320      * exist and that no record type C will follow. */
 321     fprintf(fd, "%24d", 0);
 322
 323     /* DEM spacial resolution.  Usually (3,3) (3,6) or (3,9)
 324      * depending on latitude */
 325     fprintf(fd, "%12.6E", 30.0);
 326     fprintf(fd, "%12.6E", 30.0);
 327
 328     /* accuracy code */
 329     fprintf(fd, "%12.6E", 1.0);
 330
 331     /* dimension of arrays to follow (1)*/
 332     fprintf(fd, "%6d", 1);
 333
 334     /* number of profiles */
 335     fprintf(fd, "%6d", 3600 / raw->ydim + 1);
 336
 337     /* pad the end */
 338     fprintf(fd, "%160s", "");
 339
 340
 341     /* Dump "B" records */
 342
 343     for ( j = 0; j <= 120; j++ ) {
 344         /* row / column id of this profile */
 345         fprintf(fd, "%6d%6d", 1, j + 1);
 346
 347         /* Number of rows and columns (elevation points) in this
 348            profile */
 349         fprintf(fd, "%6d%6d", 3600 / raw->xdim + 1, 1);
 350
 351         /* Ground planimetric coordinates (arc-seconds) of the first
 352          * elevation in the profile */
 353         fprintf(fd, "%20.15fD+06", ilon * 3600.0 / 1000000.0);
 354         fprintf(fd, "%20.15fD+06", (ilat * 3600.0 + j * raw->ydim) / 1000000.0);
 355
 356         /* Elevation of local datum for the profile.  Always zero for
 357          * 1-degree DEM, the reference is mean sea level. */
 358         fprintf(fd, "%6.1f", 0.0);
 359         fprintf(fd, "%18s", "");
 360
 361         /* Minimum and maximum elevations for the profile. */
 362         fprintf(fd, "   %20.15E", 0.0);
 363         fprintf(fd, "   %20.15E", 0.0);
 364
 365         /* One (usually) dimensional array (1,prof_num_cols) of
 366            elevations */
 367         for ( i = 0; i <= 120; i++ ) {
 368             fprintf(fd, "%6.0f", raw->edge[j][i]);
 369         }
 370     }
 371
 372     fprintf(fd, "\n");
 373
 374     fclose(fd);
 375 }
 376
 377
 378 /* Read a horizontal strip of (1 vertical degree) from the raw DEM
 379  * file specified by the upper latitude of the stripe specified in
 380  * degrees.  The output the individual ASCII format DEM tiles.  */
 381 void rawProcessStrip( fgRAWDEM *raw, int lat_degrees, char *path ) {
 382     int lat, yrange;
 383     int i, j, index, row, col;
 384     int min, max;
 385     int span, num_degrees, tile_width;
 386     int xstart, xend;
 387
 388     /* convert to arcsec */
 389     lat = lat_degrees * 3600;
 390
 391     printf("Max Latitude = %d arcsec\n", lat);
 392
 393     /* validity check ... */
 394     if ( (lat > raw->rooty) ||
 395          (lat < (raw->rooty - raw->nrows * raw->ydim + 1)) ) {
 396         printf("Latitude out of range for this DEM file\n");
 397         return;
 398     }
 399
 400     printf ("Reading strip starting at %d (top and working down)\n", lat);
 401
 402     /* advance to the correct latitude */
 403     rawAdvancePosition(raw, (raw->rooty - lat));
 404
 405     /* printf("short = %d\n", sizeof(short)); */
 406
 407     yrange = 3600 / raw->ydim;
 408
 409     for ( i = 0; i < yrange; i++ ) {
 410         index = yrange - i - 1;
 411         /* printf("About to read into row %d\n", index); */
 412         rawReadNextRow(raw, index);
 413
 414         for ( j = 0; j < raw->ncols; j++ ) {
 415             if ( raw->strip[index][j] == -9999 ) {
 416                 /* map ocean to 0 for now */
 417                 raw->strip[index][j] = 0;
 418             }
 419         }
 420     }
 421
 422     /* extract individual tiles from the strip */
 423     span = raw->ncols * raw->xdim;
 424     num_degrees = span / 3600;
 425     tile_width = raw->ncols / num_degrees;
 426     printf("span = %d  num_degrees = %d  width = %d\n",
 427            span, num_degrees, tile_width);
 428
 429     for ( i = 0; i < num_degrees; i++ ) {
 430         xstart = i * tile_width;
 431         xend = xstart + 120;
 432
 433         min = 10000; max = -10000;
 434         for ( row = 0; row < yrange; row++ ) {
 435             for ( col = xstart; col < xend; col++ ) {
 436                 /* Copy from strip to pixel centered tile.  Yep,
 437                  * row/col are reversed here.  raw->strip is backwards
 438                  * for convenience.  I am converting to [x,y] now. */
 439                 raw->center[col-xstart][row] = raw->strip[row][col];
 440
 441                 if ( raw->strip[row][col] < min) {
 442                     min = raw->strip[row][col];
 443                 }
 444
 445                 if ( raw->strip[row][col] > max) {
 446                     max = raw->strip[row][col];
 447                 }
 448             }
 449         }
 450
 451         raw->tmp_min = min;
 452         raw->tmp_max = max;
 453
 454         /* Convert from pixel centered to pixel edge values */
 455         rawConvertCenter2Edge(raw);
 456
 457         /* Dump out the ascii format DEM file */
 458         rawDumpAsciiDEM(raw, path, (raw->rootx / 3600) + i, lat_degrees - 1);
 459     }
 460 }
 461
 462