3 // dem.c -- DEM management class
5 // Written by Curtis Olson, started March 1998.
7 // Copyright (C) 1998 Curtis L. Olson - curt@me.umn.edu
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.
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.
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.
24 // (Log is kept at end of this file)
31 #include <ctype.h> // isspace()
32 #include <stdlib.h> // atoi()
33 #include <math.h> // rint()
36 #include <sys/stat.h> // stat()
39 # include <unistd.h> // stat()
43 // #include <zlib/zlib.h>
44 #include <Misc/fgstream.hxx>
45 #include <Misc/strutils.hxx>
46 #include <Include/compiler.h>
47 FG_USING_NAMESPACE(std);
50 #include "leastsqs.hxx"
52 #include <Include/fg_constants.h>
55 #define MAX_EX_NODES 10000
60 # define MKDIR(a) mkdir(a)
62 # define MKDIR(a) mkdir(a,S_IRWXU) // I am just guessing at this flag (NHV)
63 # endif // __BORLANDC__
67 fgDEM::fgDEM( void ) {
68 // printf("class fgDEM CONstructor called.\n");
69 dem_data = new float[DEM_SIZE_1][DEM_SIZE_1];
70 output_data = new float[DEM_SIZE_1][DEM_SIZE_1];
76 // return the file path name ( foo/bar/file.ext = foo/bar )
77 static void extract_path ( const char *in, char *base) {
84 while ( (i >= 0) && (in[i] != '/') ) {
92 // Make a subdirectory
93 static int my_mkdir (const char *dir) {
99 result = stat (dir, &stat_buf);
103 result = stat (dir, &stat_buf);
105 printf ("problem creating %s\n", dir);
107 printf ("%s created\n", dir);
110 printf ("%s already exists\n", dir);
120 int fgDEM::open ( const string& file ) {
121 // open input file (or read from stdin)
123 printf("Loading DEM data file: stdin\n");
125 // fd = gzdopen(STDIN_FILENO, "r");
126 printf("Not yet ported ...\n");
129 in = new fg_gzifstream( file );
131 cout << "Cannot open " << file << endl;
134 cout << "Loading DEM data file: " << file << endl;
142 int fgDEM::close () {
143 // the fg_gzifstream doesn't seem to have a close()
151 // return next token from input stream
152 string fgDEM::next_token() {
155 in->stream() >> token;
157 // cout << " returning " + token + "\n";
163 // return next integer from input stream
164 int fgDEM::next_int() {
167 in->stream() >> result;
173 // return next double from input stream
174 double fgDEM::next_double() {
177 in->stream() >> result;
183 // return next exponential num from input stream
184 double fgDEM::next_exp() {
190 token = next_token();
193 const char* p = token.c_str();
197 for ( ; *p != 0; ++p )
205 return ::atof( buf );
207 sscanf(token.c_str(), "%lfD%d", &mantissa, &exp);
209 // cout << " Mantissa = " << mantissa << " Exp = " << exp << "\n";
213 for ( i = 1; i <= exp; i++ ) {
216 } else if ( exp < 0 ) {
217 for ( i = -1; i >= exp; i-- ) {
222 return( (int)rint(mantissa * (double)acc) );
227 // read and parse DEM "A" record
228 int fgDEM::read_a_record() {
235 // get the name field (144 characters)
236 for ( i = 0; i < 144; i++ ) {
241 // clean off the trailing whitespace
243 cout << " Quad name field: " << name << endl;
245 // DEM level code, 3 reflects processing by DMA
247 cout << " DEM level code = " << inum << "\n";
253 // Pattern code, 1 indicates a regular elevation pattern
255 cout << " Pattern code = " << inum << "\n";
257 // Planimetric reference system code, 0 indicates geographic
258 // coordinate system.
260 cout << " Planimetric reference code = " << inum << "\n";
264 cout << " Zone code = " << inum << "\n";
266 // Map projection parameters (ignored)
267 for ( i = 0; i < 15; i++ ) {
269 // printf("%d: %f\n",i,dnum);
272 // Units code, 3 represents arc-seconds as the unit of measure for
273 // ground planimetric coordinates throughout the file.
276 cout << " Unknown (X,Y) units code = " << inum << "!\n";
280 // Units code; 2 represents meters as the unit of measure for
281 // elevation coordinates throughout the file.
284 cout << " Unknown (Z) units code = " << inum << "!\n";
288 // Number (n) of sides in the polygon which defines the coverage of
289 // the DEM file (usually equal to 4).
292 cout << " Unknown polygon dimension = " << inum << "!\n";
296 // Ground coordinates of bounding box in arc-seconds
297 dem_x1 = originx = next_exp();
298 dem_y1 = originy = next_exp();
299 cout << " Origin = (" << originx << "," << originy << ")\n";
310 // Minimum/maximum elevations in meters
313 cout << " Elevation range " << dem_z1 << " to " << dem_z2 << "\n";
315 // Counterclockwise angle from the primary axis of ground
316 // planimetric referenced to the primary axis of the DEM local
318 token = next_token();
320 // Accuracy code; 0 indicates that a record of accuracy does not
321 // exist and that no record type C will follow.
323 // DEM spacial resolution. Usually (3,3,1) (3,6,1) or (3,9,1)
324 // depending on latitude
326 // I will eventually have to do something with this for data at
327 // higher latitudes */
328 token = next_token();
329 cout << " accuracy & spacial resolution string = " << token << endl;
331 cout << " length = " << i << "\n";
334 inum = atoi( token.substr( 0, i - 36 ) );
335 row_step = atof( token.substr( i - 36, 12 ) );
336 col_step = atof( token.substr( i - 24, 12 ) );
337 //token.substr( 25, 12 )
339 ptr = token.c_str() + i - 12;
340 cout << " last field = " << ptr << " = " << atof(ptr) << "\n";
344 col_step = atof(ptr);
345 cout << " last field = " << ptr << " = " << col_step << "\n";
349 row_step = atof(ptr);
350 cout << " last field = " << ptr << " = " << row_step << "\n";
353 // accuracy code = atod(token)
357 cout << " Accuracy code = " << inum << "\n";
359 cout << " column step = " << col_step <<
360 " row step = " << row_step << "\n";
362 // dimension of arrays to follow (1)
363 token = next_token();
365 // number of profiles
366 dem_num_profiles = cols = next_int();
367 cout << " Expecting " << dem_num_profiles << " profiles\n";
373 // read and parse DEM "B" record
374 void fgDEM::read_b_record( ) {
378 // row / column id of this profile
379 prof_row = next_int();
380 prof_col = next_int();
381 // printf("col id = %d row id = %d\n", prof_col, prof_row);
383 // Number of columns and rows (elevations) in this profile
384 prof_num_rows = rows = next_int();
385 prof_num_cols = next_int();
386 // printf(" profile num rows = %d\n", prof_num_rows);
388 // Ground planimetric coordinates (arc-seconds) of the first
389 // elevation in the profile
390 prof_x1 = next_exp();
391 prof_y1 = next_exp();
392 // printf(" Starting at %.2f %.2f\n", prof_x1, prof_y1);
394 // Elevation of local datum for the profile. Always zero for
395 // 1-degree DEM, the reference is mean sea level.
396 token = next_token();
398 // Minimum and maximum elevations for the profile.
399 token = next_token();
400 token = next_token();
402 // One (usually) dimensional array (prof_num_cols,1) of elevations
403 for ( i = 0; i < prof_num_rows; i++ ) {
404 prof_data = next_int();
405 dem_data[cur_col][i] = (float)prof_data;
411 int fgDEM::parse( ) {
416 if ( !read_a_record() ) {
420 for ( i = 0; i < dem_num_profiles; i++ ) {
421 // printf("Ready to read next b record\n");
425 if ( cur_col % 100 == 0 ) {
426 cout << " loaded " << cur_col << " profiles of data\n";
430 cout << " Done parsing\n";
436 // return the current altitude based on mesh data. We should rewrite
437 // this to interpolate exact values, but for now this is good enough
438 double fgDEM::interpolate_altitude( double lon, double lat ) {
439 // we expect incoming (lon,lat) to be in arcsec for now
441 double xlocal, ylocal, dx, dy, zA, zB, elev;
442 int x1, x2, x3, y1, y2, y3;
446 /* determine if we are in the lower triangle or the upper triangle
454 then calculate our end points
457 xlocal = (lon - originx) / col_step;
458 ylocal = (lat - originy) / row_step;
460 xindex = (int)(xlocal);
461 yindex = (int)(ylocal);
463 // printf("xindex = %d yindex = %d\n", xindex, yindex);
465 if ( xindex + 1 == cols ) {
469 if ( yindex + 1 == rows ) {
473 if ( (xindex < 0) || (xindex + 1 >= cols) ||
474 (yindex < 0) || (yindex + 1 >= rows) ) {
478 dx = xlocal - xindex;
479 dy = ylocal - yindex;
483 // printf(" Lower triangle\n");
487 z1 = dem_data[x1][y1];
491 z2 = dem_data[x2][y2];
495 z3 = dem_data[x3][y3];
497 // printf(" dx = %.2f dy = %.2f\n", dx, dy);
498 // printf(" (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1);
499 // printf(" (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2);
500 // printf(" (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3);
502 zA = dx * (z2 - z1) + z1;
503 zB = dx * (z3 - z1) + z1;
505 // printf(" zA = %.2f zB = %.2f\n", zA, zB);
507 if ( dx > FG_EPSILON ) {
508 elev = dy * (zB - zA) / dx + zA;
514 // printf(" Upper triangle\n");
518 z1 = dem_data[x1][y1];
522 z2 = dem_data[x2][y2];
526 z3 = dem_data[x3][y3];
528 // printf(" dx = %.2f dy = %.2f\n", dx, dy);
529 // printf(" (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1);
530 // printf(" (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2);
531 // printf(" (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3);
533 zA = dy * (z2 - z1) + z1;
534 zB = dy * (z3 - z1) + z1;
536 // printf(" zA = %.2f zB = %.2f\n", zA, zB );
537 // printf(" xB - xA = %.2f\n", col_step * dy / row_step);
539 if ( dy > FG_EPSILON ) {
540 elev = dx * (zB - zA) / dy + zA;
550 // Use least squares to fit a simpler data set to dem data
551 void fgDEM::fit( double error, fgBUCKET *p ) {
552 double x[DEM_SIZE_1], y[DEM_SIZE_1];
553 double m, b, ave_error, max_error;
555 int n, row, start, end;
556 int colmin, colmax, rowmin, rowmax;
558 // FILE *dem, *fit, *fit1;
560 printf("Initializing output mesh structure\n");
563 // determine dimensions
564 colmin = p->x * ( (cols - 1) / 8);
565 colmax = colmin + ( (cols - 1) / 8);
566 rowmin = p->y * ( (rows - 1) / 8);
567 rowmax = rowmin + ( (rows - 1) / 8);
568 printf("Fitting region = %d,%d to %d,%d\n", colmin, rowmin, colmax, rowmax);
570 // include the corners explicitly
571 outputmesh_set_pt(colmin, rowmin, dem_data[colmin][rowmin]);
572 outputmesh_set_pt(colmin, rowmax, dem_data[colmin][rowmax]);
573 outputmesh_set_pt(colmax, rowmax, dem_data[colmax][rowmax]);
574 outputmesh_set_pt(colmax, rowmin, dem_data[colmax][rowmin]);
576 printf("Beginning best fit procedure\n");
578 for ( row = rowmin; row <= rowmax; row++ ) {
579 // fit = fopen("fit.dat", "w");
580 // fit1 = fopen("fit1.dat", "w");
584 // printf(" fitting row = %d\n", row);
586 while ( start < colmax ) {
590 x[(end - start) - 1] = 0.0 + ( start * col_step );
591 y[(end - start) - 1] = dem_data[start][row];
593 while ( (end <= colmax) && good_fit ) {
594 n = (end - start) + 1;
595 // printf("Least square of first %d points\n", n);
596 x[end - start] = 0.0 + ( end * col_step );
597 y[end - start] = dem_data[end][row];
598 least_squares(x, y, n, &m, &b);
599 ave_error = least_squares_error(x, y, n, m, b);
600 max_error = least_squares_max_error(x, y, n, m, b);
603 printf("%d - %d ave error = %.2f max error = %.2f y = %.2f*x + %.2f\n",
604 start, end, ave_error, max_error, m, b);
606 f = fopen("gnuplot.dat", "w");
607 for ( j = 0; j <= end; j++) {
608 fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ),
611 for ( j = start; j <= end; j++) {
612 fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ),
617 printf("Please hit return: "); gets(junk);
620 if ( max_error > error ) {
628 // error exceeded the threshold, back up
629 end -= 2; // back "end" up to the last good enough fit
630 n--; // back "n" up appropriately too
632 // we popped out of the above loop while still within
633 // the error threshold, so we must be at the end of
638 least_squares(x, y, n, &m, &b);
639 ave_error = least_squares_error(x, y, n, m, b);
640 max_error = least_squares_max_error(x, y, n, m, b);
644 printf("%d - %d ave error = %.2f max error = %.2f y = %.2f*x + %.2f\n",
645 start, end, ave_error, max_error, m, b);
648 fprintf(fit1, "%.2f %.2f\n", x[0], m * x[0] + b);
649 fprintf(fit1, "%.2f %.2f\n", x[end-start], m * x[end-start] + b);
652 if ( start > colmin ) {
653 // skip this for the first line segment
655 outputmesh_set_pt(start, row, (lasty + cury) / 2);
656 // fprintf(fit, "%.2f %.2f\n", x[0], (lasty + cury) / 2);
659 lasty = m * x[end-start] + b;
667 dem = fopen("gnuplot.dat", "w");
668 for ( j = 0; j < DEM_SIZE_1; j++) {
669 fprintf(dem, "%.2f %.2f\n", 0.0 + ( j * col_step ),
675 // NOTICE, this is for testing only. This instance of
676 // output_nodes should be removed. It should be called only
677 // once at the end once all the nodes have been generated.
678 // newmesh_output_nodes(&nm, "mesh.node");
679 // printf("Please hit return: "); gets(junk);
682 // outputmesh_output_nodes(fg_root, p);
686 // Initialize output mesh structure
687 void fgDEM::outputmesh_init( void ) {
690 for ( j = 0; j < DEM_SIZE_1; j++ ) {
691 for ( i = 0; i < DEM_SIZE_1; i++ ) {
692 output_data[i][j] = -9999.0;
698 // Get the value of a mesh node
699 double fgDEM::outputmesh_get_pt( int i, int j ) {
700 return ( output_data[i][j] );
704 // Set the value of a mesh node
705 void fgDEM::outputmesh_set_pt( int i, int j, double value ) {
706 // printf("Setting data[%d][%d] = %.2f\n", i, j, value);
707 output_data[i][j] = value;
711 // Write out a node file that can be used by the "triangle" program.
712 // Check for an optional "index.node.ex" file in case there is a .poly
713 // file to go along with this node file. Include these nodes first
714 // since they are referenced by position from the .poly file.
715 void fgDEM::outputmesh_output_nodes( const string& fg_root, fgBUCKET *p ) {
716 double exnodes[MAX_EX_NODES][3];
717 struct stat stat_buf;
719 char base_path[256], file[256], exfile[256];
726 int colmin, colmax, rowmin, rowmax;
727 int i, j, count, excount, result;
729 // determine dimensions
730 colmin = p->x * ( (cols - 1) / 8);
731 colmax = colmin + ( (cols - 1) / 8);
732 rowmin = p->y * ( (rows - 1) / 8);
733 rowmax = rowmin + ( (rows - 1) / 8);
734 cout << " dumping region = " << colmin << "," << rowmin << " to " <<
735 colmax << "," << rowmax << "\n";
737 // generate the base directory
738 fgBucketGenBasePath(p, base_path);
739 cout << "fg_root = " << fg_root << " Base Path = " << base_path << endl;
740 dir = fg_root + "/Scenery/" + base_path;
741 cout << "Dir = " << dir << endl;
743 // stat() directory and create if needed
745 result = stat(dir.c_str(), &stat_buf);
746 if ( result != 0 && errno == ENOENT ) {
747 cout << "Creating directory\n";
751 command = "mkdir -p " + dir + "\n";
752 system( command.c_str() );
756 // Cygwin crashes when trying to output to node file
757 // explicitly making directory structure seems OK on Win95
759 extract_path (base_path, tmp_path);
761 dir = fg_root + "/Scenery";
762 if (my_mkdir ( dir.c_str() )) { exit (-1); }
764 dir = fg_root + "/Scenery/" + tmp_path;
765 if (my_mkdir ( dir.c_str() )) { exit (-1); }
767 dir = fg_root + "/Scenery/" + base_path;
768 if (my_mkdir ( dir.c_str() )) { exit (-1); }
773 // assume directory exists
776 // get index and generate output file name
777 index = fgBucketGenIndex(p);
778 sprintf(file, "%s/%ld.node", dir.c_str(), index);
780 // get (optional) extra node file name (in case there is matching
782 strcpy(exfile, file);
783 strcat(exfile, ".ex");
785 // load extra nodes if they exist
787 if ( (fd = fopen(exfile, "r")) != NULL ) {
789 fscanf(fd, "%d %d %d %d", &excount, &junki, &junki, &junki);
791 if ( excount > MAX_EX_NODES - 1 ) {
792 printf("Error, too many 'extra' nodes, increase array size\n");
795 printf(" Expecting %d 'extra' nodes\n", excount);
798 for ( i = 1; i <= excount; i++ ) {
799 fscanf(fd, "%d %lf %lf %lf\n", &junki,
800 &exnodes[i][0], &exnodes[i][1], &exnodes[i][2]);
801 printf("(extra) %d %.2f %.2f %.2f\n",
802 i, exnodes[i][0], exnodes[i][1], exnodes[i][2]);
807 printf("Creating node file: %s\n", file);
808 fd = fopen(file, "w");
810 // first count regular nodes to generate header
812 for ( j = rowmin; j <= rowmax; j++ ) {
813 for ( i = colmin; i <= colmax; i++ ) {
814 if ( output_data[i][j] > -9000.0 ) {
818 // printf(" count = %d\n", count);
820 fprintf(fd, "%d 2 1 0\n", count + excount);
822 // now write out extra node data
823 for ( i = 1; i <= excount; i++ ) {
824 fprintf(fd, "%d %.2f %.2f %.2f\n",
825 i, exnodes[i][0], exnodes[i][1], exnodes[i][2]);
828 // write out actual node data
830 for ( j = rowmin; j <= rowmax; j++ ) {
831 for ( i = colmin; i <= colmax; i++ ) {
832 if ( output_data[i][j] > -9000.0 ) {
833 fprintf(fd, "%d %.2f %.2f %.2f\n",
835 originx + (double)i * col_step,
836 originy + (double)j * row_step,
840 // printf(" count = %d\n", count);
847 fgDEM::~fgDEM( void ) {
848 // printf("class fgDEM DEstructor called.\n");
850 delete [] output_data;
855 // Revision 1.19 1998/10/22 21:59:19 curt
856 // Fixed a couple subtle bugs that resulted from some of my c++ conversions.
857 // One bug could cause a segfault on certain input, and the other bug could
858 // cause the whole procedure to go balistic and generate huge files (also only
859 // on rare input combinations.)
861 // Revision 1.18 1998/10/18 01:17:09 curt
864 // Revision 1.17 1998/10/16 19:08:12 curt
865 // Portability updates from Bernie Bright.
867 // Revision 1.16 1998/10/02 21:41:39 curt
870 // Revision 1.15 1998/09/21 20:53:59 curt
871 // minor tweaks to clean a few additional things up after the rewrite.
873 // Revision 1.14 1998/09/19 17:59:45 curt
874 // Use c++ streams (fg_gzifstream). Also converted many character arrays to
877 // Revision 1.13 1998/09/09 16:24:04 curt
878 // Fixed a bug in the handling of exclude files which was causing
879 // a crash by calling fclose() on an invalid file handle.
881 // Revision 1.12 1998/08/24 20:03:31 curt
882 // Eliminated a possible memory overrun error.
883 // Use the proper free() rather than the incorrect delete().
885 // Revision 1.11 1998/07/20 12:46:11 curt
886 // When outputing to a .node file, first check for an optional
887 // "index.node.ex" file in case there is a .poly file to go along with this
888 // node file. Include these nodes first since they are referenced by position
889 // from the .poly file. This is my first pass at adding an area "cutout"
890 // feature to the terrain generation pipeline.
892 // Revision 1.10 1998/07/13 20:58:02 curt
895 // Revision 1.9 1998/07/13 15:29:49 curt
896 // Added #ifdef HAVE_CONFIG_H
898 // Revision 1.8 1998/07/04 00:47:18 curt
899 // typedef'd struct fgBUCKET.
901 // Revision 1.7 1998/06/05 18:14:39 curt
902 // Abort out early when reading the "A" record if it doesn't look like
903 // a proper DEM file.
905 // Revision 1.6 1998/05/02 01:49:21 curt
906 // Fixed a bug where the wrong variable was being initialized.
908 // Revision 1.5 1998/04/25 15:00:32 curt
909 // Changed "r" to "rb" in gzopen() options. This fixes bad behavior in win32.
911 // Revision 1.4 1998/04/22 13:14:46 curt
912 // Fixed a bug in zlib usage.
914 // Revision 1.3 1998/04/18 03:53:05 curt
915 // Added zlib support.
917 // Revision 1.2 1998/04/14 02:43:27 curt
918 // Used "new" to auto-allocate large DEM parsing arrays in class constructor.
920 // Revision 1.1 1998/04/08 22:57:22 curt
921 // Adopted Gnu automake/autoconf system.
923 // Revision 1.3 1998/04/06 21:09:41 curt
924 // Additional win32 support.
925 // Fixed a bad bug in dem file parsing that was causing the output to be
926 // flipped about x = y.
928 // Revision 1.2 1998/03/23 20:35:41 curt
929 // Updated to use FG_EPSILON
931 // Revision 1.1 1998/03/19 02:54:47 curt
932 // Reorganized into a class lib called fgDEM.
934 // Revision 1.1 1998/03/19 01:46:28 curt