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