DEM/dem.cxx

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