diff --git a/DemRaw2ascii/main.c b/DemRaw2ascii/main.c index 08ccd42ec..94aedb635 100644 --- a/DemRaw2ascii/main.c +++ b/DemRaw2ascii/main.c @@ -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. + * */ diff --git a/DemRaw2ascii/rawdem.c b/DemRaw2ascii/rawdem.c index 4004a9eed..88727d73d 100644 --- a/DemRaw2ascii/rawdem.c +++ b/DemRaw2ascii/rawdem.c @@ -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. * diff --git a/DemRaw2ascii/rawdem.h b/DemRaw2ascii/rawdem.h index 1bdb8e16c..8877b2eaa 100644 --- a/DemRaw2ascii/rawdem.h +++ b/DemRaw2ascii/rawdem.h @@ -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. + * */