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