1 /* rawdem.c -- library of routines for processing raw dem files (30 arcsec)
3 * Written by Curtis Olson, started February 1998.
5 * Copyright (C) 1998 Curtis L. Olson - curt@me.umn.edu
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.
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.
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.
22 * (Log is kept at end of this file)
30 #include <math.h> /* rint() */
32 #include <stdlib.h> /* atoi() atof() */
33 #include <string.h> /* swab() */
35 #include <sys/types.h> /* open() */
38 #include <unistd.h> /* close() */
43 /* Read the DEM header to determine various key parameters for this
45 void rawReadDemHdr( fgRAWDEM *raw, char *hdr_file ) {
47 char line[256], key[256], value[256];
51 if ( (hdr = fopen(hdr_file, "r")) == NULL ) {
52 printf("Error opening DEM header file: %s\n", hdr_file);
56 /* process each line */
57 while ( (fgets(line, 256, hdr) != NULL) ) {
58 /* printf("%s", line); */
63 while ( (line[i] != ' ') && (i < len) ) {
69 /* skip middle space */
70 while ( (line[i] == ' ') && (i < len) ) {
76 while ( (line[i] != '\n') && (i < len) ) {
77 value[i-offset] = line[i];
80 value[i-offset] = '\0';
81 /* printf("key='%s' value='%s'\n", key, value); */
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 ) {
90 raw->ulxmap = (int)rint(tmp * 3600.0); /* convert to arcsec */
92 # error Port me rint()
94 } else if ( strcmp(key, "ULYMAP") == 0 ) {
97 raw->ulymap = (int)rint(tmp * 3600.0); /* convert to arcsec */
99 # error Port me rint()
101 } else if ( strcmp(key, "XDIM") == 0 ) {
104 raw->xdim = (int)rint(tmp * 3600.0); /* convert to arcsec */
106 # error Port me rint()
108 } else if ( strcmp(key, "YDIM") == 0 ) {
111 raw->ydim = (int)rint(tmp * 3600.0); /* convert to arcsec */
113 # error Port me rint()
120 raw->rootx = raw->ulxmap - (raw->xdim / 2);
121 raw->rooty = raw->ulymap + (raw->ydim / 2);
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,
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);
139 /* Close a raw DEM file. */
140 void rawCloseDemFile( fgRAWDEM *raw ) {
145 /* Advance file pointer position to correct latitude (row) */
146 void rawAdvancePosition( fgRAWDEM *raw, int arcsec ) {
149 offset = 2 * raw->ncols * ( arcsec / raw->ydim );
151 if ( (result = lseek(raw->fd, offset, SEEK_SET)) == -1 ) {
152 printf("Error lseek filed trying to offset by %ld\n", offset);
156 printf("Successful seek ahead of %ld bytes\n", result);
160 /* Read the next row of data */
161 void rawReadNextRow( fgRAWDEM *raw, int index ) {
162 char buf[MAX_COLS_X_2];
165 if ( raw->ncols > MAX_ROWS ) {
166 printf("Error, buf needs to be bigger in rawReadNextRow()\n");
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); */
174 /* reverse byte order */
175 /* it would be nice to test in advance some how if we need to do
177 /* swab(frombuf, tobuf, 2 * raw->ncols); */
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];
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:
197 void rawConvertCenter2Edge( fgRAWDEM *raw ) {
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];
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;
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] +
220 raw->center[i-1][j] ) / 4;
226 /* Dump out the ascii format DEM file */
227 void rawDumpAsciiDEM( fgRAWDEM *raw, char *path, int ilon, int ilat ) {
231 char lon_sign, lat_sign;
235 /* Generate output file name */
253 sprintf(outfile, "%s/%c%03d%c%03d.dem", path, lon_sign, lon, lat_sign, lat);
255 printf("outfile = %s\n", outfile);
257 if ( (fd = fopen(outfile, "w")) == NULL ) {
258 printf("Error opening output file = %s\n", outfile);
262 /* Dump the "A" record */
264 /* print descriptive header (144 characters) */
265 sprintf(tmp, "%s - Generated from a 30 arcsec binary DEM", outfile);
266 fprintf(fd, "%-144s", tmp);
268 /* DEM level code, 3 reflects processing by DMA */
269 fprintf(fd, "%6d", 1);
271 /* Pattern code, 1 indicates a regular elevation pattern */
272 fprintf(fd, "%6d", 1);
274 /* Planimetric reference system code, 0 indicates geographic
275 * coordinate system. */
276 fprintf(fd, "%6d", 0);
279 fprintf(fd, "%6d", 0);
281 /* Map projection parameters (ignored) */
282 for ( i = 0; i < 15; i++ ) {
283 fprintf(fd, "%6.1f%18s", 0.0, "");
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);
290 /* Units code; 2 represents meters as the unit of measure for
291 * elevation coordinates throughout the file. */
292 fprintf(fd, "%6d", 2);
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);
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);
302 fprintf(fd, "%20.15fD+06", ilon * 3600.0 / 1000000.0);
303 fprintf(fd, "%20.15fD+06", (ilat+1) * 3600.0 / 1000000.0);
305 fprintf(fd, "%20.15fD+06", (ilon+1) * 3600.0 / 1000000.0);
306 fprintf(fd, "%20.15fD+06", (ilat+1) * 3600.0 / 1000000.0);
308 fprintf(fd, "%20.15fD+06", (ilon+1) * 3600.0 / 1000000.0);
309 fprintf(fd, "%20.15fD+06", (ilat) * 3600.0 / 1000000.0);
311 /* Minimum/maximum elevations in meters */
312 fprintf(fd, " %20.15E", (double)raw->tmp_min);
313 fprintf(fd, " %20.15E", (double)raw->tmp_max);
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);
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);
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);
330 fprintf(fd, "%12.6E", 1.0);
332 /* dimension of arrays to follow (1)*/
333 fprintf(fd, "%6d", 1);
335 /* number of profiles */
336 fprintf(fd, "%6d", 3600 / raw->ydim + 1);
339 fprintf(fd, "%160s", "");
342 /* Dump "B" records */
344 for ( j = 0; j <= 120; j++ ) {
345 /* row / column id of this profile */
346 fprintf(fd, "%6d%6d", 1, j + 1);
348 /* Number of rows and columns (elevation points) in this
350 fprintf(fd, "%6d%6d", 3600 / raw->xdim + 1, 1);
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);
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", "");
362 /* Minimum and maximum elevations for the profile. */
363 fprintf(fd, " %20.15E", 0.0);
364 fprintf(fd, " %20.15E", 0.0);
366 /* One (usually) dimensional array (1,prof_num_cols) of
368 for ( i = 0; i <= 120; i++ ) {
369 fprintf(fd, "%6.0f", raw->edge[j][i]);
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 ) {
384 int i, j, index, row, col;
386 int span, num_degrees, tile_width;
389 /* convert to arcsec */
390 lat = lat_degrees * 3600;
392 printf("Max Latitude = %d arcsec\n", lat);
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");
401 printf ("Reading strip starting at %d (top and working down)\n", lat);
403 /* advance to the correct latitude */
404 rawAdvancePosition(raw, (raw->rooty - lat));
406 /* printf("short = %d\n", sizeof(short)); */
408 yrange = 3600 / raw->ydim;
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);
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;
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);
430 for ( i = 0; i < num_degrees; i++ ) {
431 xstart = i * tile_width;
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];
442 if ( raw->strip[row][col] < min) {
443 min = raw->strip[row][col];
446 if ( raw->strip[row][col] > max) {
447 max = raw->strip[row][col];
455 /* Convert from pixel centered to pixel edge values */
456 rawConvertCenter2Edge(raw);
458 /* Dump out the ascii format DEM file */
459 rawDumpAsciiDEM(raw, path, (raw->rootx / 3600) + i, lat_degrees - 1);
465 /* Revision 1.6 1998/04/27 03:32:03 curt
466 /* Wrapped rint()'s in #ifdef HAVE_RINT
468 * Revision 1.5 1998/04/18 03:59:46 curt
469 * Incorporated into gnu automake/autoconf system.
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.
476 * Revision 1.3 1998/03/03 13:10:29 curt
477 * Close to a working version.
479 * Revision 1.2 1998/03/03 02:04:01 curt
480 * Starting DEM Ascii format output routine.
482 * Revision 1.1 1998/03/02 23:31:01 curt