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