/*****************************************************************************\ 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: gtopo30_mask * * * * Programmer(s): Jan C. Depner * * * * Date Written: September 2006 * * * * Purpose: Reads the uncompressed GTOPO30 and/or SRTM30 files * * (*.DEM) and creates a world (or as much as is * * covered) land mask. * * * * Arguments: argv[1] - output file name * * * * Caveats: ALL of the GTOPO30/SRTM30 DEM files must be in the * * directory that you are running from (CWD). This * * program should really only be run after building * * the original srtm compressed land mask (.clm) file * * in order to fill in the areas not covered by * * the srtm3 data (90S to 56S and 60N to 90N). * * * \***************************************************************************/ /* 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 "gtopo30_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, rows, cols, nodata, lat_offset, lon_offset, header_size = 0, addr, start_lon; NV_INT16 *row; NV_FLOAT64 xdim = 0.0, ydim = 0.0; FILE *fp, *ofp; NV_CHAR string[512], ofile[512], ltstring[20], lonhem, varin[1024], info[1024], version[128], zversion[128], created[128]; NV_U_BYTE zero = 0, mapbuf[4], water[4], land[4], in_buf[1800], out_buf[2000]; 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 the file doesn't exist we must create it (not a normal situation). */ if ((ofp = fopen (ofile, "rb+")) == NULL) { if ((ofp = fopen (ofile, "wb+")) == NULL) { perror (ofile); exit (-1); } /* Write the (minimalist) ASCII header. */ t = time (&t); cur_tm = gmtime (&t); header_size = HEADER_SIZE; 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); } } } /* If it exists, open it and read the header. */ else { while (fgets (varin, sizeof (varin), ofp)) { 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); } } /* Loop over the entire world, south to north. */ for (i = -90 ; i < 90 ; i++) { shift_lat = i + 90; /* Figure out what lat band we're using. */ if (i < -60) { sprintf (ltstring, "S60"); lat_offset = 30 - shift_lat - 1; } else if (i < -10) { sprintf (ltstring, "S10"); lat_offset = 80 - shift_lat - 1; } else if (i < 40) { sprintf (ltstring, "N40"); lat_offset = 130 - shift_lat - 1; } else { sprintf (ltstring, "N90"); lat_offset = 180 - shift_lat - 1; } /* Loop over the entire world, west to east. */ for (j = -180 ; j < 180 ; j++) { shift_lon = j + 180; /* Read the 64800 cell map to see if the cell is "undefined" (2). */ fseek (ofp, header_size + (shift_lat * 360 + shift_lon) * 4, SEEK_SET); fread (mapbuf, 4, 1, ofp); addr = bit_unpack (mapbuf, 0, 32); /* Only add data to "undefined" cells. */ if (addr == 2) { /* Figure out what lon band we're using. */ if (i < -60) { start_lon = shift_lon / 60; start_lon *= 60; } else { start_lon = shift_lon / 40; start_lon *= 40; } if (start_lon < 180) { lonhem = 'W'; } else { lonhem = 'E'; } lon_offset = (shift_lon - start_lon) * 120; /* Build the file name. */ sprintf (string, "%1c%03d%s.HDR", lonhem, abs (start_lon - 180), ltstring); /* If we can open the header file, read it. */ if ((fp = fopen (string, "rb")) != NULL) { while (fgets (string, sizeof (string), fp) != NULL) { if (strstr (string, "NROWS")) sscanf (string, "NROWS %d", &rows); if (strstr (string, "NCOLS")) sscanf (string, "NCOLS %d", &cols); if (strstr (string, "NODATA")) sscanf (string, "NODATA %d", &nodata); if (strstr (string, "XDIM")) sscanf (string, "XDIM %lf", &xdim); if (strstr (string, "YDIM")) sscanf (string, "YDIM %lf", &ydim); } fclose (fp); row = (NV_INT16 *) calloc (cols, sizeof (NV_INT16)); if (row == NULL) { perror ("Allocating row"); exit (-1); } /* We opened the .HDR file so we should be able to open the .DEM file. */ sprintf (string, "%1c%03d%s.DEM", lonhem, abs (start_lon - 180), ltstring); if ((fp = fopen (string, "rb")) == NULL) { perror (string); exit (-1); } hit_land = 0; hit_water = 0; prev_value = -1; pos = 0; /* Find the cell we need. */ fseek (fp, lat_offset * 120 * cols * sizeof (NV_INT16), SEEK_SET); /* Loop through the cell. */ for (k = 0 ; k < 120 ; k++) { /* Read a row. */ fread (row, cols * sizeof (NV_INT16), 1, fp); /* Loop through the row. */ for (m = lon_offset ; m < lon_offset + 120 ; m++) { /* If our system is little endian 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, "nodata" means water. */ if (row[m] != nodata) { 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++; } } free (row); /* 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); } /* All water cell. */ else if (!hit_land && hit_water) { fseek (ofp, HEADER_SIZE + (shift_lat * 360 + shift_lon) * 4, SEEK_SET); fwrite (water, 4, 1, ofp); } /* Mixed land and water. */ else { if (!hit_water && !hit_land) { fprintf (stderr, "WTF, over!\n"); exit (-1); } fseek (ofp, 0, SEEK_END); pos = ftell (ofp); /* Compress at maximum level. */ size = 2000; n = compress2 (out_buf, (uLongf *) &size, in_buf, 1800, 9); /* Write the type of data and size of the buffer to the file. */ bit_pack (mapbuf, 0, 3, 2); 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); } } } percent = (NV_INT32) (((float) (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); }