/*****************************************************************************\ 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: srtm3_mask * * * * Programmer(s): Jan C. Depner * * * * Date Written: September 2006 * * * * Purpose: Reads the uncompressed SRTM3 files (*.hgt) and * * creates world (or as much as is covered) land mask. * * * * Arguments: argv[1] - output file name * * * * Caveats: ALL of the SRTM3 hgt files must be in the directory * * that you are running from (CWD). * * * \***************************************************************************/ /* Description of the compressed land mask (.clm) 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 * 4 bytes, binary, 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 - 1's and 0's (woo hoo) 3 bits - resolution, 0 = one second mask, 1 = 3 second mask, 2 = 30 second mask, others TBD 29 bits - size of the zlib level 9 compressed block (SB) SB bytes - data The data is stored as a series of single bits for water (0) and land (1). Each bit represents a one second, three second, or thirty second cell in the block. The block is a one-degree square. It will be 3600 X 3600, 1200 X 1200, 120 X 120, or 60 X 60 depending on the resolution. It is ordered in the same fashion as the srtm3 data, that is, west to east starting in the northwest corner and moving southward. 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 "srtm3_mask_version.h" #define HEADER_SIZE 16384 #define NINT(a) ((a)<0.0 ? (NV_INT32) ((a) - 0.5) : (NV_INT32) ((a) + 0.5)) /* 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, char *argv[]) { NV_INT32 i, j, k, m, shift_lat, shift_lon, percent = 0, old_percent = -1, size,n, hit_land, hit_water, pos, swap = 1, prev_value; NV_INT16 row[1201]; FILE *fp, *ofp; NV_CHAR string[512], ofile[512], lathem, lonhem; NV_U_BYTE zero = 0, mapbuf[4], water[4], land[4], in_buf[180000], out_buf[200000]; 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); } /* Are we on a big endian machine? */ if (big_endian ()) swap = 0;; /* Set the water and land flags for the map area. */ bit_pack (water, 0, 32, 0); bit_pack (land, 0, 32, 1); /* Open the output file. */ strcpy (ofile, argv[1]); if (strcmp (&ofile[strlen (ofile) - 4], ".clm")) strcat (ofile, ".clm"); if ((ofp = fopen (ofile, "wb")) == NULL) { perror (ofile); exit (-1); } /* Write the (minimalist) ASCII header. */ t = time (&t); cur_tm = gmtime (&t); fprintf (ofp, "[HEADER SIZE] = %d\n", HEADER_SIZE); fprintf (ofp, "[VERSION] = %s\n", VERSION); fprintf (ofp, "[ZLIB VERSION] = %s\n", zlibVersion ()); fprintf (ofp, "[CREATION DATE] = %s", asctime (cur_tm)); fprintf (ofp, "[END OF HEADER]\n"); /* Zero out the remainder of the header. */ j = ftell (ofp); for (i = j ; i < HEADER_SIZE ; i++) fwrite (&zero, 1, 1, ofp); /* Set the default for all map addresses to 2 (undefined). */ for (i = -90 ; i < 90 ; i++) { for (j = -180 ; j < 180 ; j++) { bit_pack (mapbuf, 0, 32, 2); fwrite (mapbuf, 4, 1, ofp); } } /* Loop through the entire world, lat then lon. */ for (i = -90 ; i < 90 ; i++) { lathem = 'N'; if (i < 0) lathem = 'S'; shift_lat = i + 90; /* Skip North 60 files since they are mostly hosed. */ if (i != 60) { for (j = -180 ; j < 180 ; j++) { lonhem = 'E'; if (j < 0) lonhem = 'W'; 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) { hit_land = 0; hit_water = 0; prev_value = -1; pos = 0; /* Note that we're only going to 1200 not 1201 because we don't need the redundant data. */ for (k = 0 ; k < 1200 ; k++) { /* Read one row (all 1201). */ fread (row, 1201 * sizeof (NV_INT16), 1, fp); for (m = 0 ; m < 1200 ; m++) { /* If we're on a little endian system we need to swap the bytes. */ if (swap) { #ifdef __GNUC__ swab (&row[m], &row[m], 2); #else _swab ((NV_CHAR *) &row[m], (NV_CHAR *) &row[m], 2); #endif } /* From what I've seen, -32768 means land but no defined elevation. */ if (row[m]) { hit_land = 1; bit_pack (in_buf, pos, 1, 1); } else { hit_water = 1; bit_pack (in_buf, pos, 1, 0); } prev_value = row[m]; pos++; } } /* All land cell. */ if (hit_land && !hit_water) { fseek (ofp, HEADER_SIZE + (shift_lat * 360 + shift_lon) * 4, SEEK_SET); fwrite (land, 4, 1, ofp); } else { /* There are no SRTM3 cells with all zeros (water). */ if ((!hit_water && !hit_land) || (hit_water & !hit_land)) { fprintf (stderr, "Hit land: %d Hit water; %d WTF, over!\n", hit_land, hit_water); exit (-1); } /* Save the address of where we're going to write the compressed block. */ fseek (ofp, 0, SEEK_END); pos = ftell (ofp); /* Compress at maximum level. */ size = 200000; n = compress2 (out_buf, (uLongf *) &size, in_buf, 180000, 9); /* Write the type of data and size of the buffer to the file. */ bit_pack (mapbuf, 0, 3, 1); bit_pack (mapbuf, 3, 29, size); fwrite (mapbuf, 4, 1, ofp); /* Write the buffer to the file. */ fwrite (out_buf, size, 1, ofp); /* Write the address of the block to the map. */ fseek (ofp, HEADER_SIZE + (shift_lat * 360 + shift_lon) * 4, SEEK_SET); bit_pack (mapbuf, 0, 32, pos); fwrite (mapbuf, 4, 1, ofp); } fclose (fp); } /* We didn't have a file so it must be all water... */ else { /* Except outside the range of the SRTM3 data and the three missing files mentioned above. */ if (i >= -56 && i < 60 && (i != 59 || j != -166) && (i != 59 || j != -167) && (i != 59 || j != 170)) { fseek (ofp, HEADER_SIZE + (shift_lat * 360 + shift_lon) * 4, SEEK_SET); fwrite (water, 4, 1, ofp); } } } } percent = (NV_INT32) (((NV_FLOAT32) (i + 90) / 180.0) * 100.0); if (percent != old_percent) { fprintf (stderr, "%03d%% processed\r", percent); old_percent = percent; } } fclose (ofp); fprintf (stderr, "100%% processed \n\n"); return (0); }