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