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