DEM/dem.cxx

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