/*****************************************************************************\ This module is public domain software that was developed by the U.S. Naval Oceanographic Office. This software 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. \*****************************************************************************/ /***************************************************************************\ * * * Module Name: dump_srtm_topo * * * * Programmer(s): Jan C. Depner * * * * Date Written: November 2006 * * * * Purpose: Reads a one-degree cell from the compressed SRTM * * files (*.cte) and dumps an ASCII YXZ file of the * * best available resolution data. * * * * Arguments: argv[1] - southwest latitude of cell * * argv[2] - southwest longitude of cell * * argv[3] - output file name * * * \***************************************************************************/ #include #include #include #include #include #include #include #include #include "nvtypes.h" #include "dump_srtm_topo_version.h" #include "read_srtm_topo.h" #include "get_area_mbr.h" #include "inside.h" #define NINT(a) ((a)<0.0L ? (NV_INT32) ((a) - 0.5L) : (NV_INT32) ((a) + 0.5L)) NV_INT32 main (NV_INT32 argc, NV_CHAR *argv[]) { NV_INT16 *array; NV_INT32 i, j, size, lat, lon, start_lat, start_lon, end_lat, end_lon, polygon_count, percent = 0, old_percent = -1; NV_FLOAT64 inc, pointlat, pointlon, polygon_x[200], polygon_y[200]; NV_F64_XYMBR mbr; FILE *fp; printf ("\n\n%s\n\n", VERSION); if (argc < 3) { fprintf (stderr, "Usage: %s AREA_FILE OUTPUT_FILE\n\n", argv[0]); fprintf (stderr, "\tWhere:\n\n"); fprintf (stderr, "\tAREA_FILE = Area file name\n\n"); fprintf (stderr, "\t\tThe area file name must have a .ARE extension\n"); fprintf (stderr, "\t\tfor ISS60 type area files. All others are assumed to be\n"); fprintf (stderr, "\t\tgeneric area files which consist of a simple list of\n"); fprintf (stderr, "\t\tpolygon points. The points may be in any of the following\n"); fprintf (stderr, "\t\tformats:\n\n"); fprintf (stderr, "\t\t\tHemisphere Degrees Minutes Seconds.decimal\n"); fprintf (stderr, "\t\t\tHemisphere Degrees Minutes.decimal\n"); fprintf (stderr, "\t\t\tHemisphere Degrees.decimal\n"); fprintf (stderr, "\t\t\tSign Degrees Minutes Seconds.decimal\n"); fprintf (stderr, "\t\t\tSign Degrees Minutes.decimal\n"); fprintf (stderr, "\t\t\tSign Degrees.decimal\n\n"); fprintf (stderr, "\t\tThe lat and lon must be entered one per line, separated by\n"); fprintf (stderr, "\t\ta comma. You do not need to repeat the first point, the\n"); fprintf (stderr, "\t\tpolygon will be closed automatically.\n\n\n"); fprintf (stderr, "\tOUTPUT_FILE = output file name\n\n"); exit (-1); } if ((fp = fopen (argv[2], "w")) == NULL) { perror (argv[2]); exit (-1); } get_area_mbr (argv[1], &polygon_count, polygon_x, polygon_y, &mbr); if (mbr.min_y < 0.0) { start_lat = (NV_INT32) (mbr.min_y - 1.0); } else { start_lat = (NV_INT32) mbr.min_y; } if (mbr.max_y < 0.0) { end_lat = (NV_INT32) mbr.max_y; } else { end_lat = (NV_INT32) (mbr.max_y + 1.0); } if (mbr.min_x < 0.0) { start_lon = (NV_INT32) (mbr.min_x - 1.0); } else { start_lon = (NV_INT32) mbr.min_x; } if (mbr.max_x < 0.0) { end_lon = (NV_INT32) mbr.max_x; } else { end_lon = (NV_INT32) (mbr.max_x + 1.0); } for (lat = start_lat ; lat < end_lat ; lat++) { for (lon = start_lon ; lon < end_lon ; lon++) { size = read_srtm_topo_one_degree (lat, lon, &array); if (size > 2) { inc = 1.0L / size; for (i = 0 ; i < size ; i++) { pointlat = ((NV_FLOAT64) lat + 1.0L) - (NV_FLOAT64) i * inc; for (j = 0 ; j < size ; j++) { pointlon = (NV_FLOAT64) lon + (NV_FLOAT64) j * inc; /* Make sure our point is inside our polygon. */ if (inside (polygon_x, polygon_y, polygon_count, pointlon, pointlat)) { /* Don't save zero values. This should be water but there's no way to tell in the SRTM data so we just say no ;-) */ if (array[i * size + j]) fprintf (fp, "%0.9f,%0.9f,%0.2f\n", pointlat, pointlon, (NV_FLOAT32) -array[i * size + j]); } } } } else { fprintf (stderr, "No data in cell at %d %d\n", lat, lon); } } percent = (NV_INT32) (((NV_FLOAT32) (lat - start_lat) / (NV_FLOAT32) (end_lat - start_lat)) * 100.0); if (percent != old_percent) { fprintf (stderr, "%03d%% processed\r", percent); old_percent = percent; } } fprintf (stderr, "100%% processed\n"); fclose (fp); cleanup_srtm_topo (); return (0); }