X-Git-Url: https://git.gag.com/?a=blobdiff_plain;f=utils%2Fsrtm2sdf.c;h=c955eb9ca15df2fce673804c735c0fecf907178b;hb=8ea1d9feaa6e6b0bf044123e08f8225db02b7dc3;hp=a12ee35534ccfdcdf3f3b3ec18b30272ef6a3816;hpb=cae76b32deb53ddbfb94b44de132a72435f56e88;p=debian%2Fsplat diff --git a/utils/srtm2sdf.c b/utils/srtm2sdf.c index a12ee35..c955eb9 100644 --- a/utils/srtm2sdf.c +++ b/utils/srtm2sdf.c @@ -1,12 +1,14 @@ -/************************************************************\ - ** Created originally by Jonathan Naylor, G4KLX ** - ** Later embellished by John Magliacane, KD2BD to ** - ** detect and handle voids found in the SRTM data ** - ************************************************************ - ** Compile like this: ** - ** cc -Wall -O3 -s -lbz2 srtm2sdf.c -o srtm2sdf ** - ** Last modification: 18-Mar-2006 ** -\************************************************************/ +/**************************************************************\ + ** Created originally by Jonathan Naylor, G4KLX. ** + ** Later embellished by John Magliacane, KD2BD to ** + ** detect and handle voids found in the SRTM data, ** + ** SRTM-3 data in .BIL and .HGT format, and high ** + ** resolution SRTM-1 one arc-second topography data. ** + ************************************************************** + ** Compile like this: ** + ** cc -Wall -O3 -s -lbz2 srtm2sdf.c -o srtm2sdf ** + ** Last modification: 12-Mar-2009 ** +\**************************************************************/ #include #include @@ -17,15 +19,19 @@ #define BZBUFFER 65536 -char sdf_filename[25], sdf_path[255], replacement_flag, opened=0; -int srtm[1201][1201], usgs[1200][1200], max_north, max_west, - min_north, min_west, merge=0, min_elevation, bzerror; +char sdf_filename[30], sdf_path[255], replacement_flag, opened=0, + hgt=0, bil=0; + +int srtm[3601][3601], usgs[1201][1201], max_north, max_west, n, + min_north, min_west, merge=0, min_elevation, bzerror, ippd, mpi; int ReadSRTM(char *filename) { - int x, y, infile, byte, bytes_read; + int x, y, infile, byte=0, bytes_read; unsigned char error, buffer[2]; - char north[3], west[4], *base=NULL; + char north[3], west[4], *base=NULL, blw_filename[255]; + double cell_size, deg_north, deg_west; + FILE *fd=NULL; if (strstr(filename, ".zip")!=NULL) { @@ -34,50 +40,122 @@ int ReadSRTM(char *filename) } - if (strstr(filename, ".hgt")==NULL) + if (strstr(filename, ".tgz")!=NULL) + { + fprintf(stderr, "*** Error: \"%s\" must be uncompressed\n",filename); + return -1; + + } + + if ((strstr(filename, ".hgt")==NULL) && (strstr(filename, ".bil")==NULL)) { - fprintf(stderr, "*** Error: \"%s\" does not have the correct extension (.hgt)\n",filename); + fprintf(stderr, "*** Error: \"%s\" does not have the correct extension (.hgt or .bil)\n",filename); return -1; } + if (strstr(filename, ".hgt")!=NULL) + hgt=1; + + if (strstr(filename, ".bil")!=NULL) + bil=1; + base=strrchr(filename, '/'); - + if (base==NULL) base=filename; else base+=1; - north[0]=base[1]; - north[1]=base[2]; - north[2]=0; + if (hgt) + { + /* We obtain coordinates from the base of the .HGT filename */ + + north[0]=base[1]; + north[1]=base[2]; + north[2]=0; - west[0]=base[4]; - west[1]=base[5]; - west[2]=base[6]; - west[3]=0; + west[0]=base[4]; + west[1]=base[5]; + west[2]=base[6]; + west[3]=0; - if ((base[0]!='N' && base[0]!='S') || (base[3]!='W' && base[3]!='E')) - { - fprintf(stderr, "*** Error: \"%s\" doesn't look like a valid SRTM filename.\n", filename); - return -1; + if ((base[0]!='N' && base[0]!='S') || (base[3]!='W' && base[3]!='E')) + { + fprintf(stderr, "*** Error: \"%s\" doesn't look like a valid .hgt SRTM filename.\n", filename); + return -1; + } + + max_west=atoi(west); + + if (base[3]=='E') + max_west=360-max_west; + + min_west=max_west-1; + + if (max_west==360) + max_west=0; + + if (base[0]=='N') + min_north=atoi(north); + else + min_north=-atoi(north); + + max_north=min_north+1; } - max_west=atoi(west); + if (bil) + { + /* We obtain .BIL file coordinates + from the corresponding .BLW file */ - if (base[3]=='E') - max_west=360-max_west; + strncpy(blw_filename,filename,250); + x=strlen(filename); - min_west=max_west-1; + if (x>3) + { + blw_filename[x-2]='l'; + blw_filename[x-1]='w'; + blw_filename[x]=0; - if (max_west==360) - max_west=0; + fd=fopen(blw_filename,"rb"); - if (base[0]=='N') - min_north=atoi(north); - else - min_north=-atoi(north); + if (fd!=NULL) + { + n=fscanf(fd,"%lf",&cell_size); - max_north=min_north+1; + if ((cell_size<0.0008) || (cell_size>0.0009)) + { + printf("\n*** .BIL file's cell size is incompatible with SPLAT!!\n"); + exit(1); + } + + n=fscanf(fd,"%lf",°_west); + n=fscanf(fd,"%lf",°_west); + n=fscanf(fd,"%lf",°_west); + + n=fscanf(fd,"%lf",°_west); + + n=fscanf(fd,"%lf",°_north); + + fclose(fd); + } + + min_north=(int)(deg_north); + max_north=min_north+1; + + if (deg_west<0.0) + deg_west=-deg_west; + else + deg_west=360.0-deg_west; + + min_west=(int)(deg_west); + + if (min_west==360) + min_west=0; + + max_west=min_west+1; + } + } infile=open(filename, O_RDONLY); @@ -87,7 +165,7 @@ int ReadSRTM(char *filename) return -1; } - read(infile,&buffer,2); + n=read(infile,&buffer,2); if ((buffer[0]=='P') && (buffer[1]=='K')) { @@ -98,7 +176,10 @@ int ReadSRTM(char *filename) lseek(infile,0L,SEEK_SET); - sprintf(sdf_filename, "%d:%d:%d:%d.sdf", min_north, max_north, min_west, max_west); + if (ippd==3600) + sprintf(sdf_filename, "%d:%d:%d:%d-hd.sdf", min_north, max_north, min_west, max_west); + else + sprintf(sdf_filename, "%d:%d:%d:%d.sdf", min_north, max_north, min_west, max_west); error=0; replacement_flag=0; @@ -106,20 +187,41 @@ int ReadSRTM(char *filename) printf("Reading %s... ", filename); fflush(stdout); - for (x=0; (x<1201 && error==0); x++) - for (y=0; (y<1201 && error==0); y++) + for (x=0; (x<=ippd && error==0); x++) + for (y=0; (y<=ippd && error==0); y++) { bytes_read=read(infile,&buffer,2); if (bytes_read==2) { - byte=buffer[1]+(buffer[0]<<8); + if (bil) + { + /* "little-endian" structure */ + + byte=buffer[0]+(buffer[1]<<8); + + if (buffer[1]&128) + byte-=0x10000; + } - if (buffer[0]&128) - byte-=0x10000; + if (hgt) + { + /* "big-endian" structure */ + + byte=buffer[1]+(buffer[0]<<8); + + if (buffer[0]&128) + byte-=0x10000; + } /* Flag problem elevations here */ + if (byte<-32768) + byte=-32768; + + if (byte>32767) + byte=32767; + if (byte<=min_elevation) replacement_flag=1; @@ -133,7 +235,6 @@ int ReadSRTM(char *filename) if (error) { fprintf(stderr,"\n*** Error: Premature EOF detected while reading \"%s\"! :-(\n",filename); - return -1; } close(infile); @@ -146,7 +247,7 @@ int LoadSDF_SDF(char *name) /* This function reads uncompressed SPLAT Data Files (.sdf) into memory. */ - int x, y, dummy; + int x, y, n, dummy; char sdf_file[255], path_plus_name[255]; FILE *infile; @@ -167,17 +268,17 @@ int LoadSDF_SDF(char *name) if (infile==NULL) return 0; - fscanf(infile,"%d", &dummy); - fscanf(infile,"%d", &dummy); - fscanf(infile,"%d", &dummy); - fscanf(infile,"%d", &dummy); + n=fscanf(infile,"%d", &dummy); + n=fscanf(infile,"%d", &dummy); + n=fscanf(infile,"%d", &dummy); + n=fscanf(infile,"%d", &dummy); printf("\nReading %s... ",path_plus_name); fflush(stdout); for (x=0; x<1200; x++) for (y=0; y<1200; y++) - fscanf(infile,"%d",&usgs[x][y]); + n=fscanf(infile,"%d",&usgs[x][y]); fclose(infile); @@ -332,21 +433,21 @@ int ReadUSGS() return (LoadSDF(usgs_filename)); } -void average_terrain(x,y,z) +void average_terrain(y,x,z) int x, y, z; { long accum; int temp=0, count, bad_value; double average; - bad_value=srtm[x][y]; + bad_value=srtm[y][x]; accum=0L; count=0; - if (x!=0) + if (y>=2) { - temp=srtm[x-1][y]; + temp=srtm[y-1][x]; if (temp>bad_value) { @@ -355,9 +456,9 @@ int x, y, z; } } - if (x!=1201) + if (y<=mpi) { - temp=srtm[x+1][y]; + temp=srtm[y+1][x]; if (temp>bad_value) { @@ -366,9 +467,9 @@ int x, y, z; } } - if ((x!=0) && (y!=1201)) + if ((y>=2) && (x<=(mpi-1))) { - temp=srtm[x-1][y+1]; + temp=srtm[y-1][x+1]; if (temp>bad_value) { @@ -377,9 +478,9 @@ int x, y, z; } } - if (y!=1201) + if (x<=(mpi-1)) { - temp=srtm[x][y+1]; + temp=srtm[y][x+1]; if (temp>bad_value) { @@ -388,9 +489,9 @@ int x, y, z; } } - if ((x!=1201) && (y!=1201)) + if ((x<=(mpi-1)) && (y<=mpi)) { - temp=srtm[x+1][y+1]; + temp=srtm[y+1][x+1]; if (temp>bad_value) { @@ -399,9 +500,9 @@ int x, y, z; } } - if ((x!=0) && (y!=0)) + if ((x>=1) && (y>=2)) { - temp=srtm[x-1][y-1]; + temp=srtm[y-1][x-1]; if (temp>bad_value) { @@ -410,9 +511,9 @@ int x, y, z; } } - if (y!=0) + if (x>=1) { - temp=srtm[x][y-1]; + temp=srtm[y][x-1]; if (temp>bad_value) { @@ -421,9 +522,9 @@ int x, y, z; } } - if ((x!=1201) && (y!=0)) + if ((y<=mpi) && (x>=1)) { - temp=srtm[x+1][y-1]; + temp=srtm[y+1][x-1]; if (temp>bad_value) { @@ -439,13 +540,24 @@ int x, y, z; } if (temp>min_elevation) - srtm[x][y]=temp; + srtm[y][x]=temp; else - srtm[x][y]=min_elevation; + srtm[y][x]=min_elevation; } void WriteSDF(char *filename) { + /* Like the HGT files, the extreme southwest corner + * provides the point of reference for the SDF file. + * The SDF file extends from min_north degrees to + * the south to min_north+(mpi/ippd) degrees to + * the north, and max_west degrees to the west to + * max_west-(mpi/ippd) degrees to the east. The + * overlapping edge redundancy present in the HGT + * and earlier USGS files is not necessary, nor + * is it present in SDF files. + */ + int x, y, byte, last_good_byte=0; FILE *outfile; @@ -456,10 +568,10 @@ void WriteSDF(char *filename) fprintf(outfile, "%d\n%d\n%d\n%d\n", max_west, min_north, min_west, max_north); - for (x=1200; x>0; x--) - for (y=1200; y>0; y--) + for (y=ippd; y>=1; y--) /* Omit the northern most edge */ + for (x=mpi; x>=0; x--) /* Omit the eastern most edge */ { - byte=srtm[x][y]; + byte=srtm[y][x]; if (byte>min_elevation) last_good_byte=byte; @@ -467,12 +579,18 @@ void WriteSDF(char *filename) if (byte and/or from string */