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