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