485 lines
13 KiB
C
485 lines
13 KiB
C
/* rawdem.c -- library of routines for processing raw dem files (30 arcsec)
|
|
*
|
|
* Written by Curtis Olson, started February 1998.
|
|
*
|
|
* Copyright (C) 1998 Curtis L. Olson - curt@me.umn.edu
|
|
*
|
|
* This program is free software; you can redistribute it and/or modify
|
|
* it under the terms of the GNU General Public License as published by
|
|
* the Free Software Foundation; either version 2 of the License, or
|
|
* (at your option) any later version.
|
|
*
|
|
* This program is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
* GNU General Public License for more details.
|
|
*
|
|
* You should have received a copy of the GNU General Public License
|
|
* along with this program; if not, write to the Free Software
|
|
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
|
|
*
|
|
* $Id$
|
|
* (Log is kept at end of this file)
|
|
*/
|
|
|
|
|
|
#ifdef HAVE_CONFIG_H
|
|
# include <config.h>
|
|
#endif
|
|
|
|
#include <math.h> /* rint() */
|
|
#include <stdio.h>
|
|
#include <stdlib.h> /* atoi() atof() */
|
|
#include <string.h> /* swab() */
|
|
|
|
#include <sys/types.h> /* open() */
|
|
#include <sys/stat.h>
|
|
#include <fcntl.h>
|
|
#include <unistd.h> /* close() */
|
|
|
|
#include "rawdem.h"
|
|
|
|
|
|
/* Read the DEM header to determine various key parameters for this
|
|
* DEM file */
|
|
void rawReadDemHdr( fgRAWDEM *raw, char *hdr_file ) {
|
|
FILE *hdr;
|
|
char line[256], key[256], value[256];
|
|
int i, len, offset;
|
|
double tmp;
|
|
|
|
if ( (hdr = fopen(hdr_file, "r")) == NULL ) {
|
|
printf("Error opening DEM header file: %s\n", hdr_file);
|
|
exit(-1);
|
|
}
|
|
|
|
/* process each line */
|
|
while ( (fgets(line, 256, hdr) != NULL) ) {
|
|
/* printf("%s", line); */
|
|
len = strlen(line);
|
|
|
|
/* extract key */
|
|
i = 0;
|
|
while ( (line[i] != ' ') && (i < len) ) {
|
|
key[i] = line[i];
|
|
i++;
|
|
}
|
|
key[i] = '\0';
|
|
|
|
/* skip middle space */
|
|
while ( (line[i] == ' ') && (i < len) ) {
|
|
i++;
|
|
}
|
|
offset = i;
|
|
|
|
/* extract value */
|
|
while ( (line[i] != '\n') && (i < len) ) {
|
|
value[i-offset] = line[i];
|
|
i++;
|
|
}
|
|
value[i-offset] = '\0';
|
|
/* printf("key='%s' value='%s'\n", key, value); */
|
|
|
|
if ( strcmp(key, "NROWS") == 0 ) {
|
|
raw->nrows = atoi(value);
|
|
} else if ( strcmp(key, "NCOLS") == 0 ) {
|
|
raw->ncols = atoi(value);
|
|
} else if ( strcmp(key, "ULXMAP") == 0 ) {
|
|
tmp = atof(value);
|
|
#ifdef HAVE_RINT
|
|
raw->ulxmap = (int)rint(tmp * 3600.0); /* convert to arcsec */
|
|
#else
|
|
# error Port me rint()
|
|
#endif
|
|
} else if ( strcmp(key, "ULYMAP") == 0 ) {
|
|
tmp = atof(value);
|
|
#ifdef HAVE_RINT
|
|
raw->ulymap = (int)rint(tmp * 3600.0); /* convert to arcsec */
|
|
#else
|
|
# error Port me rint()
|
|
#endif
|
|
} else if ( strcmp(key, "XDIM") == 0 ) {
|
|
tmp = atof(value);
|
|
#ifdef HAVE_RINT
|
|
raw->xdim = (int)rint(tmp * 3600.0); /* convert to arcsec */
|
|
#else
|
|
# error Port me rint()
|
|
#endif
|
|
} else if ( strcmp(key, "YDIM") == 0 ) {
|
|
tmp = atof(value);
|
|
#ifdef HAVE_RINT
|
|
raw->ydim = (int)rint(tmp * 3600.0); /* convert to arcsec */
|
|
#else
|
|
# error Port me rint()
|
|
#endif
|
|
} else {
|
|
/* ignore for now */
|
|
}
|
|
}
|
|
|
|
raw->rootx = raw->ulxmap - (raw->xdim / 2);
|
|
raw->rooty = raw->ulymap + (raw->ydim / 2);
|
|
|
|
printf("%d %d %d %d %d %d %d %d\n", raw->nrows, raw->ncols,
|
|
raw->ulxmap, raw->ulymap, raw->rootx, raw->rooty, raw->xdim,
|
|
raw->ydim);
|
|
}
|
|
|
|
|
|
/* Open a raw DEM file. */
|
|
void rawOpenDemFile( fgRAWDEM *raw, char *raw_dem_file ) {
|
|
printf("Opening Raw DEM file: %s\n", raw_dem_file);
|
|
if ( (raw->fd = open(raw_dem_file ,O_RDONLY)) == -1 ) {
|
|
printf("Error opening Raw DEM file: %s\n", raw_dem_file);
|
|
exit(-1);
|
|
}
|
|
}
|
|
|
|
|
|
/* Close a raw DEM file. */
|
|
void rawCloseDemFile( fgRAWDEM *raw ) {
|
|
close(raw->fd);
|
|
}
|
|
|
|
|
|
/* Advance file pointer position to correct latitude (row) */
|
|
void rawAdvancePosition( fgRAWDEM *raw, int arcsec ) {
|
|
long offset, result;
|
|
|
|
offset = 2 * raw->ncols * ( arcsec / raw->ydim );
|
|
|
|
if ( (result = lseek(raw->fd, offset, SEEK_SET)) == -1 ) {
|
|
printf("Error lseek filed trying to offset by %ld\n", offset);
|
|
exit(-1);
|
|
}
|
|
|
|
printf("Successful seek ahead of %ld bytes\n", result);
|
|
}
|
|
|
|
|
|
/* Read the next row of data */
|
|
void rawReadNextRow( fgRAWDEM *raw, int index ) {
|
|
char buf[MAX_COLS_X_2];
|
|
int i, result;
|
|
|
|
if ( raw->ncols > MAX_ROWS ) {
|
|
printf("Error, buf needs to be bigger in rawReadNextRow()\n");
|
|
exit(-1);
|
|
}
|
|
|
|
/* printf("Attempting to read %d bytes\n", 2 * raw->ncols); */
|
|
result = read(raw->fd, buf, 2 * raw->ncols);
|
|
/* printf("Read %d bytes\n", result); */
|
|
|
|
/* reverse byte order */
|
|
/* it would be nice to test in advance some how if we need to do
|
|
* this */
|
|
/* swab(frombuf, tobuf, 2 * raw->ncols); */
|
|
|
|
for ( i = 0; i < raw->ncols; i++ ) {
|
|
/* printf("hi = %d lo = %d\n", buf[2*i], buf[2*i + 1]); */
|
|
raw->strip[index][i] = ( (buf[2*i] + 1) << 8 ) + buf[2*i + 1];
|
|
}
|
|
}
|
|
|
|
|
|
/* Convert from pixel centered values to pixel corner values. This is
|
|
accomplished by taking the average of the closes center nodes. In
|
|
the following diagram "x" marks the data point location:
|
|
|
|
+-----+ x-----x
|
|
| | | |
|
|
| x | ===> | |
|
|
| | | |
|
|
+-----+ x-----x
|
|
|
|
*/
|
|
void rawConvertCenter2Edge( fgRAWDEM *raw ) {
|
|
int i, j;
|
|
|
|
/* derive corner nodes */
|
|
raw->edge[0][0] = raw->center[0][0];
|
|
raw->edge[120][0] = raw->center[119][0];
|
|
raw->edge[120][120] = raw->center[119][119];
|
|
raw->edge[0][120] = raw->center[0][119];
|
|
|
|
/* derive edge nodes */
|
|
for ( i = 1; i < 120; i++ ) {
|
|
raw->edge[i][0] = (raw->center[i-1][0] + raw->center[i][0]) / 2.0;
|
|
raw->edge[i][120] = (raw->center[i-1][119] + raw->center[i][119]) / 2.0;
|
|
raw->edge[0][i] = (raw->center[0][i-1] + raw->center[0][i]) / 2.0;
|
|
raw->edge[120][i] = (raw->center[119][i-1] + raw->center[119][i]) / 2.0;
|
|
}
|
|
|
|
/* derive internal nodes */
|
|
for ( j = 1; j < 120; j++ ) {
|
|
for ( i = 1; i < 120; i++ ) {
|
|
raw->edge[i][j] = ( raw->center[i-1][j-1] +
|
|
raw->center[i] [j-1] +
|
|
raw->center[i] [j] +
|
|
raw->center[i-1][j] ) / 4;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
/* Dump out the ascii format DEM file */
|
|
void rawDumpAsciiDEM( fgRAWDEM *raw, char *path, int ilon, int ilat ) {
|
|
char outfile[256];
|
|
char tmp[256];
|
|
int lon, lat;
|
|
char lon_sign, lat_sign;
|
|
int i, j;
|
|
FILE *fd;
|
|
|
|
/* Generate output file name */
|
|
|
|
if ( ilon >= 0 ) {
|
|
lon = ilon;
|
|
lon_sign = 'e';
|
|
} else {
|
|
lon = -ilon;
|
|
lon_sign = 'w';
|
|
}
|
|
|
|
if ( ilat >= 0 ) {
|
|
lat = ilat;
|
|
lat_sign = 'n';
|
|
} else {
|
|
lat = -ilat;
|
|
lat_sign = 's';
|
|
}
|
|
|
|
sprintf(outfile, "%s/%c%03d%c%03d.dem", path, lon_sign, lon, lat_sign, lat);
|
|
|
|
printf("outfile = %s\n", outfile);
|
|
|
|
if ( (fd = fopen(outfile, "w")) == NULL ) {
|
|
printf("Error opening output file = %s\n", outfile);
|
|
exit(-1);
|
|
}
|
|
|
|
/* Dump the "A" record */
|
|
|
|
/* print descriptive header (144 characters) */
|
|
sprintf(tmp, "%s - Generated from a 30 arcsec binary DEM", outfile);
|
|
fprintf(fd, "%-144s", tmp);
|
|
|
|
/* DEM level code, 3 reflects processing by DMA */
|
|
fprintf(fd, "%6d", 1);
|
|
|
|
/* Pattern code, 1 indicates a regular elevation pattern */
|
|
fprintf(fd, "%6d", 1);
|
|
|
|
/* Planimetric reference system code, 0 indicates geographic
|
|
* coordinate system. */
|
|
fprintf(fd, "%6d", 0);
|
|
|
|
/* Zone code */
|
|
fprintf(fd, "%6d", 0);
|
|
|
|
/* Map projection parameters (ignored) */
|
|
for ( i = 0; i < 15; i++ ) {
|
|
fprintf(fd, "%6.1f%18s", 0.0, "");
|
|
}
|
|
|
|
/* Units code, 3 represents arc-seconds as the unit of measure for
|
|
* ground planimetric coordinates throughout the file. */
|
|
fprintf(fd, "%6d", 3);
|
|
|
|
/* Units code; 2 represents meters as the unit of measure for
|
|
* elevation coordinates throughout the file. */
|
|
fprintf(fd, "%6d", 2);
|
|
|
|
/* Number (n) of sides in the polygon which defines the coverage of
|
|
* the DEM file (usually equal to 4). */
|
|
fprintf(fd, "%6d", 4);
|
|
|
|
/* Ground coordinates of bounding box in arc-seconds */
|
|
fprintf(fd, "%20.15fD+06", ilon * 3600.0 / 1000000.0);
|
|
fprintf(fd, "%20.15fD+06", ilat * 3600.0 / 1000000.0);
|
|
|
|
fprintf(fd, "%20.15fD+06", ilon * 3600.0 / 1000000.0);
|
|
fprintf(fd, "%20.15fD+06", (ilat+1) * 3600.0 / 1000000.0);
|
|
|
|
fprintf(fd, "%20.15fD+06", (ilon+1) * 3600.0 / 1000000.0);
|
|
fprintf(fd, "%20.15fD+06", (ilat+1) * 3600.0 / 1000000.0);
|
|
|
|
fprintf(fd, "%20.15fD+06", (ilon+1) * 3600.0 / 1000000.0);
|
|
fprintf(fd, "%20.15fD+06", (ilat) * 3600.0 / 1000000.0);
|
|
|
|
/* Minimum/maximum elevations in meters */
|
|
fprintf(fd, " %20.15E", (double)raw->tmp_min);
|
|
fprintf(fd, " %20.15E", (double)raw->tmp_max);
|
|
|
|
/* Counterclockwise angle from the primary axis of ground
|
|
* planimetric referenced to the primary axis of the DEM local
|
|
* reference system. */
|
|
fprintf(fd, "%6.1f", 0.0);
|
|
|
|
/* Accuracy code; 0 indicates that a record of accuracy does not
|
|
* exist and that no record type C will follow. */
|
|
fprintf(fd, "%24d", 0);
|
|
|
|
/* DEM spacial resolution. Usually (3,3) (3,6) or (3,9)
|
|
* depending on latitude */
|
|
fprintf(fd, "%12.6E", 30.0);
|
|
fprintf(fd, "%12.6E", 30.0);
|
|
|
|
/* accuracy code */
|
|
fprintf(fd, "%12.6E", 1.0);
|
|
|
|
/* dimension of arrays to follow (1)*/
|
|
fprintf(fd, "%6d", 1);
|
|
|
|
/* number of profiles */
|
|
fprintf(fd, "%6d", 3600 / raw->ydim + 1);
|
|
|
|
/* pad the end */
|
|
fprintf(fd, "%160s", "");
|
|
|
|
|
|
/* Dump "B" records */
|
|
|
|
for ( j = 0; j <= 120; j++ ) {
|
|
/* row / column id of this profile */
|
|
fprintf(fd, "%6d%6d", 1, j + 1);
|
|
|
|
/* Number of rows and columns (elevation points) in this
|
|
profile */
|
|
fprintf(fd, "%6d%6d", 3600 / raw->xdim + 1, 1);
|
|
|
|
/* Ground planimetric coordinates (arc-seconds) of the first
|
|
* elevation in the profile */
|
|
fprintf(fd, "%20.15fD+06", ilon * 3600.0 / 1000000.0);
|
|
fprintf(fd, "%20.15fD+06", (ilat * 3600.0 + j * raw->ydim) / 1000000.0);
|
|
|
|
/* Elevation of local datum for the profile. Always zero for
|
|
* 1-degree DEM, the reference is mean sea level. */
|
|
fprintf(fd, "%6.1f", 0.0);
|
|
fprintf(fd, "%18s", "");
|
|
|
|
/* Minimum and maximum elevations for the profile. */
|
|
fprintf(fd, " %20.15E", 0.0);
|
|
fprintf(fd, " %20.15E", 0.0);
|
|
|
|
/* One (usually) dimensional array (1,prof_num_cols) of
|
|
elevations */
|
|
for ( i = 0; i <= 120; i++ ) {
|
|
fprintf(fd, "%6.0f", raw->edge[j][i]);
|
|
}
|
|
}
|
|
|
|
fprintf(fd, "\n");
|
|
|
|
fclose(fd);
|
|
}
|
|
|
|
|
|
/* Read a horizontal strip of (1 vertical degree) from the raw DEM
|
|
* file specified by the upper latitude of the stripe specified in
|
|
* degrees. The output the individual ASCII format DEM tiles. */
|
|
void rawProcessStrip( fgRAWDEM *raw, int lat_degrees, char *path ) {
|
|
int lat, yrange;
|
|
int i, j, index, row, col;
|
|
int min, max;
|
|
int span, num_degrees, tile_width;
|
|
int xstart, xend;
|
|
|
|
/* convert to arcsec */
|
|
lat = lat_degrees * 3600;
|
|
|
|
printf("Max Latitude = %d arcsec\n", lat);
|
|
|
|
/* validity check ... */
|
|
if ( (lat > raw->rooty) ||
|
|
(lat < (raw->rooty - raw->nrows * raw->ydim + 1)) ) {
|
|
printf("Latitude out of range for this DEM file\n");
|
|
return;
|
|
}
|
|
|
|
printf ("Reading strip starting at %d (top and working down)\n", lat);
|
|
|
|
/* advance to the correct latitude */
|
|
rawAdvancePosition(raw, (raw->rooty - lat));
|
|
|
|
/* printf("short = %d\n", sizeof(short)); */
|
|
|
|
yrange = 3600 / raw->ydim;
|
|
|
|
for ( i = 0; i < yrange; i++ ) {
|
|
index = yrange - i - 1;
|
|
/* printf("About to read into row %d\n", index); */
|
|
rawReadNextRow(raw, index);
|
|
|
|
for ( j = 0; j < raw->ncols; j++ ) {
|
|
if ( raw->strip[index][j] == -9999 ) {
|
|
/* map ocean to 0 for now */
|
|
raw->strip[index][j] = 0;
|
|
}
|
|
}
|
|
}
|
|
|
|
/* extract individual tiles from the strip */
|
|
span = raw->ncols * raw->xdim;
|
|
num_degrees = span / 3600;
|
|
tile_width = raw->ncols / num_degrees;
|
|
printf("span = %d num_degrees = %d width = %d\n",
|
|
span, num_degrees, tile_width);
|
|
|
|
for ( i = 0; i < num_degrees; i++ ) {
|
|
xstart = i * tile_width;
|
|
xend = xstart + 120;
|
|
|
|
min = 10000; max = -10000;
|
|
for ( row = 0; row < yrange; row++ ) {
|
|
for ( col = xstart; col < xend; col++ ) {
|
|
/* Copy from strip to pixel centered tile. Yep,
|
|
* row/col are reversed here. raw->strip is backwards
|
|
* for convenience. I am converting to [x,y] now. */
|
|
raw->center[col-xstart][row] = raw->strip[row][col];
|
|
|
|
if ( raw->strip[row][col] < min) {
|
|
min = raw->strip[row][col];
|
|
}
|
|
|
|
if ( raw->strip[row][col] > max) {
|
|
max = raw->strip[row][col];
|
|
}
|
|
}
|
|
}
|
|
|
|
raw->tmp_min = min;
|
|
raw->tmp_max = max;
|
|
|
|
/* Convert from pixel centered to pixel edge values */
|
|
rawConvertCenter2Edge(raw);
|
|
|
|
/* Dump out the ascii format DEM file */
|
|
rawDumpAsciiDEM(raw, path, (raw->rootx / 3600) + i, lat_degrees - 1);
|
|
}
|
|
}
|
|
|
|
|
|
/* $Log$
|
|
/* Revision 1.6 1998/04/27 03:32:03 curt
|
|
/* Wrapped rint()'s in #ifdef HAVE_RINT
|
|
/*
|
|
* Revision 1.5 1998/04/18 03:59:46 curt
|
|
* Incorporated into gnu automake/autoconf system.
|
|
*
|
|
* Revision 1.4 1998/04/06 21:09:43 curt
|
|
* Additional win32 support.
|
|
* Fixed a bad bug in dem file parsing that was causing the output to be
|
|
* flipped about x = y.
|
|
*
|
|
* Revision 1.3 1998/03/03 13:10:29 curt
|
|
* Close to a working version.
|
|
*
|
|
* Revision 1.2 1998/03/03 02:04:01 curt
|
|
* Starting DEM Ascii format output routine.
|
|
*
|
|
* Revision 1.1 1998/03/02 23:31:01 curt
|
|
* Initial revision.
|
|
*
|
|
*/
|