git.gag.com Git - debian/splat/blob - utils/srtm2sdf.c

   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 \**************************************************************/
  12
  13 #include <stdio.h>
  14 #include <stdlib.h>
  15 #include <string.h>
  16 #include <fcntl.h>
  17 #include <unistd.h>
  18 #include <bzlib.h>
  19
  20 #define BZBUFFER 65536
  21
  22 char    sdf_filename[30], sdf_path[255], replacement_flag, opened=0,
  23         hgt=0, bil=0;
  24
  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;
  27
  28 int ReadSRTM(char *filename)
  29 {
  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;
  34         FILE *fd=NULL;
  35
  36         if (strstr(filename, ".zip")!=NULL)
  37         {
  38                 fprintf(stderr, "*** Error: \"%s\" must be uncompressed\n",filename);
  39                 return -1;
  40
  41         }
  42
  43         if (strstr(filename, ".tgz")!=NULL)
  44         {
  45                 fprintf(stderr, "*** Error: \"%s\" must be uncompressed\n",filename);
  46                 return -1;
  47
  48         }
  49
  50         if ((strstr(filename, ".hgt")==NULL) && (strstr(filename, ".bil")==NULL))
  51         {
  52                 fprintf(stderr, "*** Error: \"%s\" does not have the correct extension (.hgt or .bil)\n",filename);
  53                 return -1;
  54         }
  55
  56         if (strstr(filename, ".hgt")!=NULL)
  57                 hgt=1;
  58
  59         if (strstr(filename, ".bil")!=NULL)
  60                 bil=1;
  61
  62         base=strrchr(filename, '/');
  63
  64         if (base==NULL)
  65                 base=filename;
  66         else
  67                 base+=1;
  68
  69         if (hgt)
  70         {
  71                 /* We obtain coordinates from the base of the .HGT filename */
  72
  73                 north[0]=base[1];
  74                 north[1]=base[2];
  75                 north[2]=0;
  76
  77                 west[0]=base[4];
  78                 west[1]=base[5];
  79                 west[2]=base[6];
  80                 west[3]=0;
  81
  82                 if ((base[0]!='N' && base[0]!='S') || (base[3]!='W' && base[3]!='E'))
  83                 {
  84                         fprintf(stderr, "*** Error: \"%s\" doesn't look like a valid .hgt SRTM filename.\n", filename);
  85                         return -1;
  86                 }
  87
  88                 max_west=atoi(west);
  89
  90                 if (base[3]=='E')
  91                         max_west=360-max_west;
  92
  93                 min_west=max_west-1;
  94
  95                 if (max_west==360)
  96                         max_west=0;
  97
  98                 if (base[0]=='N')
  99                         min_north=atoi(north);
 100                 else
 101                         min_north=-atoi(north);
 102
 103                 max_north=min_north+1;
 104         }
 105
 106         if (bil)
 107         {
 108                 /* We obtain .BIL file coordinates
 109                    from the corresponding .BLW file */
 110
 111                 strncpy(blw_filename,filename,250);
 112                 x=strlen(filename);
 113
 114                 if (x>3)
 115                 {
 116                         blw_filename[x-2]='l';
 117                         blw_filename[x-1]='w';
 118                         blw_filename[x]=0;
 119
 120                         fd=fopen(blw_filename,"rb");
 121
 122                         if (fd!=NULL)
 123                         {
 124                                 n=fscanf(fd,"%lf",&cell_size);
 125
 126                                 if ((cell_size<0.0008) || (cell_size>0.0009))
 127                                 {
 128                                         printf("\n*** .BIL file's cell size is incompatible with SPLAT!!\n");
 129                                         exit(1);
 130                                 }
 131
 132                                 n=fscanf(fd,"%lf",&deg_west);
 133                                 n=fscanf(fd,"%lf",&deg_west);
 134                                 n=fscanf(fd,"%lf",&deg_west);
 135
 136                                 n=fscanf(fd,"%lf",&deg_west);
 137
 138                                 n=fscanf(fd,"%lf",&deg_north);
 139
 140                                 fclose(fd);
 141                         }
 142
 143                         min_north=(int)(deg_north);
 144                         max_north=min_north+1;
 145
 146                         if (deg_west<0.0)
 147                                 deg_west=-deg_west;
 148                         else
 149                                 deg_west=360.0-deg_west;
 150
 151                         min_west=(int)(deg_west);
 152
 153                         if (min_west==360)
 154                                 min_west=0;
 155
 156                         max_west=min_west+1;
 157                 }
 158         }
 159
 160         infile=open(filename, O_RDONLY);
 161
 162         if (infile==0)
 163         {
 164                 fprintf(stderr, "*** Error: Cannot open \"%s\"\n", filename);
 165                 return -1;
 166         }
 167
 168         n=read(infile,&buffer,2);
 169
 170         if ((buffer[0]=='P') && (buffer[1]=='K'))
 171         {
 172                 fprintf(stderr, "*** Error: \"%s\" still appears to be compressed!\n",filename);
 173                 close(infile);
 174                 return -1;
 175         }
 176
 177         lseek(infile,0L,SEEK_SET);
 178
 179         if (ippd==3600)
 180                 sprintf(sdf_filename, "%d:%d:%d:%d-hd.sdf", min_north, max_north, min_west, max_west);
 181         else
 182                 sprintf(sdf_filename, "%d:%d:%d:%d.sdf", min_north, max_north, min_west, max_west);
 183
 184         error=0;
 185         replacement_flag=0;
 186
 187         printf("Reading %s... ", filename);
 188         fflush(stdout);
 189
 190         for (x=0; (x<=ippd && error==0); x++)
 191                 for (y=0; (y<=ippd && error==0); y++)
 192                 {
 193                         bytes_read=read(infile,&buffer,2);
 194
 195                         if (bytes_read==2)
 196                         {
 197                                 if (bil)
 198                                 {
 199                                         /* "little-endian" structure */
 200
 201                                         byte=buffer[0]+(buffer[1]<<8);
 202
 203                                         if (buffer[1]&128)
 204                                                 byte-=0x10000;
 205                                 }
 206
 207                                 if (hgt)
 208                                 {
 209                                         /* "big-endian" structure */
 210
 211                                         byte=buffer[1]+(buffer[0]<<8);
 212
 213                                         if (buffer[0]&128)
 214                                                 byte-=0x10000;
 215                                 }
 216
 217                                 /* Flag problem elevations here */
 218
 219                                 if (byte<-32768)
 220                                         byte=-32768;
 221
 222                                 if (byte>32767)
 223                                         byte=32767;
 224
 225                                 if (byte<=min_elevation)
 226                                         replacement_flag=1;
 227
 228                                 srtm[x][y]=byte;
 229                         }
 230
 231                         else
 232                                 error=1;
 233                 }
 234
 235         if (error)
 236         {
 237                 fprintf(stderr,"\n*** Error: Premature EOF detected while reading \"%s\"!  :-(\n",filename);
 238         }
 239
 240         close(infile);
 241
 242         return 0;
 243 }
 244
 245 int LoadSDF_SDF(char *name)
 246 {
 247         /* This function reads uncompressed
 248            SPLAT Data Files (.sdf) into memory. */
 249
 250         int x, y, n, dummy;
 251         char sdf_file[255], path_plus_name[255];
 252         FILE *infile;
 253
 254         for (x=0; name[x]!='.' && name[x]!=0 && x<250; x++)
 255                 sdf_file[x]=name[x];
 256
 257         sdf_file[x]='.';
 258         sdf_file[x+1]='s';
 259         sdf_file[x+2]='d';
 260         sdf_file[x+3]='f';
 261         sdf_file[x+4]=0;
 262
 263         strncpy(path_plus_name,sdf_path,255);
 264         strncat(path_plus_name,sdf_file,255);
 265
 266         infile=fopen(path_plus_name,"rb");
 267
 268         if (infile==NULL)
 269                 return 0;
 270
 271         n=fscanf(infile,"%d", &dummy);
 272         n=fscanf(infile,"%d", &dummy);
 273         n=fscanf(infile,"%d", &dummy);
 274         n=fscanf(infile,"%d", &dummy);
 275
 276         printf("\nReading %s... ",path_plus_name);
 277         fflush(stdout);
 278
 279         for (x=0; x<1200; x++)
 280                 for (y=0; y<1200; y++)
 281                         n=fscanf(infile,"%d",&usgs[x][y]);
 282
 283         fclose(infile);
 284
 285         return 1;
 286 }
 287
 288 char *BZfgets(BZFILE *bzfd, unsigned length)
 289 {
 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. */
 297
 298         static int x, y, nBuf;
 299         static char buffer[BZBUFFER+1], output[BZBUFFER+1];
 300         char done=0;
 301
 302         if (opened!=1 && bzerror==BZ_OK)
 303         {
 304                 /* First time through.  Initialize everything! */
 305
 306                 x=0;
 307                 y=0;
 308                 nBuf=0;
 309                 opened=1;
 310                 output[0]=0;
 311         }
 312
 313         do
 314         {
 315                 if (x==nBuf && bzerror!=BZ_STREAM_END && bzerror==BZ_OK && opened)
 316                 {
 317                         /* Uncompress data into a static buffer */
 318
 319                         nBuf=BZ2_bzRead(&bzerror, bzfd, buffer, BZBUFFER);
 320                         buffer[nBuf]=0;
 321                         x=0;
 322                 }
 323
 324                 /* Build a string from buffer contents */
 325
 326                 output[y]=buffer[x];
 327
 328                 if (output[y]=='\n' || output[y]==0 || y==(int)length-1)
 329                 {
 330                         output[y+1]=0;
 331                         done=1;
 332                         y=0;
 333                 }
 334
 335                 else
 336                         y++;
 337                 x++;
 338
 339         } while (done==0);
 340
 341         if (output[0]==0)
 342                 opened=0;
 343
 344         return (output);
 345 }
 346
 347 int LoadSDF_BZ(char *name)
 348 {
 349         /* This function reads .bz2 compressed
 350            SPLAT Data Files into memory. */
 351
 352         int x, y, dummy;
 353         char sdf_file[255], path_plus_name[255];
 354         FILE *fd;
 355         BZFILE *bzfd;
 356
 357         for (x=0; name[x]!='.' && name[x]!=0 && x<247; x++)
 358                 sdf_file[x]=name[x];
 359
 360         sdf_file[x]='.';
 361         sdf_file[x+1]='s';
 362         sdf_file[x+2]='d';
 363         sdf_file[x+3]='f';
 364         sdf_file[x+4]='.';
 365         sdf_file[x+5]='b';
 366         sdf_file[x+6]='z';
 367         sdf_file[x+7]='2';
 368         sdf_file[x+8]=0;
 369
 370         strncpy(path_plus_name,sdf_path,255);
 371         strncat(path_plus_name,sdf_file,255);
 372
 373         fd=fopen(path_plus_name,"rb");
 374         bzfd=BZ2_bzReadOpen(&bzerror,fd,0,0,NULL,0);
 375
 376         if (fd!=NULL && bzerror==BZ_OK)
 377         {
 378                 printf("\nReading %s... ",path_plus_name);
 379                 fflush(stdout);
 380
 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);
 385
 386                 for (x=0; x<1200; x++)
 387                         for (y=0; y<1200; y++)
 388                                 sscanf(BZfgets(bzfd,20),"%d",&usgs[x][y]);
 389
 390                 fclose(fd);
 391
 392                 BZ2_bzReadClose(&bzerror,bzfd);
 393
 394                 return 1;
 395         }
 396
 397         else
 398                 return 0;
 399 }
 400
 401 char LoadSDF(char *name)
 402 {
 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. */
 408
 409         int return_value=-1;
 410
 411         /* Try to load an uncompressed SDF first. */
 412
 413         return_value=LoadSDF_SDF(name);
 414
 415         /* If that fails, try loading a compressed SDF. */
 416
 417         if (return_value==0 || return_value==-1)
 418                 return_value=LoadSDF_BZ(name);
 419
 420         return return_value;
 421 }
 422
 423 int ReadUSGS()
 424 {
 425         char usgs_filename[15];
 426
 427         /* usgs_filename is a minimal filename ("40:41:74:75").
 428            Full path and extentions are added later though
 429            subsequent function calls. */
 430
 431         sprintf(usgs_filename, "%d:%d:%d:%d", min_north, max_north, min_west, max_west);
 432
 433         return (LoadSDF(usgs_filename));
 434 }
 435
 436 void average_terrain(y,x,z)
 437 int x, y, z;
 438 {
 439         long accum;
 440         int temp=0, count, bad_value;
 441         double average;
 442
 443         bad_value=srtm[y][x];
 444
 445         accum=0L;
 446         count=0;
 447
 448         if (y>=2)
 449         {
 450                 temp=srtm[y-1][x];
 451
 452                 if (temp>bad_value)
 453                 {
 454                         accum+=temp;
 455                         count++;
 456                 }
 457         }
 458
 459         if (y<=mpi)
 460         {
 461                 temp=srtm[y+1][x];
 462
 463                 if (temp>bad_value)
 464                 {
 465                         accum+=temp;
 466                         count++;
 467                 }
 468         }
 469
 470         if ((y>=2) && (x<=(mpi-1)))
 471         {
 472                 temp=srtm[y-1][x+1];
 473
 474                 if (temp>bad_value)
 475                 {
 476                         accum+=temp;
 477                         count++;
 478                 }
 479         }
 480
 481         if (x<=(mpi-1))
 482         {
 483                 temp=srtm[y][x+1];
 484
 485                 if (temp>bad_value)
 486                 {
 487                         accum+=temp;
 488                         count++;
 489                 }
 490         }
 491
 492         if ((x<=(mpi-1)) && (y<=mpi))
 493         {
 494                 temp=srtm[y+1][x+1];
 495
 496                 if (temp>bad_value)
 497                 {
 498                         accum+=temp;
 499                         count++;
 500                 }
 501         }
 502
 503         if ((x>=1) && (y>=2))
 504         {
 505                 temp=srtm[y-1][x-1];
 506
 507                 if (temp>bad_value)
 508                 {
 509                         accum+=temp;
 510                         count++;
 511                 }
 512         }
 513
 514         if (x>=1)
 515         {
 516                 temp=srtm[y][x-1];
 517
 518                 if (temp>bad_value)
 519                 {
 520                         accum+=temp;
 521                         count++;
 522                 }
 523         }
 524
 525         if ((y<=mpi) && (x>=1))
 526         {
 527                 temp=srtm[y+1][x-1];
 528
 529                 if (temp>bad_value)
 530                 {
 531                         accum+=temp;
 532                         count++;
 533                 }
 534         }
 535
 536         if (count!=0)
 537         {
 538                 average=(((double)accum)/((double)count));
 539                 temp=(int)(average+0.5);
 540         }
 541
 542         if (temp>min_elevation)
 543                 srtm[y][x]=temp;
 544         else
 545                 srtm[y][x]=min_elevation;
 546 }
 547
 548 void WriteSDF(char *filename)
 549 {
 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.
 559          */
 560
 561         int x, y, byte, last_good_byte=0;
 562         FILE *outfile;
 563
 564         printf("\nWriting %s... ", filename);
 565         fflush(stdout);
 566
 567         outfile=fopen(filename,"wb");
 568
 569         fprintf(outfile, "%d\n%d\n%d\n%d\n", max_west, min_north, min_west, max_north);
 570
 571         for (y=ippd; y>=1; y--)         /* Omit the northern most edge */
 572                 for (x=mpi; x>=0; x--) /* Omit the eastern most edge  */
 573                 {
 574                         byte=srtm[y][x];
 575
 576                         if (byte>min_elevation)
 577                                 last_good_byte=byte;
 578
 579                         if (byte<min_elevation)
 580                         {
 581                                 if (merge)
 582                                 {
 583                                         if (ippd==3600)
 584                                                 fprintf(outfile,"%d\n",usgs[1200-(y/3)][1199-(x/3)]);
 585                                         else
 586                                                 fprintf(outfile,"%d\n",usgs[1200-y][1199-x]);
 587                                 }
 588
 589                                 else
 590                                 {
 591                                         average_terrain(y,x,last_good_byte);
 592                                         fprintf(outfile,"%d\n",srtm[y][x]);
 593                                 }
 594                         }
 595                         else
 596                                 fprintf(outfile,"%d\n",byte);
 597                 }
 598
 599         printf("Done!\n");
 600
 601         fclose(outfile);
 602 }
 603
 604 int main(int argc, char **argv)
 605 {
 606         int x, y, z=0;
 607         char *env=NULL, string[255], *s=NULL;
 608         FILE *fd;
 609
 610         if (strstr(argv[0], "srtm2sdf-hd")!=NULL)
 611         {
 612                 ippd=3600;      /* High Definition (1 arc-sec) Mode */
 613                 strncpy(string,"srtm2sdf-hd\0",12);
 614         }
 615
 616         else
 617         {
 618                 ippd=1200;      /* Standard Definition (3 arc-sec) Mode */
 619                 strncpy(string,"srtm2sdf\0",9);
 620         }
 621
 622         mpi=ippd-1;             /* Maximum pixel index per degree */
 623
 624         if (argc==1 || (argc==2 && strncmp(argv[1],"-h",2)==0))
 625         {
 626                 if (ippd==1200)
 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");
 628
 629                 if (ippd==3600)
 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");
 631
 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);
 639
 640                 return 1;
 641         }
 642
 643         y=argc-1;
 644
 645         min_elevation=0;
 646
 647         for (x=1, z=0; x<=y; x++)
 648         {
 649                 if (strcmp(argv[x],"-d")==0)
 650                 {
 651                         z=x+1;
 652
 653                         if (z<=y && argv[z][0] && argv[z][0]!='-')
 654                                 strncpy(sdf_path,argv[z],253);
 655                 }
 656
 657                 if (strcmp(argv[x],"-n")==0)
 658                 {
 659                         z=x+1;
 660
 661                         if (z<=y && argv[z][0])
 662                         {
 663                                 sscanf(argv[z],"%d",&min_elevation);
 664
 665                                 if (min_elevation<-32767)
 666                                         min_elevation=0;
 667                         }
 668                 }
 669         }
 670
 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. */
 675
 676         if (sdf_path[0]==0)
 677         {
 678                 env=getenv("HOME");
 679
 680                 sprintf(string,"%s/.splat_path",env);
 681
 682                 fd=fopen(string,"r");
 683
 684                 if (fd!=NULL)
 685                 {
 686                         s=fgets(string,253,fd);
 687
 688                         /* Remove <CR> and/or <LF> from string */
 689
 690                         for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0 && x<253; x++);
 691                         string[x]=0;
 692
 693                         strncpy(sdf_path,string,253);
 694
 695                         fclose(fd);
 696                 }
 697         }
 698
 699         /* Ensure a trailing '/' is present in sdf_path */
 700
 701         if (sdf_path[0])
 702         {
 703                 x=strlen(sdf_path);
 704
 705                 if (sdf_path[x-1]!='/' && x!=0)
 706                 {
 707                         sdf_path[x]='/';
 708                         sdf_path[x+1]=0;
 709                 }
 710         }
 711
 712         if (ReadSRTM(argv[z+1])==0)
 713         {
 714                 if (replacement_flag && sdf_path[0])
 715                         merge=ReadUSGS();
 716
 717                 WriteSDF(sdf_filename);
 718         }
 719
 720         return 0;
 721 }
 722