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