1
0
Fork 0

Close to a working version.

This commit is contained in:
curt 1998-03-03 13:10:28 +00:00
parent 64114fd6e7
commit 4b577c8bab
3 changed files with 138 additions and 113 deletions

View file

@ -63,7 +63,7 @@ int main(int argc, char **argv) {
rawOpenDemFile(&raw, dem_file);
for ( i = 41; i <= 90; i++ ) {
rawReadStrip(&raw, i);
rawProcessStrip(&raw, i, output_dir);
}
/* close the raw data file */
@ -74,7 +74,10 @@ int main(int argc, char **argv) {
/* $Log$
/* Revision 1.1 1998/03/02 23:31:01 curt
/* Initial revision.
/* Revision 1.2 1998/03/03 13:10:28 curt
/* Close to a working version.
/*
* Revision 1.1 1998/03/02 23:31:01 curt
* Initial revision.
*
*/

View file

@ -204,152 +204,163 @@ void rawConvertCenter2Edge( fgRAWDEM *raw ) {
/* Dump out the ascii format DEM file */
void rawDumpAsciiDEM( fgRAWDEM *raw ) {
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 */
/* get the name field (144 characters) */
for ( i = 0; i < 144; i++ ) {
name[i] = fgetc(fd);
}
name[i+1] = '\0';
/* clean off the whitespace at the end */
for ( i = strlen(name)-2; i > 0; i-- ) {
if ( !isspace(name[i]) ) {
i=0;
} else {
name[i] = '\0';
}
}
printf(" Quad name field: %s\n", name);
/* get quadrangle id (now part of previous section */
/* next_token(fd, dem_quadrangle); */
/* printf(" Quadrangle = %s\n", dem_quadrangle); */
/* 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 */
inum = next_int(fd);
printf(" DEM level code = %d\n", inum);
fprintf(fd, "%6d", 1);
/* Pattern code, 1 indicates a regular elevation pattern */
inum = next_int(fd);
printf(" Pattern code = %d\n", inum);
fprintf(fd, "%6d", 1);
/* Planimetric reference system code, 0 indicates geographic
* coordinate system. */
inum = next_int(fd);
printf(" Planimetric reference code = %d\n", inum);
fprintf(fd, "%6d", 0);
/* Zone code */
inum = next_int(fd);
printf(" Zone code = %d\n", inum);
fprintf(fd, "%6d", 0);
/* Map projection parameters (ignored) */
for ( i = 0; i < 15; i++ ) {
dnum = next_double(fd);
/* printf("%d: %f\n",i,dnum); */
fprintf(fd, "%6.1f%18s", 0.0, "");
}
/* Units code, 3 represents arc-seconds as the unit of measure for
/* Units code, 3 represents arc-seconds as the unit of measure for
* ground planimetric coordinates throughout the file. */
inum = next_int(fd);
if ( inum != 3 ) {
printf(" Unknown (X,Y) units code = %d!\n", inum);
exit(-1);
}
fprintf(fd, "%6d", 3);
/* Units code; 2 represents meters as the unit of measure for
* elevation coordinates throughout the file. */
inum = next_int(fd);
if ( inum != 2 ) {
printf(" Unknown (Z) units code = %d!\n", inum);
exit(-1);
}
fprintf(fd, "%6d", 2);
/* Number (n) of sides in the polygon which defines the coverage of
* the DEM file (usually equal to 4). */
inum = next_int(fd);
if ( inum != 4 ) {
printf(" Unknown polygon dimension = %d!\n", inum);
exit(-1);
}
fprintf(fd, "%6d", 4);
/* Ground coordinates of bounding box in arc-seconds */
dem_x1 = m->originx = next_exp(fd);
dem_y1 = m->originy = next_exp(fd);
printf(" Origin = (%.2f,%.2f)\n", m->originx, m->originy);
fprintf(fd, "%20.15fD+06", ilon * 3600.0 / 1000000.0);
fprintf(fd, "%20.15fD+06", ilat * 3600.0 / 1000000.0);
dem_x2 = next_exp(fd);
dem_y2 = next_exp(fd);
fprintf(fd, "%20.15fD+06", ilon * 3600.0 / 1000000.0);
fprintf(fd, "%20.15fD+06", (ilat+1) * 3600.0 / 1000000.0);
dem_x3 = next_exp(fd);
dem_y3 = next_exp(fd);
fprintf(fd, "%20.15fD+06", (ilon+1) * 3600.0 / 1000000.0);
fprintf(fd, "%20.15fD+06", (ilat+1) * 3600.0 / 1000000.0);
dem_x4 = next_exp(fd);
dem_y4 = next_exp(fd);
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 */
dem_z1 = next_exp(fd);
dem_z2 = next_exp(fd);
printf(" Elevation range %.4f %.4f\n", dem_z1, dem_z2);
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. */
next_token(fd, token);
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. */
/* DEM spacial resolution. Usually (3,3,1) (3,6,1) or (3,9,1)
fprintf(fd, "%24d", 0);
/* DEM spacial resolution. Usually (3,3) (3,6) or (3,9)
* depending on latitude */
/* I will eventually have to do something with this for data at
* higher latitudes */
next_token(fd, token);
printf(" accuracy & spacial resolution string = %s\n", token);
i = strlen(token);
printf(" length = %d\n", i);
fprintf(fd, "%12.6E", 30.0);
fprintf(fd, "%12.6E", 30.0);
ptr = token + i - 12;
printf(" last field = %s = %.2f\n", ptr, atof(ptr));
ptr[0] = '\0';
ptr = ptr - 12;
m->col_step = atof(ptr);
printf(" last field = %s = %.2f\n", ptr, m->row_step);
ptr[0] = '\0';
ptr = ptr - 12;
m->row_step = atof(ptr);
printf(" last field = %s = %.2f\n", ptr, m->col_step);
ptr[0] = '\0';
/* accuracy code = atod(token) */
inum = atoi(token);
printf(" Accuracy code = %d\n", inum);
printf(" column step = %.2f row step = %.2f\n",
m->col_step, m->row_step);
/* accuracy code */
fprintf(fd, "%12.6E", 1.0);
/* dimension of arrays to follow (1)*/
next_token(fd, token);
fprintf(fd, "%6d", 1);
/* number of profiles */
dem_num_profiles = m->cols = m->rows = next_int(fd);
printf(" Expecting %d profiles\n", dem_num_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 columns and rows (elevations) 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 (prof_num_cols,1) 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 */
void rawReadStrip( fgRAWDEM *raw, int lat_degrees ) {
* 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 = 0, max = 0;
int min, max;
int span, num_degrees, tile_width;
int xstart, xend;
@ -384,17 +395,8 @@ void rawReadStrip( fgRAWDEM *raw, int lat_degrees ) {
/* map ocean to 0 for now */
raw->strip[index][j] = 0;
}
if ( raw->strip[index][j] < min) {
min = raw->strip[index][j];
}
if ( raw->strip[index][j] > max) {
max = raw->strip[index][j];
}
}
}
printf("min = %d max = %d\n", min, max);
/* extract individual tiles from the strip */
span = raw->ncols * raw->xdim;
@ -407,28 +409,43 @@ void rawReadStrip( fgRAWDEM *raw, int lat_degrees ) {
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);
rawDumpAsciiDEM(raw, path, (raw->rootx / 3600) + i, lat_degrees - 1);
}
}
/* $Log$
/* Revision 1.2 1998/03/03 02:04:01 curt
/* Starting DEM Ascii format output routine.
/* 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.
*

View file

@ -41,6 +41,8 @@ typedef struct {
int rooty; /* Y coord of upper left *edge* of DEM region in degrees */
int xdim; /* X dimension of a pixel */
int ydim; /* Y dimension of a pixel */
int tmp_min; /* current 1x1 degree tile minimum */
int tmp_max; /* current 1x1 degree tile maximum */
/* file ptr */
int fd; /* Raw DEM file descriptor */
@ -66,15 +68,18 @@ void rawCloseDemFile( fgRAWDEM *raw );
/* Read a horizontal strip of (1 vertical degree) from the raw DEM
* file specified by the upper latitude of the stripe specified in
* degrees */
void rawReadStrip( fgRAWDEM *raw, int lat );
* degrees. The output the individual ASCII format DEM tiles. */
void rawProcessStrip( fgRAWDEM *raw, int lat_degrees, char *path );
#endif /* _RAWDEM_H */
/* $Log$
/* Revision 1.1 1998/03/02 23:31:02 curt
/* Initial revision.
/* Revision 1.2 1998/03/03 13:10:30 curt
/* Close to a working version.
/*
* Revision 1.1 1998/03/02 23:31:02 curt
* Initial revision.
*
*/