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 ** SRTM-3 data in .BIL and .HGT format, and high **
6 ** resolution SRTM-1 one arc-second topography data. **
7 **************************************************************
8 ** Compile like this: **
9 ** cc -Wall -O3 -s -lbz2 srtm2sdf.c -o srtm2sdf **
10 ** Last modification: 12-Mar-2009 **
11 \**************************************************************/
20 #define BZBUFFER 65536
22 char sdf_filename[30], sdf_path[255], replacement_flag, opened=0,
25 int srtm[3601][3601], usgs[1201][1201], max_north, max_west, n,
26 min_north, min_west, merge=0, min_elevation, bzerror, ippd, mpi;
28 int ReadSRTM(char *filename)
30 int x, y, infile, byte=0, bytes_read;
31 unsigned char error, buffer[2];
32 char north[3], west[4], *base=NULL, blw_filename[255];
33 double cell_size, deg_north, deg_west;
36 if (strstr(filename, ".zip")!=NULL)
38 fprintf(stderr, "*** Error: \"%s\" must be uncompressed\n",filename);
43 if (strstr(filename, ".tgz")!=NULL)
45 fprintf(stderr, "*** Error: \"%s\" must be uncompressed\n",filename);
50 if ((strstr(filename, ".hgt")==NULL) && (strstr(filename, ".bil")==NULL))
52 fprintf(stderr, "*** Error: \"%s\" does not have the correct extension (.hgt or .bil)\n",filename);
56 if (strstr(filename, ".hgt")!=NULL)
59 if (strstr(filename, ".bil")!=NULL)
62 base=strrchr(filename, '/');
71 /* We obtain coordinates from the base of the .HGT filename */
82 if ((base[0]!='N' && base[0]!='S') || (base[3]!='W' && base[3]!='E'))
84 fprintf(stderr, "*** Error: \"%s\" doesn't look like a valid .hgt SRTM filename.\n", filename);
91 max_west=360-max_west;
99 min_north=atoi(north);
101 min_north=-atoi(north);
103 max_north=min_north+1;
108 /* We obtain .BIL file coordinates
109 from the corresponding .BLW file */
111 strncpy(blw_filename,filename,250);
116 blw_filename[x-2]='l';
117 blw_filename[x-1]='w';
120 fd=fopen(blw_filename,"rb");
124 n=fscanf(fd,"%lf",&cell_size);
126 if ((cell_size<0.0008) || (cell_size>0.0009))
128 printf("\n*** .BIL file's cell size is incompatible with SPLAT!!\n");
132 n=fscanf(fd,"%lf",°_west);
133 n=fscanf(fd,"%lf",°_west);
134 n=fscanf(fd,"%lf",°_west);
136 n=fscanf(fd,"%lf",°_west);
138 n=fscanf(fd,"%lf",°_north);
143 min_north=(int)(deg_north);
144 max_north=min_north+1;
149 deg_west=360.0-deg_west;
151 min_west=(int)(deg_west);
160 infile=open(filename, O_RDONLY);
164 fprintf(stderr, "*** Error: Cannot open \"%s\"\n", filename);
168 n=read(infile,&buffer,2);
170 if ((buffer[0]=='P') && (buffer[1]=='K'))
172 fprintf(stderr, "*** Error: \"%s\" still appears to be compressed!\n",filename);
177 lseek(infile,0L,SEEK_SET);
180 sprintf(sdf_filename, "%d:%d:%d:%d-hd.sdf", min_north, max_north, min_west, max_west);
182 sprintf(sdf_filename, "%d:%d:%d:%d.sdf", min_north, max_north, min_west, max_west);
187 printf("Reading %s... ", filename);
190 for (x=0; (x<=ippd && error==0); x++)
191 for (y=0; (y<=ippd && error==0); y++)
193 bytes_read=read(infile,&buffer,2);
199 /* "little-endian" structure */
201 byte=buffer[0]+(buffer[1]<<8);
209 /* "big-endian" structure */
211 byte=buffer[1]+(buffer[0]<<8);
217 /* Flag problem elevations here */
225 if (byte<=min_elevation)
237 fprintf(stderr,"\n*** Error: Premature EOF detected while reading \"%s\"! :-(\n",filename);
245 int LoadSDF_SDF(char *name)
247 /* This function reads uncompressed
248 SPLAT Data Files (.sdf) into memory. */
251 char sdf_file[255], path_plus_name[255];
254 for (x=0; name[x]!='.' && name[x]!=0 && x<250; x++)
263 strncpy(path_plus_name,sdf_path,255);
264 strncat(path_plus_name,sdf_file,255);
266 infile=fopen(path_plus_name,"rb");
271 n=fscanf(infile,"%d", &dummy);
272 n=fscanf(infile,"%d", &dummy);
273 n=fscanf(infile,"%d", &dummy);
274 n=fscanf(infile,"%d", &dummy);
276 printf("\nReading %s... ",path_plus_name);
279 for (x=0; x<1200; x++)
280 for (y=0; y<1200; y++)
281 n=fscanf(infile,"%d",&usgs[x][y]);
288 char *BZfgets(BZFILE *bzfd, unsigned length)
290 /* This function returns at most one less than 'length' number
291 of characters from a bz2 compressed file whose file descriptor
292 is pointed to by *bzfd. In operation, a buffer is filled with
293 uncompressed data (size = BZBUFFER), which is then parsed
294 and doled out as NULL terminated character strings every time
295 this function is invoked. A NULL string indicates an EOF
296 or error condition. */
298 static int x, y, nBuf;
299 static char buffer[BZBUFFER+1], output[BZBUFFER+1];
302 if (opened!=1 && bzerror==BZ_OK)
304 /* First time through. Initialize everything! */
315 if (x==nBuf && bzerror!=BZ_STREAM_END && bzerror==BZ_OK && opened)
317 /* Uncompress data into a static buffer */
319 nBuf=BZ2_bzRead(&bzerror, bzfd, buffer, BZBUFFER);
324 /* Build a string from buffer contents */
328 if (output[y]=='\n' || output[y]==0 || y==(int)length-1)
347 int LoadSDF_BZ(char *name)
349 /* This function reads .bz2 compressed
350 SPLAT Data Files into memory. */
353 char sdf_file[255], path_plus_name[255];
357 for (x=0; name[x]!='.' && name[x]!=0 && x<247; x++)
370 strncpy(path_plus_name,sdf_path,255);
371 strncat(path_plus_name,sdf_file,255);
373 fd=fopen(path_plus_name,"rb");
374 bzfd=BZ2_bzReadOpen(&bzerror,fd,0,0,NULL,0);
376 if (fd!=NULL && bzerror==BZ_OK)
378 printf("\nReading %s... ",path_plus_name);
381 sscanf(BZfgets(bzfd,255),"%d",&dummy);
382 sscanf(BZfgets(bzfd,255),"%d",&dummy);
383 sscanf(BZfgets(bzfd,255),"%d",&dummy);
384 sscanf(BZfgets(bzfd,255),"%d",&dummy);
386 for (x=0; x<1200; x++)
387 for (y=0; y<1200; y++)
388 sscanf(BZfgets(bzfd,20),"%d",&usgs[x][y]);
392 BZ2_bzReadClose(&bzerror,bzfd);
401 char LoadSDF(char *name)
403 /* This function loads the requested SDF file from the filesystem.
404 First, it tries to invoke the LoadSDF_SDF() function to load an
405 uncompressed SDF file (since uncompressed files load slightly
406 faster). Failing that, it tries to load a compressed SDF file
407 by invoking the LoadSDF_BZ() function. */
411 /* Try to load an uncompressed SDF first. */
413 return_value=LoadSDF_SDF(name);
415 /* If that fails, try loading a compressed SDF. */
417 if (return_value==0 || return_value==-1)
418 return_value=LoadSDF_BZ(name);
425 char usgs_filename[15];
427 /* usgs_filename is a minimal filename ("40:41:74:75").
428 Full path and extentions are added later though
429 subsequent function calls. */
431 sprintf(usgs_filename, "%d:%d:%d:%d", min_north, max_north, min_west, max_west);
433 return (LoadSDF(usgs_filename));
436 void average_terrain(y,x,z)
440 int temp=0, count, bad_value;
443 bad_value=srtm[y][x];
470 if ((y>=2) && (x<=(mpi-1)))
492 if ((x<=(mpi-1)) && (y<=mpi))
503 if ((x>=1) && (y>=2))
525 if ((y<=mpi) && (x>=1))
538 average=(((double)accum)/((double)count));
539 temp=(int)(average+0.5);
542 if (temp>min_elevation)
545 srtm[y][x]=min_elevation;
548 void WriteSDF(char *filename)
550 /* Like the HGT files, the extreme southwest corner
551 * provides the point of reference for the SDF file.
552 * The SDF file extends from min_north degrees to
553 * the south to min_north+(mpi/ippd) degrees to
554 * the north, and max_west degrees to the west to
555 * max_west-(mpi/ippd) degrees to the east. The
556 * overlapping edge redundancy present in the HGT
557 * and earlier USGS files is not necessary, nor
558 * is it present in SDF files.
561 int x, y, byte, last_good_byte=0;
564 printf("\nWriting %s... ", filename);
567 outfile=fopen(filename,"wb");
569 fprintf(outfile, "%d\n%d\n%d\n%d\n", max_west, min_north, min_west, max_north);
571 for (y=ippd; y>=1; y--) /* Omit the northern most edge */
572 for (x=mpi; x>=0; x--) /* Omit the eastern most edge */
576 if (byte>min_elevation)
579 if (byte<min_elevation)
584 fprintf(outfile,"%d\n",usgs[1200-(y/3)][1199-(x/3)]);
586 fprintf(outfile,"%d\n",usgs[1200-y][1199-x]);
591 average_terrain(y,x,last_good_byte);
592 fprintf(outfile,"%d\n",srtm[y][x]);
596 fprintf(outfile,"%d\n",byte);
604 int main(int argc, char **argv)
607 char *env=NULL, string[255], *s=NULL;
610 if (strstr(argv[0], "srtm2sdf-hd")!=NULL)
612 ippd=3600; /* High Definition (1 arc-sec) Mode */
613 strncpy(string,"srtm2sdf-hd\0",12);
618 ippd=1200; /* Standard Definition (3 arc-sec) Mode */
619 strncpy(string,"srtm2sdf\0",9);
622 mpi=ippd-1; /* Maximum pixel index per degree */
624 if (argc==1 || (argc==2 && strncmp(argv[1],"-h",2)==0))
627 fprintf(stderr, "\nsrtm2sdf: Generates SPLAT! elevation data files from unzipped\nSRTM-3 elevation data files, and replaces SRTM data voids with\nelevation data from older usgs2sdf derived SDF files.\n\n");
630 fprintf(stderr, "\nsrtm2sdf-hd: Generates SPLAT! elevation data files from unzipped\nSRTM-1 elevation data files, and replaces SRTM data voids with\naverages, or elevation data from older usgs2sdf derived SDF files.\n\n");
632 fprintf(stderr, "\tAvailable Options...\n\n");
633 fprintf(stderr, "\t-d directory path of usgs2sdf derived SDF files\n\t (overrides path in ~/.splat_path file)\n\n");
634 fprintf(stderr, "\t-n elevation limit (in meters) below which SRTM data is\n\t replaced by USGS-derived .sdf data (default = 0 meters)\n\n");
635 fprintf(stderr, "Examples: %s N40W074.hgt\n",string);
636 fprintf(stderr, " %s -d /cdrom/sdf N40W074.hgt\n",string);
637 fprintf(stderr, " %s -d /dev/null N40W074.hgt (prevents data replacement)\n",string);
638 fprintf(stderr, " %s -n -5 N40W074.hgt\n\n",string);
647 for (x=1, z=0; x<=y; x++)
649 if (strcmp(argv[x],"-d")==0)
653 if (z<=y && argv[z][0] && argv[z][0]!='-')
654 strncpy(sdf_path,argv[z],253);
657 if (strcmp(argv[x],"-n")==0)
661 if (z<=y && argv[z][0])
663 sscanf(argv[z],"%d",&min_elevation);
665 if (min_elevation<-32767)
671 /* If no SDF path was specified on the command line (-d), check
672 for a path specified in the $HOME/.splat_path file. If the
673 file is not found, then sdf_path[] remains NULL, and a data
674 merge will not be attempted if voids are found in the SRTM file. */
680 sprintf(string,"%s/.splat_path",env);
682 fd=fopen(string,"r");
686 s=fgets(string,253,fd);
688 /* Remove <CR> and/or <LF> from string */
690 for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0 && x<253; x++);
693 strncpy(sdf_path,string,253);
699 /* Ensure a trailing '/' is present in sdf_path */
705 if (sdf_path[x-1]!='/' && x!=0)
712 if (ReadSRTM(argv[z+1])==0)
714 if (replacement_flag && sdf_path[0])
717 WriteSDF(sdf_filename);