/*****************************************************************************\ 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: srtm1_topo * * * * Programmer(s): Jan C. Depner * * * * Date Written: October 2006 * * * * Purpose: Reads the uncompressed SRTM1 files (*.hgt) and * * creates world (or as much as is covered) * * compressed topoographic elevation (.cte) file. * * * * Arguments: argv[1] - output file name * * * * Caveats: ALL of the SRTM1 hgt files must be in the directory * * that you are running from (CWD). This program * * creates a "huge" file. This is in actuality a * * directory containing a number of 2GB files. See * * huge_io.c for more detail. * * * \***************************************************************************/ /* Description of the compressed topographic elevation (.cte) file format (look Ma, no endians!) Header - 16384 bytes, ASCII [HEADER SIZE] = 16384 [CREATION DATE] = [VERSION] = [ZLIB VERSION] = [END OF HEADER] One-degree map - 64800 * 36 bits, double precision integer, stored as characters. Records start at 90S,180W and proceed west to east then south to north (that is, the second record is for 90S,179W and the 361st record is for 89S,180W). Record contains 0 if all water, 1 if all land, 2 if undefined, or address if both land and water. Data - 3 bits - Not used (originally planned for resolution in mixed files) 30 bits - size of the zlib level 9 compressed block (SB) 31 bits - size of the uncompressed block SB bytes - data Inside the compressed block the data has already been delta coded and bit packed. The format of the data stored in the compressed block is as follows: 16 bits - signed short, starting value, stored as characters 16 bits - signed short, bias value, stored as characters 4 bits - number of bits used to store delta coding offsets (NB) NB bits - first offset NB bits - offset from first NB bits - offset from second . . . NB bits - last offset (may be the 12,960,000th, 1,440,000th, or 14,400th offset depending on resolution) Undefined values (-32768) will be stored as (int) pow (2.0L, NB) - 1 and will not be used in the delta chain. The deltas are computed by subtracting the previous valid value from the current value. The data is traversed west to east for one row, then east to west for the next row, then west to east, etc. I call this the snake dance (named after the lines at Disney World ;-) After the delta is computed the bias is added to it. The compression is compliments of the ZLIB compression library which can be found at http://www.zlib.net/. Many thanks to Jean-loup Gailly, Mark Adler, and all others associated with that effort. For some unknown reason the majority of the N60 files in srtm3 are totally hosed up so I'm going to mark them as undefined (2). Also, N59E170, N59W167, and N59W166 seem to be missing so these are forced to undefined as well. */ #include #include #include #include #include #include #include #include #include "nvtypes.h" #include "bit_pack.h" #include "huge_io.h" #include "srtm1_topo_version.h" #define HEADER_SIZE 16384 #define NINT(a) ((a)<0.0L ? (NV_INT32) ((a) - 0.5L) : (NV_INT32) ((a) + 0.5L)) #define LOG2 0.3010299956639812L #define MAP_BYTES (64800 * 36) / 8 /* A function to determine if your system is big or litle endian. Returns a 0 if it's little endian or a 3 if it's big endian. */ static NV_INT32 big_endian () { union { NV_INT32 word; NV_U_BYTE byte[4]; } a; a.word = 0x00010203; return ((NV_INT32) a.byte[3]); } NV_INT32 main (NV_INT32 argc, NV_CHAR *argv[]) { NV_INT32 i, j, k, m, shift_lat, shift_lon, percent = 0, old_percent = -1, n, hit_land, hit_water, pos, swap = 1, prev_value, first, ndx, total_bits, total_bytes, out_bytes, ofh; FILE *fp; NV_INT16 **row, *delta, bias, start_val = 0, last_val = 0, diff, min_diff, max_diff, num_bits, null_val; NV_CHAR string[512], ofile[512], lathem, lonhem, header_block[HEADER_SIZE]; NV_INT64 lpos; NV_U_BYTE water[4], head[8], *in_buf, *out_buf, map[MAP_BYTES]; time_t t; struct tm *cur_tm; printf ("\n\n%s\n\n", VERSION); if (argc < 2) { fprintf (stderr, "Usage: %s OUTPUT_FILENAME\n\n", argv[0]); exit (-1); } /* Hmmm, had problems with statically allocated arrays of this size. */ delta = (NV_INT16 *) calloc (3600 * 3600, sizeof (NV_INT16)); if (delta == NULL) { perror ("Allocating delta"); exit (-1); } row = (NV_INT16 **) calloc (3600, sizeof (NV_INT16 *)); if (row == NULL) { perror ("Allocating row"); exit (-1); } for (i = 0 ; i < 3600 ; i++) { row[i] = (NV_INT16 *) calloc (3601, sizeof (NV_INT16)); if (row[i] == NULL) { perror ("Allocating row[i]"); exit (-1); } } /* Are we on a big endian machine? */ if (big_endian ()) swap = 0;; /* Set the water flag for the map area. */ bit_pack (water, 0, 32, 0); /* Open the output file. */ strcpy (ofile, argv[1]); if (strcmp (&ofile[strlen (ofile) - 4], ".cte")) strcat (ofile, ".cte"); if ((ofh = hfopen (ofile, "wb")) == -1) { perror (ofile); exit (-1); } /* Write the (minimalist) ASCII header. */ memset (header_block, 0, HEADER_SIZE); t = time (&t); cur_tm = gmtime (&t); sprintf (header_block, "[HEADER SIZE] = %d\n", HEADER_SIZE); sprintf (&header_block[strlen (header_block)], "[VERSION] = %s\n", VERSION); sprintf (&header_block[strlen (header_block)], "[ZLIB VERSION] = %s\n", zlibVersion ()); sprintf (&header_block[strlen (header_block)], "[CREATION DATE] = %s", asctime (cur_tm)); sprintf (&header_block[strlen (header_block)], "[END OF HEADER]\n"); hfwrite (header_block, HEADER_SIZE, 1, ofh); /* Set the default for all map addresses to 2 (undefined). Also, write the blank map area out to the file to save the space. */ for (i = 0 ; i < 180 ; i++) { for (j = 0 ; j < 360 ; j++) { double_bit_pack (map, (i * 360 + j) * 36, 36, 2); } } hfseek (ofh, HEADER_SIZE, SEEK_SET); hfwrite (map, MAP_BYTES, 1, ofh); /* Loop through -56 to 60 lat. */ for (i = -56 ; i < 60 ; i++) { lathem = 'N'; if (i < 0) lathem = 'S'; /* Shift into the 0 to 180 world. */ shift_lat = i + 90; /* Loop through the entire range of longitudes. */ for (j = -180 ; j < 180 ; j++) { lonhem = 'E'; if (j < 0) lonhem = 'W'; /* Shift into the 0 to 360 world. */ shift_lon = j + 180; /* Build the file name. */ sprintf (string, "%1c%02d%1c%03d.hgt", lathem, abs (i), lonhem, abs (j)); /* If we can open the file, read it. */ if ((fp = fopen (string, "rb")) != NULL) { prev_value = -1; pos = 0; /* Note that we're only going to 3600 not 3601 because we don't need the redundant data. */ for (k = 0 ; k < 3600 ; k++) { /* Read one row (all 3601). */ fread (row[k], 3601 * sizeof (NV_INT16), 1, fp); for (m = 0 ; m < 3600 ; m++) { /* If we're on a little endian system we need to swap the bytes. */ if (swap) { #ifdef __GNUC__ swab (&row[k][m], &row[k][m], 2); #else _swab ((NV_CHAR *) &row[k][m], (NV_CHAR *) &row[k][m], 2); #endif } } } first = 1; ndx = 0; min_diff = 32767; max_diff = -32768; hit_land = 0; hit_water = 0; /* Loop through the cell. */ for (k = 0 ; k < 3600 ; k++) { /* Snake dance test. */ if (!(k % 2)) { /* West to east. */ for (m = 0 ; m < 3600 ; m++) { if (row[k][m] != -32768) { /* Set the "hit" flags. */ if (row[k][m]) { hit_land = 1; } else { hit_water = 1; } /* First time through set the last_val to start_val. */ if (first) { start_val = row[k][m]; last_val = start_val; first = 0; } diff = row[k][m] - last_val; if (diff > max_diff) max_diff = diff; if (diff < min_diff) min_diff = diff; last_val = row[k][m]; delta[ndx++] = diff; } else { /* From what I've seen, -32768 means land but no defined elevation. */ hit_land = 1; delta[ndx++] = -32768; } } } else { /* East to west. */ for (m = 3599 ; m >= 0 ; m--) { if (row[k][m] != -32768) { /* Set the "hit" flags. */ if (row[k][m]) { hit_land = 1; } else { hit_water = 1; } /* First time through set the last_val to start_val. */ if (first) { start_val = row[k][m]; last_val = start_val; first = 0; } diff = row[k][m] - last_val; if (diff > max_diff) max_diff = diff; if (diff < min_diff) min_diff = diff; last_val = row[k][m]; delta[ndx++] = diff; } else { /* From what I've seen, -32768 means land but no defined elevation. */ hit_land = 1; delta[ndx++] = -32768; } } } } /* All water cell. */ if (!hit_land && hit_water) { double_bit_pack (map, (shift_lat * 360 + shift_lon) * 36, 36, 0); } /* Land and water or all land cell. */ else { /* The bias is the negative of the minimum difference. We bias so we don't have to play with sign extension in the bit unpacking. */ bias = -min_diff; /* Add two to the max difference to allow room for the null value. */ max_diff += (bias + 2); /* Compute the number of bits needed to store a delta. */ num_bits = NINT (log10 ((NV_FLOAT64) max_diff) / LOG2 + 0.5L); null_val = NINT (pow (2.0L, (NV_FLOAT64) num_bits)) - 1; /* Compute the total bytes needed to store the block. */ total_bits = 16 + 16 + 4 + 3600 * 3600 * num_bits; total_bytes = total_bits / 8; if (total_bits % 8) total_bytes++; /* Allocate the uncompressed memory block. */ in_buf = (NV_U_BYTE *) calloc (total_bytes, 1); if (in_buf == NULL) { perror ("Allocating in_buf"); exit (-1); } /* Allocate the compressed memory block. */ out_bytes = total_bytes + NINT ((NV_FLOAT32) total_bytes * 0.10) + 12; out_buf = (NV_U_BYTE *) calloc (out_bytes, 1); if (out_buf == NULL) { perror ("Allocating out_buf"); exit (-1); } /* Pack the internal header. */ pos = 0; bit_pack (in_buf, pos, 16, start_val); pos += 16; bit_pack (in_buf, pos, 16, bias); pos += 16; bit_pack (in_buf, pos, 4, num_bits); pos += 4; /* Pack the deltas. */ for (k = 0 ; k < 3600 * 3600 ; k++) { if (delta[k] == -32768) { delta[k] = null_val; } else { delta[k] += bias; } bit_pack (in_buf, pos, num_bits, delta[k]); pos += num_bits; } /* Compress at maximum level. */ n = compress2 (out_buf, (uLongf *) &out_bytes, in_buf, total_bytes, 9); if (n) { fprintf (stderr, "Error %d compressing record\n", n); exit (-1); } free (in_buf); /* Pack the header. */ pos = 0; bit_pack (head, pos, 3, 0); pos += 3; bit_pack (head, pos, 30, out_bytes); pos += 30; bit_pack (head, pos, 31, total_bytes); /* Get the address where we're going to write the compressed block. */ hfseek (ofh, 0, SEEK_END); lpos = hftell (ofh); /* Write the header to the file. */ hfwrite (head, 8, 1, ofh); /* Write the buffer to the file. */ hfwrite (out_buf, out_bytes, 1, ofh); free (out_buf); /* Save the address of the block to the map. */ double_bit_pack (map, (shift_lat * 360 + shift_lon) * 36, 36, lpos); } fclose (fp); } } percent = (NV_INT32) (((NV_FLOAT32) (i + 56) / 116.0) * 100.0); if (percent != old_percent) { fprintf (stderr, "%03d%% processed\r", percent); old_percent = percent; } } /* Write the map. */ hfseek (ofh, HEADER_SIZE, SEEK_SET); hfwrite (map, MAP_BYTES, 1, ofh); hfclose (ofh); fprintf (stderr, "100%% processed \n\n"); return (0); }