--- /dev/null
+/************************************************************\
+ ** 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 **
+\************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <fcntl.h>
+#include <unistd.h>
+#include <bzlib.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;
+
+int ReadSRTM(char *filename)
+{
+ int x, y, infile, byte, bytes_read;
+ unsigned char error, buffer[2];
+ char north[3], west[4], *base=NULL;
+
+ if (strstr(filename, ".zip")!=NULL)
+ {
+ fprintf(stderr, "*** Error: \"%s\" must be uncompressed\n",filename);
+ return -1;
+
+ }
+
+ if (strstr(filename, ".hgt")==NULL)
+ {
+ fprintf(stderr, "*** Error: \"%s\" does not have the correct extension (.hgt)\n",filename);
+ return -1;
+ }
+
+ base=strrchr(filename, '/');
+
+ if (base==NULL)
+ base=filename;
+ else
+ base+=1;
+
+ 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;
+
+ 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;
+ }
+
+ 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;
+
+ infile=open(filename, O_RDONLY);
+
+ if (infile==0)
+ {
+ fprintf(stderr, "*** Error: Cannot open \"%s\"\n", filename);
+ return -1;
+ }
+
+ read(infile,&buffer,2);
+
+ if ((buffer[0]=='P') && (buffer[1]=='K'))
+ {
+ fprintf(stderr, "*** Error: \"%s\" still appears to be compressed!\n",filename);
+ close(infile);
+ return -1;
+ }
+
+ lseek(infile,0L,SEEK_SET);
+
+ 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++)
+ {
+ bytes_read=read(infile,&buffer,2);
+
+ if (bytes_read==2)
+ {
+ byte=buffer[1]+(buffer[0]<<8);
+
+ if (buffer[0]&128)
+ byte-=0x10000;
+
+ /* Flag problem elevations here */
+
+ if (byte<=min_elevation)
+ replacement_flag=1;
+
+ srtm[x][y]=byte;
+ }
+
+ else
+ error=1;
+ }
+
+ if (error)
+ {
+ fprintf(stderr,"\n*** Error: Premature EOF detected while reading \"%s\"! :-(\n",filename);
+ return -1;
+ }
+
+ close(infile);
+
+ return 0;
+}
+
+int LoadSDF_SDF(char *name)
+{
+ /* This function reads uncompressed
+ SPLAT Data Files (.sdf) into memory. */
+
+ int x, y, dummy;
+ char sdf_file[255], path_plus_name[255];
+ FILE *infile;
+
+ for (x=0; name[x]!='.' && name[x]!=0 && x<250; x++)
+ sdf_file[x]=name[x];
+
+ sdf_file[x]='.';
+ sdf_file[x+1]='s';
+ sdf_file[x+2]='d';
+ sdf_file[x+3]='f';
+ sdf_file[x+4]=0;
+
+ strncpy(path_plus_name,sdf_path,255);
+ strncat(path_plus_name,sdf_file,255);
+
+ infile=fopen(path_plus_name,"rb");
+
+ if (infile==NULL)
+ return 0;
+
+ fscanf(infile,"%d", &dummy);
+ fscanf(infile,"%d", &dummy);
+ fscanf(infile,"%d", &dummy);
+ 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]);
+
+ fclose(infile);
+
+ return 1;
+}
+
+char *BZfgets(BZFILE *bzfd, unsigned length)
+{
+ /* This function returns at most one less than 'length' number
+ of characters from a bz2 compressed file whose file descriptor
+ is pointed to by *bzfd. In operation, a buffer is filled with
+ uncompressed data (size = BZBUFFER), which is then parsed
+ and doled out as NULL terminated character strings every time
+ this function is invoked. A NULL string indicates an EOF
+ or error condition. */
+
+ static int x, y, nBuf;
+ static char buffer[BZBUFFER+1], output[BZBUFFER+1];
+ char done=0;
+
+ if (opened!=1 && bzerror==BZ_OK)
+ {
+ /* First time through. Initialize everything! */
+
+ x=0;
+ y=0;
+ nBuf=0;
+ opened=1;
+ output[0]=0;
+ }
+
+ do
+ {
+ if (x==nBuf && bzerror!=BZ_STREAM_END && bzerror==BZ_OK && opened)
+ {
+ /* Uncompress data into a static buffer */
+
+ nBuf=BZ2_bzRead(&bzerror, bzfd, buffer, BZBUFFER);
+ buffer[nBuf]=0;
+ x=0;
+ }
+
+ /* Build a string from buffer contents */
+
+ output[y]=buffer[x];
+
+ if (output[y]=='\n' || output[y]==0 || y==(int)length-1)
+ {
+ output[y+1]=0;
+ done=1;
+ y=0;
+ }
+
+ else
+ y++;
+ x++;
+
+ } while (done==0);
+
+ if (output[0]==0)
+ opened=0;
+
+ return (output);
+}
+
+int LoadSDF_BZ(char *name)
+{
+ /* This function reads .bz2 compressed
+ SPLAT Data Files into memory. */
+
+ int x, y, dummy;
+ char sdf_file[255], path_plus_name[255];
+ FILE *fd;
+ BZFILE *bzfd;
+
+ for (x=0; name[x]!='.' && name[x]!=0 && x<247; x++)
+ sdf_file[x]=name[x];
+
+ sdf_file[x]='.';
+ sdf_file[x+1]='s';
+ sdf_file[x+2]='d';
+ sdf_file[x+3]='f';
+ sdf_file[x+4]='.';
+ sdf_file[x+5]='b';
+ sdf_file[x+6]='z';
+ sdf_file[x+7]='2';
+ sdf_file[x+8]=0;
+
+ strncpy(path_plus_name,sdf_path,255);
+ strncat(path_plus_name,sdf_file,255);
+
+ fd=fopen(path_plus_name,"rb");
+ bzfd=BZ2_bzReadOpen(&bzerror,fd,0,0,NULL,0);
+
+ if (fd!=NULL && bzerror==BZ_OK)
+ {
+ printf("\nReading %s... ",path_plus_name);
+ fflush(stdout);
+
+ sscanf(BZfgets(bzfd,255),"%d",&dummy);
+ sscanf(BZfgets(bzfd,255),"%d",&dummy);
+ sscanf(BZfgets(bzfd,255),"%d",&dummy);
+ sscanf(BZfgets(bzfd,255),"%d",&dummy);
+
+ for (x=0; x<1200; x++)
+ for (y=0; y<1200; y++)
+ sscanf(BZfgets(bzfd,20),"%d",&usgs[x][y]);
+
+ fclose(fd);
+
+ BZ2_bzReadClose(&bzerror,bzfd);
+
+ return 1;
+ }
+
+ else
+ return 0;
+}
+
+char LoadSDF(char *name)
+{
+ /* This function loads the requested SDF file from the filesystem.
+ First, it tries to invoke the LoadSDF_SDF() function to load an
+ uncompressed SDF file (since uncompressed files load slightly
+ faster). Failing that, it tries to load a compressed SDF file
+ by invoking the LoadSDF_BZ() function. */
+
+ int return_value=-1;
+
+ /* Try to load an uncompressed SDF first. */
+
+ return_value=LoadSDF_SDF(name);
+
+ /* If that fails, try loading a compressed SDF. */
+
+ if (return_value==0 || return_value==-1)
+ return_value=LoadSDF_BZ(name);
+
+ return return_value;
+}
+
+int ReadUSGS()
+{
+ char usgs_filename[15];
+
+ /* usgs_filename is a minimal filename ("40:41:74:75").
+ Full path and extentions are added later though
+ subsequent function calls. */
+
+ sprintf(usgs_filename, "%d:%d:%d:%d", min_north, max_north, min_west, max_west);
+
+ return (LoadSDF(usgs_filename));
+}
+
+void average_terrain(x,y,z)
+int x, y, z;
+{
+ long accum;
+ int temp=0, count, bad_value;
+ double average;
+
+ bad_value=srtm[x][y];
+
+ accum=0L;
+ count=0;
+
+ if (x!=0)
+ {
+ temp=srtm[x-1][y];
+
+ if (temp>bad_value)
+ {
+ accum+=temp;
+ count++;
+ }
+ }
+
+ if (x!=1201)
+ {
+ temp=srtm[x+1][y];
+
+ if (temp>bad_value)
+ {
+ accum+=temp;
+ count++;
+ }
+ }
+
+ if ((x!=0) && (y!=1201))
+ {
+ temp=srtm[x-1][y+1];
+
+ if (temp>bad_value)
+ {
+ accum+=temp;
+ count++;
+ }
+ }
+
+ if (y!=1201)
+ {
+ temp=srtm[x][y+1];
+
+ if (temp>bad_value)
+ {
+ accum+=temp;
+ count++;
+ }
+ }
+
+ if ((x!=1201) && (y!=1201))
+ {
+ temp=srtm[x+1][y+1];
+
+ if (temp>bad_value)
+ {
+ accum+=temp;
+ count++;
+ }
+ }
+
+ if ((x!=0) && (y!=0))
+ {
+ temp=srtm[x-1][y-1];
+
+ if (temp>bad_value)
+ {
+ accum+=temp;
+ count++;
+ }
+ }
+
+ if (y!=0)
+ {
+ temp=srtm[x][y-1];
+
+ if (temp>bad_value)
+ {
+ accum+=temp;
+ count++;
+ }
+ }
+
+ if ((x!=1201) && (y!=0))
+ {
+ temp=srtm[x+1][y-1];
+
+ if (temp>bad_value)
+ {
+ accum+=temp;
+ count++;
+ }
+ }
+
+ if (count!=0)
+ {
+ average=(((double)accum)/((double)count));
+ temp=(int)(average+0.5);
+ }
+
+ if (temp>min_elevation)
+ srtm[x][y]=temp;
+ else
+ srtm[x][y]=min_elevation;
+}
+
+void WriteSDF(char *filename)
+{
+ int x, y, byte, last_good_byte=0;
+ FILE *outfile;
+
+ printf("\nWriting %s... ", filename);
+ fflush(stdout);
+
+ outfile=fopen(filename,"wb");
+
+ 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--)
+ {
+ byte=srtm[x][y];
+
+ if (byte>min_elevation)
+ last_good_byte=byte;
+
+ if (byte<min_elevation)
+ {
+ if (merge)
+ fprintf(outfile,"%d\n",usgs[1200-x][1200-y]);
+ else
+ {
+ average_terrain(x,y,last_good_byte);
+ fprintf(outfile,"%d\n",srtm[x][y]);
+ }
+ }
+ else
+ fprintf(outfile,"%d\n",byte);
+ }
+
+ printf("Done!\n");
+
+ fclose(outfile);
+}
+
+int main(int argc, char **argv)
+{
+ int x, y, z=0;
+ char *env=NULL, string[255];
+ FILE *fd;
+
+ 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");
+
+ 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");
+
+ return 1;
+ }
+
+ y=argc-1;
+
+ min_elevation=0;
+
+ for (x=1, z=0; x<=y; x++)
+ {
+ if (strcmp(argv[x],"-d")==0)
+ {
+ z=x+1;
+
+ if (z<=y && argv[z][0] && argv[z][0]!='-')
+ strncpy(sdf_path,argv[z],253);
+ }
+
+ if (strcmp(argv[x],"-n")==0)
+ {
+ z=x+1;
+
+ if (z<=y && argv[z][0])
+ {
+ sscanf(argv[z],"%d",&min_elevation);
+
+ if (min_elevation<-32768)
+ min_elevation=0;
+ }
+ }
+ }
+
+ /* If no SDF path was specified on the command line (-d), check
+ for a path specified in the $HOME/.splat_path file. If the
+ file is not found, then sdf_path[] remains NULL, and a data
+ merge will not be attempted if voids are found in the SRTM file. */
+
+ if (sdf_path[0]==0)
+ {
+ env=getenv("HOME");
+
+ sprintf(string,"%s/.splat_path",env);
+
+ fd=fopen(string,"r");
+
+ if (fd!=NULL)
+ {
+ fgets(string,253,fd);
+
+ /* Remove <CR> and/or <LF> from string */
+
+ for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0 && x<253; x++);
+ string[x]=0;
+
+ strncpy(sdf_path,string,253);
+
+ fclose(fd);
+ }
+ }
+
+ /* Ensure a trailing '/' is present in sdf_path */
+
+ if (sdf_path[0])
+ {
+ x=strlen(sdf_path);
+
+ if (sdf_path[x-1]!='/' && x!=0)
+ {
+ sdf_path[x]='/';
+ sdf_path[x+1]=0;
+ }
+ }
+
+ if (ReadSRTM(argv[z+1])==0)
+ {
+ if (replacement_flag && sdf_path[0])
+ merge=ReadUSGS();
+
+ WriteSDF(sdf_filename);
+ }
+
+ return 0;
+}
+
projects / debian / splat / blobdiff