-/************************************************************\
- ** 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 <stdio.h>
#include <stdlib.h>
#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)
{
}
- 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);
return -1;
}
- read(infile,&buffer,2);
+ n=read(infile,&buffer,2);
if ((buffer[0]=='P') && (buffer[1]=='K'))
{
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;
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;
if (error)
{
fprintf(stderr,"\n*** Error: Premature EOF detected while reading \"%s\"! :-(\n",filename);
- return -1;
}
close(infile);
/* 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;
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);
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)
{
}
}
- if (x!=1201)
+ if (y<=mpi)
{
- temp=srtm[x+1][y];
+ temp=srtm[y+1][x];
if (temp>bad_value)
{
}
}
- 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)
{
}
}
- if (y!=1201)
+ if (x<=(mpi-1))
{
- temp=srtm[x][y+1];
+ temp=srtm[y][x+1];
if (temp>bad_value)
{
}
}
- 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)
{
}
}
- 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)
{
}
}
- if (y!=0)
+ if (x>=1)
{
- temp=srtm[x][y-1];
+ temp=srtm[y][x-1];
if (temp>bad_value)
{
}
}
- 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)
{
}
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;
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;
if (byte<min_elevation)
{
if (merge)
- fprintf(outfile,"%d\n",usgs[1200-x][1200-y]);
+ {
+ if (ippd==3600)
+ fprintf(outfile,"%d\n",usgs[1200-(y/3)][1199-(x/3)]);
+ else
+ fprintf(outfile,"%d\n",usgs[1200-y][1199-x]);
+ }
+
else
- {
- average_terrain(x,y,last_good_byte);
- fprintf(outfile,"%d\n",srtm[x][y]);
- }
+ {
+ average_terrain(y,x,last_good_byte);
+ fprintf(outfile,"%d\n",srtm[y][x]);
+ }
}
else
fprintf(outfile,"%d\n",byte);
int main(int argc, char **argv)
{
int x, y, z=0;
- char *env=NULL, string[255];
+ char *env=NULL, string[255], *s=NULL;
FILE *fd;
+ if (strstr(argv[0], "srtm2sdf-hd")!=NULL)
+ {
+ ippd=3600; /* High Definition (1 arc-sec) Mode */
+ strncpy(string,"srtm2sdf-hd\0",12);
+ }
+
+ else
+ {
+ ippd=1200; /* Standard Definition (3 arc-sec) Mode */
+ strncpy(string,"srtm2sdf\0",9);
+ }
+
+ mpi=ippd-1; /* Maximum pixel index per degree */
+
if (argc==1 || (argc==2 && strncmp(argv[1],"-h",2)==0))
{
- 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");
- fprintf(stderr, "\tAvailable Options...\n\n");
- fprintf(stderr, "\t-d directory path of existing .sdf files\n\t (overrides path in ~/.splat_path file)\n\n");
- fprintf(stderr, "\t-n elevation (in meters) below which SRTM data\n\t is replaced by existing SDF data (default=0)\n\n");
+ if (ippd==1200)
+ 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");
+
+ if (ippd==3600)
+ 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");
- fprintf(stderr, "Examples: srtm2sdf N40W074.hgt\n");
- fprintf(stderr, " srtm2sdf -d /cdrom/sdf N40W074.hgt\n");
- fprintf(stderr, " srtm2sdf -d /dev/null N40W074.hgt (prevents data replacement)\n");
- fprintf(stderr, " srtm2sdf -n -5 N40W074.hgt\n\n");
+ fprintf(stderr, "\tAvailable Options...\n\n");
+ fprintf(stderr, "\t-d directory path of usgs2sdf derived SDF files\n\t (overrides path in ~/.splat_path file)\n\n");
+ 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");
+ fprintf(stderr, "Examples: %s N40W074.hgt\n",string);
+ fprintf(stderr, " %s -d /cdrom/sdf N40W074.hgt\n",string);
+ fprintf(stderr, " %s -d /dev/null N40W074.hgt (prevents data replacement)\n",string);
+ fprintf(stderr, " %s -n -5 N40W074.hgt\n\n",string);
return 1;
}
{
sscanf(argv[z],"%d",&min_elevation);
- if (min_elevation<-32768)
+ if (min_elevation<-32767)
min_elevation=0;
}
}
if (fd!=NULL)
{
- fgets(string,253,fd);
+ s=fgets(string,253,fd);
/* Remove <CR> and/or <LF> from string */