Merge commit 'upstream/1.3.0'
[debian/splat] / utils / srtm2sdf.c
index 0d29b41da17d080ed8f77aad2c12c48267b36707..c955eb9ca15df2fce673804c735c0fecf907178b 100644 (file)
@@ -1,12 +1,13 @@
 /**************************************************************\
- **     Created originally by Jonathan Naylor, G4KLX         **
+ **     Created originally by Jonathan Naylor, G4KLX.        **
  **     Later embellished by John Magliacane, KD2BD to       **
- **     detect and handle voids found in the SRTM data       **
- **  and to handle SRTM data in .BIL as well as .HGT format. **
+ **     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: 23-Sep-2007              **
+ **              Last modification: 12-Mar-2009              **
 \**************************************************************/ 
 
 #include <stdio.h>
 
 #define BZBUFFER 65536
 
-char   sdf_filename[25], sdf_path[255], replacement_flag, opened=0,
+char   sdf_filename[30], sdf_path[255], replacement_flag, opened=0,
        hgt=0, bil=0;
 
-int    srtm[1201][1201], usgs[1200][1200], max_north, max_west,
-       min_north, min_west, merge=0, min_elevation, bzerror;
+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)
 {
@@ -120,7 +121,7 @@ int ReadSRTM(char *filename)
 
                        if (fd!=NULL)
                        {
-                               fscanf(fd,"%lf",&cell_size);
+                               n=fscanf(fd,"%lf",&cell_size);
 
                                if ((cell_size<0.0008) || (cell_size>0.0009))
                                {
@@ -128,13 +129,13 @@ int ReadSRTM(char *filename)
                                        exit(1);
                                }
 
-                               fscanf(fd,"%lf",&deg_west);
-                               fscanf(fd,"%lf",&deg_west);
-                               fscanf(fd,"%lf",&deg_west);
+                               n=fscanf(fd,"%lf",&deg_west);
+                               n=fscanf(fd,"%lf",&deg_west);
+                               n=fscanf(fd,"%lf",&deg_west);
 
-                               fscanf(fd,"%lf",&deg_west);
+                               n=fscanf(fd,"%lf",&deg_west);
 
-                               fscanf(fd,"%lf",&deg_north);
+                               n=fscanf(fd,"%lf",&deg_north);
 
                                fclose(fd);
                        }
@@ -164,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'))
        {
@@ -175,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;
@@ -183,8 +187,8 @@ 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);
 
@@ -243,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;
 
@@ -264,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);
 
@@ -429,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)
                {
@@ -452,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)
                {
@@ -463,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)
                {
@@ -474,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)
                {
@@ -485,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)
                {
@@ -496,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)
                {
@@ -507,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)
                {
@@ -518,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)
                {
@@ -536,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;
 
@@ -553,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;
@@ -564,12 +579,18 @@ void WriteSDF(char *filename)
                        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);
@@ -583,20 +604,38 @@ void WriteSDF(char *filename)
 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");
 
-               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");
+               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, "\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;
        }
@@ -644,7 +683,7 @@ int main(int argc, char **argv)
 
                if (fd!=NULL)
                {
-                       fgets(string,253,fd);
+                       s=fgets(string,253,fd);
 
                        /* Remove <CR> and/or <LF> from string */