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 ) {
 208     char outfile[256];
 209     /* Generate output file name */
 210
 211
 212     /* Dump the "A" record */
 213
 214     /* get the name field (144 characters) */
 215     for ( i = 0; i < 144; i++ ) {
 216         name[i] = fgetc(fd);
 217     }
 218     name[i+1] = '\0';
 219
 220     /* clean off the whitespace at the end */
 221     for ( i = strlen(name)-2; i > 0; i-- ) {
 222         if ( !isspace(name[i]) ) {
 223             i=0;
 224         } else {
 225             name[i] = '\0';
 226         }
 227     }
 228     printf("    Quad name field: %s\n", name);
 229
 230     /* get quadrangle id (now part of previous section  */
 231     /* next_token(fd, dem_quadrangle); */
 232     /* printf("    Quadrangle = %s\n", dem_quadrangle); */
 233
 234     /* DEM level code, 3 reflects processing by DMA */
 235     inum = next_int(fd);
 236     printf("    DEM level code = %d\n", inum);
 237
 238     /* Pattern code, 1 indicates a regular elevation pattern */
 239     inum = next_int(fd);
 240     printf("    Pattern code = %d\n", inum);
 241
 242     /* Planimetric reference system code, 0 indicates geographic
 243      * coordinate system. */
 244     inum = next_int(fd);
 245     printf("    Planimetric reference code = %d\n", inum);
 246
 247     /* Zone code */
 248     inum = next_int(fd);
 249     printf("    Zone code = %d\n", inum);
 250
 251     /* Map projection parameters (ignored) */
 252     for ( i = 0; i < 15; i++ ) {
 253         dnum = next_double(fd);
 254         /* printf("%d: %f\n",i,dnum); */
 255     }
 256
 257     /* Units code, 3 represents arc-seconds as the unit of measure for
 258      * ground planimetric coordinates throughout the file. */
 259     inum = next_int(fd);
 260     if ( inum != 3 ) {
 261         printf("    Unknown (X,Y) units code = %d!\n", inum);
 262         exit(-1);
 263     }
 264
 265     /* Units code; 2 represents meters as the unit of measure for
 266      * elevation coordinates throughout the file. */
 267     inum = next_int(fd);
 268     if ( inum != 2 ) {
 269         printf("    Unknown (Z) units code = %d!\n", inum);
 270         exit(-1);
 271     }
 272
 273     /* Number (n) of sides in the polygon which defines the coverage of
 274      * the DEM file (usually equal to 4). */
 275     inum = next_int(fd);
 276     if ( inum != 4 ) {
 277         printf("    Unknown polygon dimension = %d!\n", inum);
 278         exit(-1);
 279     }
 280
 281     /* Ground coordinates of bounding box in arc-seconds */
 282     dem_x1 = m->originx = next_exp(fd);
 283     dem_y1 = m->originy = next_exp(fd);
 284     printf("    Origin = (%.2f,%.2f)\n", m->originx, m->originy);
 285
 286     dem_x2 = next_exp(fd);
 287     dem_y2 = next_exp(fd);
 288
 289     dem_x3 = next_exp(fd);
 290     dem_y3 = next_exp(fd);
 291
 292     dem_x4 = next_exp(fd);
 293     dem_y4 = next_exp(fd);
 294
 295     /* Minimum/maximum elevations in meters */
 296     dem_z1 = next_exp(fd);
 297     dem_z2 = next_exp(fd);
 298     printf("    Elevation range %.4f %.4f\n", dem_z1, dem_z2);
 299
 300     /* Counterclockwise angle from the primary axis of ground
 301      * planimetric referenced to the primary axis of the DEM local
 302      * reference system. */
 303     next_token(fd, token);
 304
 305     /* Accuracy code; 0 indicates that a record of accuracy does not
 306      * exist and that no record type C will follow. */
 307     /* DEM spacial resolution.  Usually (3,3,1) (3,6,1) or (3,9,1)
 308      * depending on latitude */
 309     /* I will eventually have to do something with this for data at
 310      * higher latitudes */
 311     next_token(fd, token);
 312     printf("    accuracy & spacial resolution string = %s\n", token);
 313     i = strlen(token);
 314     printf("    length = %d\n", i);
 315
 316     ptr = token + i - 12;
 317     printf("    last field = %s = %.2f\n", ptr, atof(ptr));
 318     ptr[0] = '\0';
 319
 320     ptr = ptr - 12;
 321     m->col_step = atof(ptr);
 322     printf("    last field = %s = %.2f\n", ptr, m->row_step);
 323     ptr[0] = '\0';
 324
 325     ptr = ptr - 12;
 326     m->row_step = atof(ptr);
 327     printf("    last field = %s = %.2f\n", ptr, m->col_step);
 328     ptr[0] = '\0';
 329
 330     /* accuracy code = atod(token) */
 331     inum = atoi(token);
 332     printf("    Accuracy code = %d\n", inum);
 333
 334     printf("    column step = %.2f  row step = %.2f\n",
 335            m->col_step, m->row_step);
 336     /* dimension of arrays to follow (1)*/
 337     next_token(fd, token);
 338
 339     /* number of profiles */
 340     dem_num_profiles = m->cols = m->rows = next_int(fd);
 341     printf("    Expecting %d profiles\n", dem_num_profiles);
 342
 343 }
 344
 345
 346 /* Read a horizontal strip of (1 vertical degree) from the raw DEM
 347  * file specified by the upper latitude of the stripe specified in
 348  * degrees */
 349 void rawReadStrip( fgRAWDEM *raw, int lat_degrees ) {
 350     int lat, yrange;
 351     int i, j, index, row, col;
 352     int min = 0, max = 0;
 353     int span, num_degrees, tile_width;
 354     int xstart, xend;
 355
 356     /* convert to arcsec */
 357     lat = lat_degrees * 3600;
 358
 359     printf("Max Latitude = %d arcsec\n", lat);
 360
 361     /* validity check ... */
 362     if ( (lat > raw->rooty) ||
 363          (lat < (raw->rooty - raw->nrows * raw->ydim + 1)) ) {
 364         printf("Latitude out of range for this DEM file\n");
 365         return;
 366     }
 367
 368     printf ("Reading strip starting at %d (top and working down)\n", lat);
 369
 370     /* advance to the correct latitude */
 371     rawAdvancePosition(raw, (raw->rooty - lat));
 372
 373     /* printf("short = %d\n", sizeof(short)); */
 374
 375     yrange = 3600 / raw->ydim;
 376
 377     for ( i = 0; i < yrange; i++ ) {
 378         index = yrange - i - 1;
 379         /* printf("About to read into row %d\n", index); */
 380         rawReadNextRow(raw, index);
 381
 382         for ( j = 0; j < raw->ncols; j++ ) {
 383             if ( raw->strip[index][j] == -9999 ) {
 384                 /* map ocean to 0 for now */
 385                 raw->strip[index][j] = 0;
 386             }
 387
 388             if ( raw->strip[index][j] < min) {
 389                 min = raw->strip[index][j];
 390             }
 391
 392             if ( raw->strip[index][j] > max) {
 393                 max = raw->strip[index][j];
 394             }
 395         }
 396     }
 397     printf("min = %d  max = %d\n", min, max);
 398
 399     /* extract individual tiles from the strip */
 400     span = raw->ncols * raw->xdim;
 401     num_degrees = span / 3600;
 402     tile_width = raw->ncols / num_degrees;
 403     printf("span = %d  num_degrees = %d  width = %d\n",
 404            span, num_degrees, tile_width);
 405
 406     for ( i = 0; i < num_degrees; i++ ) {
 407         xstart = i * tile_width;
 408         xend = xstart + 120;
 409
 410         for ( row = 0; row < yrange; row++ ) {
 411             for ( col = xstart; col < xend; col++ ) {
 412                 /* Copy from strip to pixel centered tile.  Yep,
 413                  * row/col are reversed here.  raw->strip is backwards
 414                  * for convenience.  I am converting to [x,y] now. */
 415                 raw->center[col-xstart][row] = raw->strip[row][col];
 416             }
 417         }
 418
 419         /* Convert from pixel centered to pixel edge values */
 420         rawConvertCenter2Edge(raw);
 421
 422         /* Dump out the ascii format DEM file */
 423         rawDumpAsciiDEM(raw);
 424     }
 425 }
 426
 427
 428 /* $Log$
 429 /* Revision 1.2  1998/03/03 02:04:01  curt
 430 /* Starting DEM Ascii format output routine.
 431 /*
 432  * Revision 1.1  1998/03/02 23:31:01  curt
 433  * Initial revision.
 434  *
 435  */