DEM/dem.c

   1 // -*- Mode: C++ -*-
   2 //
   3 // dem.c -- DEM management class
   4 //
   5 // Written by Curtis Olson, started March 1998.
   6 //
   7 // Copyright (C) 1998  Curtis L. Olson  - curt@me.umn.edu
   8 //
   9 // This program is free software; you can redistribute it and/or
  10 // modify it under the terms of the GNU General Public License as
  11 // published by the Free Software Foundation; either version 2 of the
  12 // License, or (at your option) any later version.
  13 //
  14 // This program is distributed in the hope that it will be useful, but
  15 // WITHOUT ANY WARRANTY; without even the implied warranty of
  16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17 // General Public License for more details.
  18 //
  19 // You should have received a copy of the GNU General Public License
  20 // along with this program; if not, write to the Free Software
  21 // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  22 //
  23 // $Id$
  24 // (Log is kept at end of this file)
  25
  26
  27 #include <ctype.h>    // isspace()
  28 #include <math.h>     // rint()
  29 #include <stdio.h>
  30 #include <stdlib.h>   // atoi()
  31 #include <sys/stat.h> // stat()
  32 #include <unistd.h>   // stat()
  33
  34 #include "dem.h"
  35 #include "leastsqs.h"
  36
  37 #include <Include/fg_constants.h>
  38
  39
  40 fgDEM::fgDEM( void ) {
  41     // printf("class fgDEM CONstructor called.\n");
  42 }
  43
  44
  45 // open a DEM file
  46 int fgDEM::open ( char *file ) {
  47     // open input file (or read from stdin)
  48     if ( strcmp(file, "-") == 0 ) {
  49         printf("Loading DEM data file: stdin\n");
  50         fd = stdin;
  51     } else {
  52         if ( (fd = fopen(file, "r")) == NULL ) {
  53             printf("Cannot open %s\n", file);
  54             return(0);
  55         }
  56         printf("Loading DEM data file: %s\n", file);
  57     }
  58
  59     return(1);
  60 }
  61
  62
  63 // close a DEM file
  64 int fgDEM::close ( void ) {
  65     fclose(fd);
  66
  67     return(1);
  68 }
  69
  70
  71 // return next token from input stream
  72 static void next_token(FILE *fd, char *token) {
  73     int result;
  74
  75     result = fscanf(fd, "%s", token);
  76
  77     if ( result == EOF ) {
  78         strcpy(token, "__END_OF_FILE__");
  79         printf("    Warning:  Reached end of file!\n");
  80     }
  81
  82     // printf("    returning %s\n", token);
  83 }
  84
  85
  86 // return next integer from input stream
  87 static int next_int(FILE *fd) {
  88     char token[80];
  89
  90     next_token(fd, token);
  91     return ( atoi(token) );
  92 }
  93
  94
  95 // return next double from input stream
  96 double next_double(FILE *fd) {
  97     char token[80];
  98
  99     next_token(fd, token);
 100     return ( atof(token) );
 101 }
 102
 103
 104 // return next exponential num from input stream
 105 int next_exp(FILE *fd) {
 106     double mantissa;
 107     int exp, acc;
 108     int i;
 109
 110     fscanf(fd, "%lfD%d", &mantissa, &exp);
 111
 112     // printf("    Mantissa = %.4f  Exp = %d\n", mantissa, exp);
 113
 114     acc = 1;
 115     if ( exp > 0 ) {
 116         for ( i = 1; i <= exp; i++ ) {
 117             acc *= 10;
 118         }
 119     } else if ( exp < 0 ) {
 120         for ( i = -1; i >= exp; i-- ) {
 121             acc /= 10;
 122         }
 123     }
 124
 125     return( (int)rint(mantissa * (double)acc) );
 126 }
 127
 128
 129 // read and parse DEM "A" record
 130 void fgDEM::read_a_record( void ) {
 131     int i, inum;
 132     double dnum;
 133     char name[144];
 134     char token[80];
 135     char *ptr;
 136
 137     // get the name field (144 characters)
 138     for ( i = 0; i < 144; i++ ) {
 139         name[i] = fgetc(fd);
 140     }
 141     name[i+1] = '\0';
 142
 143     // clean off the whitespace at the end
 144     for ( i = strlen(name)-2; i > 0; i-- ) {
 145         if ( !isspace(name[i]) ) {
 146             i=0;
 147         } else {
 148             name[i] = '\0';
 149         }
 150     }
 151     printf("    Quad name field: %s\n", name);
 152
 153     // DEM level code, 3 reflects processing by DMA
 154     inum = next_int(fd);
 155     printf("    DEM level code = %d\n", inum);
 156
 157     // Pattern code, 1 indicates a regular elevation pattern
 158     inum = next_int(fd);
 159     printf("    Pattern code = %d\n", inum);
 160
 161     // Planimetric reference system code, 0 indicates geographic
 162     // coordinate system.
 163     inum = next_int(fd);
 164     printf("    Planimetric reference code = %d\n", inum);
 165
 166     // Zone code
 167     inum = next_int(fd);
 168     printf("    Zone code = %d\n", inum);
 169
 170     // Map projection parameters (ignored)
 171     for ( i = 0; i < 15; i++ ) {
 172         dnum = next_double(fd);
 173         // printf("%d: %f\n",i,dnum);
 174     }
 175
 176     // Units code, 3 represents arc-seconds as the unit of measure for
 177     // ground planimetric coordinates throughout the file.
 178     inum = next_int(fd);
 179     if ( inum != 3 ) {
 180         printf("    Unknown (X,Y) units code = %d!\n", inum);
 181         exit(-1);
 182     }
 183
 184     // Units code; 2 represents meters as the unit of measure for
 185     // elevation coordinates throughout the file.
 186     inum = next_int(fd);
 187     if ( inum != 2 ) {
 188         printf("    Unknown (Z) units code = %d!\n", inum);
 189         exit(-1);
 190     }
 191
 192     // Number (n) of sides in the polygon which defines the coverage of
 193     // the DEM file (usually equal to 4).
 194     inum = next_int(fd);
 195     if ( inum != 4 ) {
 196         printf("    Unknown polygon dimension = %d!\n", inum);
 197         exit(-1);
 198     }
 199
 200     // Ground coordinates of bounding box in arc-seconds
 201     dem_x1 = originx = next_exp(fd);
 202     dem_y1 = originy = next_exp(fd);
 203     printf("    Origin = (%.2f,%.2f)\n", originx, originy);
 204
 205     dem_x2 = next_exp(fd);
 206     dem_y2 = next_exp(fd);
 207
 208     dem_x3 = next_exp(fd);
 209     dem_y3 = next_exp(fd);
 210
 211     dem_x4 = next_exp(fd);
 212     dem_y4 = next_exp(fd);
 213
 214     // Minimum/maximum elevations in meters
 215     dem_z1 = next_exp(fd);
 216     dem_z2 = next_exp(fd);
 217     printf("    Elevation range %.4f %.4f\n", dem_z1, dem_z2);
 218
 219     // Counterclockwise angle from the primary axis of ground
 220     // planimetric referenced to the primary axis of the DEM local
 221     // reference system.
 222     next_token(fd, token);
 223
 224     // Accuracy code; 0 indicates that a record of accuracy does not
 225     // exist and that no record type C will follow.
 226
 227     // DEM spacial resolution.  Usually (3,3,1) (3,6,1) or (3,9,1)
 228     // depending on latitude
 229
 230     // I will eventually have to do something with this for data at
 231     // higher latitudes */
 232     next_token(fd, token);
 233     printf("    accuracy & spacial resolution string = %s\n", token);
 234     i = strlen(token);
 235     printf("    length = %d\n", i);
 236
 237     ptr = token + i - 12;
 238     printf("    last field = %s = %.2f\n", ptr, atof(ptr));
 239     ptr[0] = '\0';
 240
 241     ptr = ptr - 12;
 242     col_step = atof(ptr);
 243     printf("    last field = %s = %.2f\n", ptr, col_step);
 244     ptr[0] = '\0';
 245
 246     ptr = ptr - 12;
 247     row_step = atof(ptr);
 248     printf("    last field = %s = %.2f\n", ptr, row_step);
 249     ptr[0] = '\0';
 250
 251     // accuracy code = atod(token)
 252     inum = atoi(token);
 253     printf("    Accuracy code = %d\n", inum);
 254
 255     printf("    column step = %.2f  row step = %.2f\n",
 256            col_step, row_step);
 257     // dimension of arrays to follow (1)
 258     next_token(fd, token);
 259
 260     // number of profiles
 261     dem_num_profiles = rows = next_int(fd);
 262     printf("    Expecting %d profiles\n", dem_num_profiles);
 263 }
 264
 265
 266 // read and parse DEM "B" record
 267 void fgDEM::read_b_record(float dem_data[DEM_SIZE_1][DEM_SIZE_1])
 268 {
 269     char token[80];
 270     int i;
 271
 272     // row / column id of this profile
 273     prof_col = next_int(fd);
 274     prof_row = next_int(fd);
 275     // printf("col id = %d  row id = %d\n", prof_col, prof_row);
 276
 277     // Number of columns and rows (elevations) in this profile
 278     prof_num_cols = cols = next_int(fd);
 279     prof_num_rows = next_int(fd);
 280     // printf("    profile num rows = %d\n", prof_num_cols);
 281
 282     // Ground planimetric coordinates (arc-seconds) of the first
 283     // elevation in the profile
 284     prof_x1 = next_exp(fd);
 285     prof_y1 = next_exp(fd);
 286     // printf("    Starting at %.2f %.2f\n", prof_x1, prof_y1);
 287
 288     // Elevation of local datum for the profile.  Always zero for
 289     // 1-degree DEM, the reference is mean sea level.
 290     next_token(fd, token);
 291
 292     // Minimum and maximum elevations for the profile.
 293     next_token(fd, token);
 294     next_token(fd, token);
 295
 296     // One (usually) dimensional array (prof_num_cols,1) of elevations
 297     for ( i = 0; i < prof_num_cols; i++ ) {
 298         prof_data = next_int(fd);
 299         dem_data[i][cur_row] = (float)prof_data;
 300     }
 301 }
 302
 303
 304 // parse dem file
 305 int fgDEM::parse( float dem_data[DEM_SIZE_1][DEM_SIZE_1] ) {
 306     int i;
 307
 308     cur_row = 0;
 309
 310     read_a_record();
 311
 312     for ( i = 0; i < dem_num_profiles; i++ ) {
 313         read_b_record( dem_data );
 314         cur_row++;
 315
 316         if ( cur_row % 100 == 0 ) {
 317             printf("    loaded %d profiles of data\n", cur_row);
 318         }
 319     }
 320
 321     printf("    Done parsing\n");
 322
 323     return(0);
 324 }
 325
 326
 327 // return the current altitude based on mesh data.  We should rewrite
 328 // this to interpolate exact values, but for now this is good enough
 329 double fgDEM::interpolate_altitude( float dem_data[DEM_SIZE_1][DEM_SIZE_1],
 330                                     double lon, double lat)
 331 {
 332     // we expect incoming (lon,lat) to be in arcsec for now
 333
 334     double xlocal, ylocal, dx, dy, zA, zB, elev;
 335     int x1, x2, x3, y1, y2, y3;
 336     float z1, z2, z3;
 337     int xindex, yindex;
 338
 339     /* determine if we are in the lower triangle or the upper triangle
 340        ______
 341        |   /|
 342        |  / |
 343        | /  |
 344        |/   |
 345        ------
 346
 347        then calculate our end points
 348      */
 349
 350     xlocal = (lon - originx) / col_step;
 351     ylocal = (lat - originy) / row_step;
 352
 353     xindex = (int)(xlocal);
 354     yindex = (int)(ylocal);
 355
 356     // printf("xindex = %d  yindex = %d\n", xindex, yindex);
 357
 358     if ( xindex + 1 == cols ) {
 359         xindex--;
 360     }
 361
 362     if ( yindex + 1 == rows ) {
 363         yindex--;
 364     }
 365
 366     if ( (xindex < 0) || (xindex + 1 >= cols) ||
 367          (yindex < 0) || (yindex + 1 >= rows) ) {
 368         return(-9999);
 369     }
 370
 371     dx = xlocal - xindex;
 372     dy = ylocal - yindex;
 373
 374     if ( dx > dy ) {
 375         // lower triangle
 376         // printf("  Lower triangle\n");
 377
 378         x1 = xindex;
 379         y1 = yindex;
 380         z1 = dem_data[x1][y1];
 381
 382         x2 = xindex + 1;
 383         y2 = yindex;
 384         z2 = dem_data[x2][y2];
 385
 386         x3 = xindex + 1;
 387         y3 = yindex + 1;
 388         z3 = dem_data[x3][y3];
 389
 390         // printf("  dx = %.2f  dy = %.2f\n", dx, dy);
 391         // printf("  (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1);
 392         // printf("  (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2);
 393         // printf("  (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3);
 394
 395         zA = dx * (z2 - z1) + z1;
 396         zB = dx * (z3 - z1) + z1;
 397
 398         // printf("  zA = %.2f  zB = %.2f\n", zA, zB);
 399
 400         if ( dx > FG_EPSILON ) {
 401             elev = dy * (zB - zA) / dx + zA;
 402         } else {
 403             elev = zA;
 404         }
 405     } else {
 406         // upper triangle
 407         // printf("  Upper triangle\n");
 408
 409         x1 = xindex;
 410         y1 = yindex;
 411         z1 = dem_data[x1][y1];
 412
 413         x2 = xindex;
 414         y2 = yindex + 1;
 415         z2 = dem_data[x2][y2];
 416
 417         x3 = xindex + 1;
 418         y3 = yindex + 1;
 419         z3 = dem_data[x3][y3];
 420
 421         // printf("  dx = %.2f  dy = %.2f\n", dx, dy);
 422         // printf("  (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1);
 423         // printf("  (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2);
 424         // printf("  (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3);
 425
 426         zA = dy * (z2 - z1) + z1;
 427         zB = dy * (z3 - z1) + z1;
 428
 429         // printf("  zA = %.2f  zB = %.2f\n", zA, zB );
 430         // printf("  xB - xA = %.2f\n", col_step * dy / row_step);
 431
 432         if ( dy > FG_EPSILON ) {
 433             elev = dx * (zB - zA) / dy    + zA;
 434         } else {
 435             elev = zA;
 436         }
 437     }
 438
 439     return(elev);
 440 }
 441
 442
 443 // Use least squares to fit a simpler data set to dem data
 444 void fgDEM::fit( float dem_data[DEM_SIZE_1][DEM_SIZE_1],
 445                  float output_data[DEM_SIZE_1][DEM_SIZE_1],
 446                  char *fg_root, double error, struct fgBUCKET *p )
 447 {
 448     double x[DEM_SIZE_1], y[DEM_SIZE_1];
 449     double m, b, ave_error, max_error;
 450     double cury, lasty;
 451     int n, row, start, end, good_fit;
 452     int colmin, colmax, rowmin, rowmax;
 453     // FILE *dem, *fit, *fit1;
 454
 455     printf("Initializing output mesh structure\n");
 456     outputmesh_init( output_data );
 457
 458     // determine dimensions
 459     colmin = p->x * ( (cols - 1) / 8);
 460     colmax = colmin + ( (cols - 1) / 8);
 461     rowmin = p->y * ( (rows - 1) / 8);
 462     rowmax = rowmin + ( (rows - 1) / 8);
 463     printf("Fitting region = %d,%d to %d,%d\n", colmin, rowmin, colmax, rowmax);
 464
 465     // include the corners explicitly
 466     outputmesh_set_pt(output_data, colmin, rowmin, dem_data[colmin][rowmin]);
 467     outputmesh_set_pt(output_data, colmin, rowmax, dem_data[colmin][rowmax]);
 468     outputmesh_set_pt(output_data, colmax, rowmax, dem_data[colmax][rowmax]);
 469     outputmesh_set_pt(output_data, colmax, rowmin, dem_data[colmax][rowmin]);
 470
 471     printf("Beginning best fit procedure\n");
 472
 473     for ( row = rowmin; row <= rowmax; row++ ) {
 474         // fit  = fopen("fit.dat",  "w");
 475         // fit1 = fopen("fit1.dat", "w");
 476
 477         start = colmin;
 478
 479         // printf("    fitting row = %d\n", row);
 480
 481         while ( start < colmax ) {
 482             end = start + 1;
 483             good_fit = 1;
 484
 485             x[(end - start) - 1] = 0.0 + ( start * col_step );
 486             y[(end - start) - 1] = dem_data[start][row];
 487
 488             while ( (end <= colmax) && good_fit ) {
 489                 n = (end - start) + 1;
 490                 // printf("Least square of first %d points\n", n);
 491                 x[end - start] = 0.0 + ( end * col_step );
 492                 y[end - start] = dem_data[end][row];
 493                 least_squares(x, y, n, &m, &b);
 494                 ave_error = least_squares_error(x, y, n, m, b);
 495                 max_error = least_squares_max_error(x, y, n, m, b);
 496
 497                 /*
 498                 printf("%d - %d  ave error = %.2f  max error = %.2f  y = %.2f*x + %.2f\n",
 499                 start, end, ave_error, max_error, m, b);
 500
 501                 f = fopen("gnuplot.dat", "w");
 502                 for ( j = 0; j <= end; j++) {
 503                     fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ),
 504                             dem_data[row][j]);
 505                 }
 506                 for ( j = start; j <= end; j++) {
 507                     fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ),
 508                             dem_data[row][j]);
 509                 }
 510                 fclose(f);
 511
 512                 printf("Please hit return: "); gets(junk);
 513                 */
 514
 515                 if ( max_error > error ) {
 516                     good_fit = 0;
 517                 }
 518
 519                 end++;
 520             }
 521
 522             if ( !good_fit ) {
 523                 // error exceeded the threshold, back up
 524                 end -= 2;  // back "end" up to the last good enough fit
 525                 n--;       // back "n" up appropriately too
 526             } else {
 527                 // we popped out of the above loop while still within
 528                 // the error threshold, so we must be at the end of
 529                 // the data set
 530                 end--;
 531             }
 532
 533             least_squares(x, y, n, &m, &b);
 534             ave_error = least_squares_error(x, y, n, m, b);
 535             max_error = least_squares_max_error(x, y, n, m, b);
 536
 537             /*
 538             printf("\n");
 539             printf("%d - %d  ave error = %.2f  max error = %.2f  y = %.2f*x + %.2f\n",
 540                    start, end, ave_error, max_error, m, b);
 541             printf("\n");
 542
 543             fprintf(fit1, "%.2f %.2f\n", x[0], m * x[0] + b);
 544             fprintf(fit1, "%.2f %.2f\n", x[end-start], m * x[end-start] + b);
 545             */
 546
 547             if ( start > colmin ) {
 548                 // skip this for the first line segment
 549                 cury = m * x[0] + b;
 550                 outputmesh_set_pt(output_data, start, row, (lasty + cury) / 2);
 551                 // fprintf(fit, "%.2f %.2f\n", x[0], (lasty + cury) / 2);
 552             }
 553
 554             lasty = m * x[end-start] + b;
 555             start = end;
 556         }
 557
 558         /*
 559         fclose(fit);
 560         fclose(fit1);
 561
 562         dem = fopen("gnuplot.dat", "w");
 563         for ( j = 0; j < DEM_SIZE_1; j++) {
 564             fprintf(dem, "%.2f %.2f\n", 0.0 + ( j * col_step ),
 565                     dem_data[j][row]);
 566         }
 567         fclose(dem);
 568         */
 569
 570         // NOTICE, this is for testing only.  This instance of
 571         // output_nodes should be removed.  It should be called only
 572         // once at the end once all the nodes have been generated.
 573         // newmesh_output_nodes(&nm, "mesh.node");
 574         // printf("Please hit return: "); gets(junk);
 575     }
 576
 577     outputmesh_output_nodes(output_data, fg_root, p);
 578 }
 579
 580
 581 // Initialize output mesh structure
 582 void fgDEM::outputmesh_init( float output_data[DEM_SIZE_1][DEM_SIZE_1] ) {
 583     int i, j;
 584
 585     for ( i = 0; i < DEM_SIZE_1; i++ ) {
 586         for ( j = 0; j < DEM_SIZE_1; j++ ) {
 587             output_data[i][j] = -9999.0;
 588         }
 589     }
 590 }
 591
 592
 593 // Get the value of a mesh node
 594 double fgDEM::outputmesh_get_pt( float output_data[DEM_SIZE_1][DEM_SIZE_1],
 595                                  int i, int j )
 596 {
 597     return ( output_data[i][j] );
 598 }
 599
 600
 601 // Set the value of a mesh node
 602 void fgDEM::outputmesh_set_pt( float output_data[DEM_SIZE_1][DEM_SIZE_1],
 603                                int i, int j, double value )
 604 {
 605     // printf("Setting data[%d][%d] = %.2f\n", i, j, value);
 606    output_data[i][j] = value;
 607 }
 608
 609
 610 // Write out a node file that can be used by the "triangle" program
 611 void fgDEM::outputmesh_output_nodes( float output_data[DEM_SIZE_1][DEM_SIZE_1],
 612                                      char *fg_root, struct fgBUCKET *p )
 613 {
 614     struct stat stat_buf;
 615     char base_path[256], dir[256], file[256];
 616     char command[256];
 617     FILE *fd;
 618     long int index;
 619     int colmin, colmax, rowmin, rowmax;
 620     int i, j, count, result;
 621
 622     // determine dimensions
 623     colmin = p->x * ( (cols - 1) / 8);
 624     colmax = colmin + ( (cols - 1) / 8);
 625     rowmin = p->y * ( (rows - 1) / 8);
 626     rowmax = rowmin + ( (rows - 1) / 8);
 627     printf("  dumping region = %d,%d to %d,%d\n",
 628            colmin, rowmin, colmax, rowmax);
 629
 630     // generate the base directory
 631     fgBucketGenBasePath(p, base_path);
 632     printf("fg_root = %s  Base Path = %s\n", fg_root, base_path);
 633     sprintf(dir, "%s/Scenery/%s", fg_root, base_path);
 634     printf("Dir = %s\n", dir);
 635
 636     // stat() directory and create if needed
 637     result = stat(dir, &stat_buf);
 638     if ( result != 0 ) {
 639         printf("Stat error need to create directory\n");
 640         sprintf(command, "mkdir -p %s\n", dir);
 641         system(command);
 642     } else {
 643         // assume directory exists
 644     }
 645
 646     // get index and generate output file name
 647     index = fgBucketGenIndex(p);
 648     sprintf(file, "%s/%ld.node", dir, index);
 649
 650     printf("Creating node file:  %s\n", file);
 651     fd = fopen(file, "w");
 652
 653     // first count nodes to generate header
 654     count = 0;
 655     for ( j = rowmin; j <= rowmax; j++ ) {
 656         for ( i = colmin; i <= colmax; i++ ) {
 657             if ( output_data[i][j] > -9000.0 ) {
 658                 count++;
 659             }
 660         }
 661         // printf("    count = %d\n", count);
 662     }
 663     fprintf(fd, "%d 2 1 0\n", count);
 664
 665     // now write out actual node data
 666     count = 1;
 667     for ( j = rowmin; j <= rowmax; j++ ) {
 668         for ( i = colmin; i <= colmax; i++ ) {
 669             if ( output_data[i][j] > -9000.0 ) {
 670                 fprintf(fd, "%d %.2f %.2f %.2f\n",
 671                         count++,
 672                         originx + (double)i * col_step,
 673                         originy + (double)j * row_step,
 674                         output_data[i][j]);
 675             }
 676         }
 677         // printf("    count = %d\n", count);
 678     }
 679
 680     fclose(fd);
 681 }
 682
 683
 684 fgDEM::~fgDEM( void ) {
 685     // printf("class fgDEM DEstructor called.\n");
 686 }
 687
 688
 689 // $Log$
 690 // Revision 1.2  1998/03/23 20:35:41  curt
 691 // Updated to use FG_EPSILON
 692 //
 693 // Revision 1.1  1998/03/19 02:54:47  curt
 694 // Reorganized into a class lib called fgDEM.
 695 //
 696 // Revision 1.1  1998/03/19 01:46:28  curt
 697 // Initial revision.
 698 //