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