1 /**************************************************************\
2 ** Created originally by Jonathan Naylor, G4KLX **
3 ** Later embellished by John Magliacane, KD2BD to **
4 ** detect and handle voids found in the SRTM data **
5 ** and to handle SRTM data in .BIL as well as .HGT format. **
6 **************************************************************
7 ** Compile like this: **
8 ** cc -Wall -O3 -s -lbz2 srtm2sdf.c -o srtm2sdf **
9 ** Last modification: 23-Sep-2007 **
10 \**************************************************************/
19 #define BZBUFFER 65536
21 char sdf_filename[25], sdf_path[255], replacement_flag, opened=0,
24 int srtm[1201][1201], usgs[1200][1200], max_north, max_west,
25 min_north, min_west, merge=0, min_elevation, bzerror;
27 int ReadSRTM(char *filename)
29 int x, y, infile, byte=0, bytes_read;
30 unsigned char error, buffer[2];
31 char north[3], west[4], *base=NULL, blw_filename[255];
32 double cell_size, deg_north, deg_west;
35 if (strstr(filename, ".zip")!=NULL)
37 fprintf(stderr, "*** Error: \"%s\" must be uncompressed\n",filename);
42 if (strstr(filename, ".tgz")!=NULL)
44 fprintf(stderr, "*** Error: \"%s\" must be uncompressed\n",filename);
49 if ((strstr(filename, ".hgt")==NULL) && (strstr(filename, ".bil")==NULL))
51 fprintf(stderr, "*** Error: \"%s\" does not have the correct extension (.hgt or .bil)\n",filename);
55 if (strstr(filename, ".hgt")!=NULL)
58 if (strstr(filename, ".bil")!=NULL)
61 base=strrchr(filename, '/');
70 /* We obtain coordinates from the base of the .HGT filename */
81 if ((base[0]!='N' && base[0]!='S') || (base[3]!='W' && base[3]!='E'))
83 fprintf(stderr, "*** Error: \"%s\" doesn't look like a valid .hgt SRTM filename.\n", filename);
90 max_west=360-max_west;
98 min_north=atoi(north);
100 min_north=-atoi(north);
102 max_north=min_north+1;
107 /* We obtain .BIL file coordinates
108 from the corresponding .BLW file */
110 strncpy(blw_filename,filename,250);
115 blw_filename[x-2]='l';
116 blw_filename[x-1]='w';
119 fd=fopen(blw_filename,"rb");
123 fscanf(fd,"%lf",&cell_size);
125 if ((cell_size<0.0008) || (cell_size>0.0009))
127 printf("\n*** .BIL file's cell size is incompatible with SPLAT!!\n");
131 fscanf(fd,"%lf",°_west);
132 fscanf(fd,"%lf",°_west);
133 fscanf(fd,"%lf",°_west);
135 fscanf(fd,"%lf",°_west);
137 fscanf(fd,"%lf",°_north);
142 min_north=(int)(deg_north);
143 max_north=min_north+1;
148 deg_west=360.0-deg_west;
150 min_west=(int)(deg_west);
159 infile=open(filename, O_RDONLY);
163 fprintf(stderr, "*** Error: Cannot open \"%s\"\n", filename);
167 read(infile,&buffer,2);
169 if ((buffer[0]=='P') && (buffer[1]=='K'))
171 fprintf(stderr, "*** Error: \"%s\" still appears to be compressed!\n",filename);
176 lseek(infile,0L,SEEK_SET);
178 sprintf(sdf_filename, "%d:%d:%d:%d.sdf", min_north, max_north, min_west, max_west);
183 printf("Reading %s... ", filename);
186 for (x=0; (x<1201 && error==0); x++)
187 for (y=0; (y<1201 && error==0); y++)
189 bytes_read=read(infile,&buffer,2);
195 /* "little-endian" structure */
197 byte=buffer[0]+(buffer[1]<<8);
205 /* "big-endian" structure */
207 byte=buffer[1]+(buffer[0]<<8);
213 /* Flag problem elevations here */
221 if (byte<=min_elevation)
233 fprintf(stderr,"\n*** Error: Premature EOF detected while reading \"%s\"! :-(\n",filename);
241 int LoadSDF_SDF(char *name)
243 /* This function reads uncompressed
244 SPLAT Data Files (.sdf) into memory. */
247 char sdf_file[255], path_plus_name[255];
250 for (x=0; name[x]!='.' && name[x]!=0 && x<250; x++)
259 strncpy(path_plus_name,sdf_path,255);
260 strncat(path_plus_name,sdf_file,255);
262 infile=fopen(path_plus_name,"rb");
267 fscanf(infile,"%d", &dummy);
268 fscanf(infile,"%d", &dummy);
269 fscanf(infile,"%d", &dummy);
270 fscanf(infile,"%d", &dummy);
272 printf("\nReading %s... ",path_plus_name);
275 for (x=0; x<1200; x++)
276 for (y=0; y<1200; y++)
277 fscanf(infile,"%d",&usgs[x][y]);
284 char *BZfgets(BZFILE *bzfd, unsigned length)
286 /* This function returns at most one less than 'length' number
287 of characters from a bz2 compressed file whose file descriptor
288 is pointed to by *bzfd. In operation, a buffer is filled with
289 uncompressed data (size = BZBUFFER), which is then parsed
290 and doled out as NULL terminated character strings every time
291 this function is invoked. A NULL string indicates an EOF
292 or error condition. */
294 static int x, y, nBuf;
295 static char buffer[BZBUFFER+1], output[BZBUFFER+1];
298 if (opened!=1 && bzerror==BZ_OK)
300 /* First time through. Initialize everything! */
311 if (x==nBuf && bzerror!=BZ_STREAM_END && bzerror==BZ_OK && opened)
313 /* Uncompress data into a static buffer */
315 nBuf=BZ2_bzRead(&bzerror, bzfd, buffer, BZBUFFER);
320 /* Build a string from buffer contents */
324 if (output[y]=='\n' || output[y]==0 || y==(int)length-1)
343 int LoadSDF_BZ(char *name)
345 /* This function reads .bz2 compressed
346 SPLAT Data Files into memory. */
349 char sdf_file[255], path_plus_name[255];
353 for (x=0; name[x]!='.' && name[x]!=0 && x<247; x++)
366 strncpy(path_plus_name,sdf_path,255);
367 strncat(path_plus_name,sdf_file,255);
369 fd=fopen(path_plus_name,"rb");
370 bzfd=BZ2_bzReadOpen(&bzerror,fd,0,0,NULL,0);
372 if (fd!=NULL && bzerror==BZ_OK)
374 printf("\nReading %s... ",path_plus_name);
377 sscanf(BZfgets(bzfd,255),"%d",&dummy);
378 sscanf(BZfgets(bzfd,255),"%d",&dummy);
379 sscanf(BZfgets(bzfd,255),"%d",&dummy);
380 sscanf(BZfgets(bzfd,255),"%d",&dummy);
382 for (x=0; x<1200; x++)
383 for (y=0; y<1200; y++)
384 sscanf(BZfgets(bzfd,20),"%d",&usgs[x][y]);
388 BZ2_bzReadClose(&bzerror,bzfd);
397 char LoadSDF(char *name)
399 /* This function loads the requested SDF file from the filesystem.
400 First, it tries to invoke the LoadSDF_SDF() function to load an
401 uncompressed SDF file (since uncompressed files load slightly
402 faster). Failing that, it tries to load a compressed SDF file
403 by invoking the LoadSDF_BZ() function. */
407 /* Try to load an uncompressed SDF first. */
409 return_value=LoadSDF_SDF(name);
411 /* If that fails, try loading a compressed SDF. */
413 if (return_value==0 || return_value==-1)
414 return_value=LoadSDF_BZ(name);
421 char usgs_filename[15];
423 /* usgs_filename is a minimal filename ("40:41:74:75").
424 Full path and extentions are added later though
425 subsequent function calls. */
427 sprintf(usgs_filename, "%d:%d:%d:%d", min_north, max_north, min_west, max_west);
429 return (LoadSDF(usgs_filename));
432 void average_terrain(x,y,z)
436 int temp=0, count, bad_value;
439 bad_value=srtm[x][y];
466 if ((x!=0) && (y!=1201))
488 if ((x!=1201) && (y!=1201))
499 if ((x!=0) && (y!=0))
521 if ((x!=1201) && (y!=0))
534 average=(((double)accum)/((double)count));
535 temp=(int)(average+0.5);
538 if (temp>min_elevation)
541 srtm[x][y]=min_elevation;
544 void WriteSDF(char *filename)
546 int x, y, byte, last_good_byte=0;
549 printf("\nWriting %s... ", filename);
552 outfile=fopen(filename,"wb");
554 fprintf(outfile, "%d\n%d\n%d\n%d\n", max_west, min_north, min_west, max_north);
556 for (x=1200; x>0; x--)
557 for (y=1200; y>0; y--)
561 if (byte>min_elevation)
564 if (byte<min_elevation)
567 fprintf(outfile,"%d\n",usgs[1200-x][1200-y]);
570 average_terrain(x,y,last_good_byte);
571 fprintf(outfile,"%d\n",srtm[x][y]);
575 fprintf(outfile,"%d\n",byte);
583 int main(int argc, char **argv)
586 char *env=NULL, string[255];
589 if (argc==1 || (argc==2 && strncmp(argv[1],"-h",2)==0))
591 fprintf(stderr, "\nsrtm2sdf: Generates SPLAT! elevation data files from unzipped SRTM-3\nelevation data files, and replaces SRTM data voids with elevation\ndata from existing SDF files.\n\n");
592 fprintf(stderr, "\tAvailable Options...\n\n");
593 fprintf(stderr, "\t-d directory path of existing .sdf files\n\t (overrides path in ~/.splat_path file)\n\n");
594 fprintf(stderr, "\t-n elevation (in meters) below which SRTM data\n\t is replaced by existing SDF data (default=0)\n\n");
596 fprintf(stderr, "Examples: srtm2sdf N40W074.hgt\n");
597 fprintf(stderr, " srtm2sdf -d /cdrom/sdf N40W074.hgt\n");
598 fprintf(stderr, " srtm2sdf -d /dev/null N40W074.hgt (prevents data replacement)\n");
599 fprintf(stderr, " srtm2sdf -n -5 N40W074.hgt\n\n");
608 for (x=1, z=0; x<=y; x++)
610 if (strcmp(argv[x],"-d")==0)
614 if (z<=y && argv[z][0] && argv[z][0]!='-')
615 strncpy(sdf_path,argv[z],253);
618 if (strcmp(argv[x],"-n")==0)
622 if (z<=y && argv[z][0])
624 sscanf(argv[z],"%d",&min_elevation);
626 if (min_elevation<-32767)
632 /* If no SDF path was specified on the command line (-d), check
633 for a path specified in the $HOME/.splat_path file. If the
634 file is not found, then sdf_path[] remains NULL, and a data
635 merge will not be attempted if voids are found in the SRTM file. */
641 sprintf(string,"%s/.splat_path",env);
643 fd=fopen(string,"r");
647 fgets(string,253,fd);
649 /* Remove <CR> and/or <LF> from string */
651 for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0 && x<253; x++);
654 strncpy(sdf_path,string,253);
660 /* Ensure a trailing '/' is present in sdf_path */
666 if (sdf_path[x-1]!='/' && x!=0)
673 if (ReadSRTM(argv[z+1])==0)
675 if (replacement_flag && sdf_path[0])
678 WriteSDF(sdf_filename);