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 }
 209
 210
 211 /* Read a horizontal strip of (1 vertical degree) from the raw DEM
 212  * file specified by the upper latitude of the stripe specified in
 213  * degrees */
 214 void rawReadStrip( fgRAWDEM *raw, int lat_degrees ) {
 215     int lat, yrange;
 216     int i, j, index, row, col;
 217     int min = 0, max = 0;
 218     int span, num_degrees, tile_width;
 219     int xstart, xend;
 220
 221     /* convert to arcsec */
 222     lat = lat_degrees * 3600;
 223
 224     printf("Max Latitude = %d arcsec\n", lat);
 225
 226     /* validity check ... */
 227     if ( (lat > raw->rooty) ||
 228          (lat < (raw->rooty - raw->nrows * raw->ydim + 1)) ) {
 229         printf("Latitude out of range for this DEM file\n");
 230         return;
 231     }
 232
 233     printf ("Reading strip starting at %d (top and working down)\n", lat);
 234
 235     /* advance to the correct latitude */
 236     rawAdvancePosition(raw, (raw->rooty - lat));
 237
 238     /* printf("short = %d\n", sizeof(short)); */
 239
 240     yrange = 3600 / raw->ydim;
 241
 242     for ( i = 0; i < yrange; i++ ) {
 243         index = yrange - i - 1;
 244         /* printf("About to read into row %d\n", index); */
 245         rawReadNextRow(raw, index);
 246
 247         for ( j = 0; j < raw->ncols; j++ ) {
 248             if ( raw->strip[index][j] == -9999 ) {
 249                 /* map ocean to 0 for now */
 250                 raw->strip[index][j] = 0;
 251             }
 252
 253             if ( raw->strip[index][j] < min) {
 254                 min = raw->strip[index][j];
 255             }
 256
 257             if ( raw->strip[index][j] > max) {
 258                 max = raw->strip[index][j];
 259             }
 260         }
 261     }
 262     printf("min = %d  max = %d\n", min, max);
 263
 264     /* extract individual tiles from the strip */
 265     span = raw->ncols * raw->xdim;
 266     num_degrees = span / 3600;
 267     tile_width = raw->ncols / num_degrees;
 268     printf("span = %d  num_degrees = %d  width = %d\n",
 269            span, num_degrees, tile_width);
 270
 271     for ( i = 0; i < num_degrees; i++ ) {
 272         xstart = i * tile_width;
 273         xend = xstart + 120;
 274
 275         for ( row = 0; row < yrange; row++ ) {
 276             for ( col = xstart; col < xend; col++ ) {
 277                 /* Copy from strip to pixel centered tile.  Yep,
 278                  * row/col are reversed here.  raw->strip is backwards
 279                  * for convenience.  I am converting to [x,y] now. */
 280                 raw->center[col-xstart][row] = raw->strip[row][col];
 281             }
 282         }
 283
 284         /* Convert from pixel centered to pixel edge values */
 285         rawConvertCenter2Edge(raw);
 286
 287         /* Dump out the ascii format DEM file */
 288         rawDumpAsciiDEM(raw);
 289     }
 290 }
 291
 292
 293 /* $Log$
 294 /* Revision 1.1  1998/03/02 23:31:01  curt
 295 /* Initial revision.
 296 /*
 297  */