/*****************************************************************************\ 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. \*****************************************************************************/ /* 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. */ #include #include #include #include #include #include #include #include #include "nvtypes.h" #include "bit_pack.h" #include "huge_io.h" #define HEADER_SIZE 16384 #define MAP_BYTES (64800 * 36) / 8 #define NINT(a) ((a)<0.0 ? (NV_INT32) ((a) - 0.5) : (NV_INT32) ((a) + 0.5)) static NV_INT32 first = 1, prev_size = -1, fh = -1; static NV_INT16 *box = NULL; static NV_U_BYTE *map; /***************************************************************************\ * * * Module Name: read_srtm1_topo_one_degree * * * * Programmer(s): Jan C. Depner * * * * Date Written: October 2006 * * * * Purpose: Reads the SRTM compressed topographic elevation * * file (*.cte) and returns a one-degree single * * dimensioned array containing the elevations in the * * same format as the srtm1 data files. The * * width/height of the array is returned. * * * * Arguments: lat - degree of latitude, S negative * * lon - degree of longitude, W negative * * array - topo array * * * * Returns: 0 for all water cell, 2 for undefined cell, or * * 3600 (width/height of array, it's square). * * * * Caveats: The array is one dimensional so the user/caller * * must index into the array accordingly. The data is * * stored in order from the northwest corner of the * * cell, west to east, then north to south so the last * * point in the returned array is the southeast * * corner. See pointlat and pointlon in the following * * example code: * * * * * * #include "read_srtm1_topo.h" * * * * NV_INT16 *array; * * NV_INT32 size; * * NV_FLOAT64 inc, pointlat, pointlon; * * * * size = read_srtm1_topo_one_degree (lat, lon, * * &array); * * * * if (size > 2) * * { * * inc = 1.0L / size; * * for (i = 0 ; i < size ; i++) * * { * * pointlat = (lat + 1.0L) - i * inc; * * for (j = 0 ; j < size ; j++) * * { * * pointlon = lon + j * inc; * * //DO SOMETHING WITH array[i * size + j] * * } * * } * * } * * * * * * You should also call cleanup_srtm1_topo after you * * are finished using the database in order to close * * the open file and free the associated memory. * * * * * * This program reads a "huge" file. This is in * * actuality a directory containing one or more files * * with a nominal max size of 2GB. See huge_io.c for * * more detail. * * * * * \***************************************************************************/ NV_INT32 read_srtm1_topo_one_degree (NV_INT32 lat, NV_INT32 lon, NV_INT16 **array) { static NV_CHAR dir[512], file[512], version[128], created[128], zversion[128]; static NV_INT32 header_size, prev_lat = -999, prev_lon = -999; NV_CHAR varin[1024], info[1024], header_block[HEADER_SIZE]; NV_U_BYTE *buf, *bit_box = NULL, head[4]; NV_INT32 i, j, shift_lat, shift_lon, resolution, pos, size = 0, csize, bsize, status, ndx; NV_INT64 address; NV_INT16 start_val, bias, null_val, num_bits, temp, last_val; /* First time through, open the file and read the header. */ if (first) { if (getenv ("SRTM_DATA") == NULL) { fprintf (stderr, "Unable to find SRTM_DATA environment variable\n"); return (-1); } strcpy (dir, getenv ("SRTM_DATA")); sprintf (file, "%s/srtm1_topo.cte", dir); if ((fh = hfopen (file, "rb")) == -1) { perror (file); return (-1); } /* Allocate the map memory. */ map = (NV_U_BYTE *) calloc (MAP_BYTES, sizeof (NV_U_BYTE)); if (map == NULL) { perror ("Allocating map memory"); exit (-1); } /* Read the header block. */ hfread (header_block, HEADER_SIZE, 1, fh); ndx = 0; while (header_block[ndx] != 0) { for (i = 0 ; i < 1024 ; i++) { if (header_block[ndx] == '\n') break; varin[i] = header_block[ndx]; ndx++; } varin[i] = 0; if (strstr (varin, "[END OF HEADER]")) break; /* Put everything to the right of the equals sign in 'info'. */ if (strchr (varin, '=') != NULL) strcpy (info, (strchr (varin, '=') + 1)); if (strstr (varin, "[VERSION]")) strcpy (version, info); if (strstr (varin, "[ZLIB VERSION]")) strcpy (zversion, info); if (strstr (varin, "[CREATION DATE]")) strcpy (created, info); if (strstr (varin, "[HEADER SIZE]")) sscanf (info, "%d", &header_size); ndx++; } if (header_size != HEADER_SIZE) { fprintf (stderr, "Header sizes do not match, WTF, over!\n"); exit (-1); } /* Move past the end of the header and read the map. */ hfseek (fh, header_size, SEEK_SET); hfread (map, MAP_BYTES, 1, fh); first = 0; } /* Shift into a 0 to 180 by 0 to 360 world. */ shift_lat = lat + 90; shift_lon = lon + 180; if (lon > 360) lon -= 360; /* Only read the data if we have changed one-degree cells. */ if (prev_lat != shift_lat || prev_lon != shift_lon) { /* Unpack the address from the map. */ address = double_bit_unpack (map, (shift_lat * 360 + shift_lon) * 36, 36); prev_lat = shift_lat; prev_lon = shift_lon; /* If the address is 0 (water) or 2 (undefined), return the address. */ if (address < header_size) return ((NV_INT32) address); /* Move to the address and read/unpack the header. */ hfseek (fh, address, SEEK_SET); hfread (head, 8, 1, fh); pos = 0; resolution = (NV_INT32) bit_unpack (head, pos, 3); pos += 3; csize = (NV_INT32) bit_unpack (head, pos, 30); pos += 30; bsize = (NV_INT32) bit_unpack (head, pos, 31); size = 3600; /* We have to set an approximate size for unpacking (see the ZLIB documentation). */ bsize += NINT ((NV_FLOAT32) bsize * 0.10) + 12; /* Allocate the uncompressed memory. */ bit_box = (NV_U_BYTE *) calloc (bsize, sizeof (NV_U_BYTE)); if (bit_box == NULL) { perror ("Allocating bit_box memory in read_srtm1_topo"); exit (-1); } /* Allocate the compressed memory. */ buf = (NV_U_BYTE *) calloc (csize, sizeof (NV_U_BYTE)); if (buf == NULL) { perror ("Allocating buf memory"); exit (-1); } /* Read the compressed data. */ hfread (buf, csize, 1, fh); /* Uncompress the data. */ status = uncompress (bit_box, (uLongf *) &bsize, buf, csize); if (status) { fprintf (stderr, "Error %d uncompressing record\n", status); exit (-1); } free (buf); /* Unpack the internal header. */ pos = 0; start_val = bit_unpack (bit_box, pos, 16); pos += 16; bias = bit_unpack (bit_box, pos, 16); pos += 16; num_bits = bit_unpack (bit_box, pos, 4); pos += 4; null_val = NINT (pow (2.0L, (NV_FLOAT64) num_bits)) - 1; /* Allocate the cell memory. */ if (box != NULL) free (box); box = (NV_INT16 *) calloc (size * size, sizeof (NV_INT16)); if (box == NULL) { perror ("Allocating box memory in read_srtm1_topo"); exit (-1); } /* Uncompress the data (delta coded snake dance). */ last_val = start_val; for (i = 0 ; i < size ; i++) { if (!(i % 2)) { for (j = 0 ; j < size ; j++) { temp = bit_unpack (bit_box, pos, num_bits); pos += num_bits; if (temp < null_val) { box[i * size + j] = last_val + temp - bias; last_val = box[i * size + j]; } else { box[i * size + j] = -32768; } } } else { for (j = size - 1 ; j >= 0 ; j--) { temp = bit_unpack (bit_box, pos, num_bits); pos += num_bits; if (temp < null_val) { box[i * size + j] = last_val + temp - bias; last_val = box[i * size + j]; } else { box[i * size + j] = -32768; } } } } free (bit_box); prev_size = size; } /* We didn't change cells so set the same size as last time. */ else { size = prev_size; } *array = box; return (size); } /***************************************************************************\ * * * Module Name: read_srtm1_topo * * * * Programmer(s): Jan C. Depner * * * * Date Written: September 2006 * * * * Purpose: Reads the SRTM compressed topographic elevation * * file (*.cte) and returns the elevation value. If * * the value is undefined at that point it will return * * -32768. For water it will return 0. * * * * Arguments: lat - latitude degrees, S negative * * lon - longitude degrees, W negative * * * * Returns: 0 = water, -32768 undefined, or elevation * * * \***************************************************************************/ NV_INT16 read_srtm1_topo (NV_FLOAT64 lat, NV_FLOAT64 lon) { static NV_INT16 *array; static NV_INT32 prev_ilat = -1, prev_ilon = -1, size = 0; static NV_FLOAT64 inc = 0.0; NV_INT32 ilat, ilon, lat_index, lon_index; ilat = (NV_INT32) lat; ilon = (NV_INT32) lon; /* No point in calling the function if we didn't change cells. */ if (ilat != prev_ilat || ilon != prev_ilon) { size = read_srtm1_topo_one_degree (ilat, ilon, &array); if (size == 0) return (0); inc = 1.0L / (NV_FLOAT64) size; } if (size == 0) return (0); if (size == 2) return (-32768); /* Get the cell index. */ lat_index = size - (NV_INT32) ((lat - (NV_FLOAT64) ilat) / inc); lon_index = (NV_INT32) ((lon - (NV_FLOAT64) ilon) / inc); return (array[lat_index * size + lon_index]); } /***************************************************************************\ * * * Module Name: cleanup_srtm1_topo * * * * Programmer(s): Jan C. Depner * * * * Date Written: October 2006 * * * * Purpose: Closes the open srtm1 topo file, frees memory, and * * sets the first flag back to 1. * * * * Arguments: None * * * * Returns: Nada * * * \***************************************************************************/ void cleanup_srtm1_topo () { hfclose (fh); if (box != NULL) free (box); if (map != NULL) free (map); box = NULL; map = NULL; first = 1; prev_size = -1; }