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()
37 #include <unistd.h> // stat()
41 // #include <zlib/zlib.h>
42 #include <Misc/fgstream.hxx>
43 #include <Misc/strutils.hxx>
46 #include "leastsqs.hxx"
48 #include <Include/fg_constants.h>
51 #define MAX_EX_NODES 10000
54 # define MKDIR(a) mkdir(a,S_IRWXU) // I am just guessing at this flag (NHV)
58 fgDEM::fgDEM( void ) {
59 // printf("class fgDEM CONstructor called.\n");
60 dem_data = new float[DEM_SIZE_1][DEM_SIZE_1];
61 output_data = new float[DEM_SIZE_1][DEM_SIZE_1];
67 // return the file path name ( foo/bar/file.ext = foo/bar )
68 static void extract_path (char *in, char *base) {
75 while ( (i >= 0) && (in[i] != '/') ) {
83 // Make a subdirectory
84 static int my_mkdir (char *dir) {
90 result = stat (dir, &stat_buf);
94 result = stat (dir, &stat_buf);
96 printf ("problem creating %s\n", dir);
98 printf ("%s created\n", dir);
101 printf ("%s already exists\n", dir);
111 int fgDEM::open ( const string& file ) {
112 // open input file (or read from stdin)
114 printf("Loading DEM data file: stdin\n");
116 // fd = gzdopen(STDIN_FILENO, "r");
117 printf("Not yet ported ...\n");
120 in = new fg_gzifstream( file );
122 cout << "Cannot open " + file + "\n";
125 cout << "Loading DEM data file: " + file + "\n";
133 int fgDEM::close () {
134 // the fg_gzifstream doesn't seem to have a close()
142 // return next token from input stream
143 string fgDEM::next_token() {
146 in->stream() >> token;
148 // cout << " returning " + token + "\n";
154 // return next integer from input stream
155 int fgDEM::next_int() {
158 in->stream() >> result;
164 // return next double from input stream
165 double fgDEM::next_double() {
168 in->stream() >> result;
174 // return next exponential num from input stream
175 int fgDEM::next_exp() {
181 token = next_token();
183 sscanf(token.c_str(), "%lfD%d", &mantissa, &exp);
185 // cout << " Mantissa = " << mantissa << " Exp = " << exp << "\n";
189 for ( i = 1; i <= exp; i++ ) {
192 } else if ( exp < 0 ) {
193 for ( i = -1; i >= exp; i-- ) {
198 return( (int)rint(mantissa * (double)acc) );
202 // read and parse DEM "A" record
203 int fgDEM::read_a_record() {
210 // get the name field (144 characters)
211 for ( i = 0; i < 144; i++ ) {
216 // clean off the trailing whitespace
218 cout << " Quad name field: " + name + "\n";
220 // DEM level code, 3 reflects processing by DMA
222 cout << " DEM level code = " << inum << "\n";
228 // Pattern code, 1 indicates a regular elevation pattern
230 cout << " Pattern code = " << inum << "\n";
232 // Planimetric reference system code, 0 indicates geographic
233 // coordinate system.
235 cout << " Planimetric reference code = " << inum << "\n";
239 cout << " Zone code = " << inum << "\n";
241 // Map projection parameters (ignored)
242 for ( i = 0; i < 15; i++ ) {
243 dnum = next_double();
244 // printf("%d: %f\n",i,dnum);
247 // Units code, 3 represents arc-seconds as the unit of measure for
248 // ground planimetric coordinates throughout the file.
251 cout << " Unknown (X,Y) units code = " << inum << "!\n";
255 // Units code; 2 represents meters as the unit of measure for
256 // elevation coordinates throughout the file.
259 cout << " Unknown (Z) units code = " << inum << "!\n";
263 // Number (n) of sides in the polygon which defines the coverage of
264 // the DEM file (usually equal to 4).
267 cout << " Unknown polygon dimension = " << inum << "!\n";
271 // Ground coordinates of bounding box in arc-seconds
272 dem_x1 = originx = next_exp();
273 dem_y1 = originy = next_exp();
274 cout << " Origin = (" << originx << "," << originy << ")\n";
285 // Minimum/maximum elevations in meters
288 cout << " Elevation range " << dem_z1 << " to " << dem_z2 << "\n";
290 // Counterclockwise angle from the primary axis of ground
291 // planimetric referenced to the primary axis of the DEM local
293 token = next_token();
295 // Accuracy code; 0 indicates that a record of accuracy does not
296 // exist and that no record type C will follow.
298 // DEM spacial resolution. Usually (3,3,1) (3,6,1) or (3,9,1)
299 // depending on latitude
301 // I will eventually have to do something with this for data at
302 // higher latitudes */
303 token = next_token();
304 cout << " accuracy & spacial resolution string = " + token + "\n";
306 cout << " length = " << i << "\n";
308 ptr = token.c_str() + i - 12;
309 cout << " last field = " << ptr << " = " << atof(ptr) << "\n";
313 col_step = atof(ptr);
314 cout << " last field = " << ptr << " = " << col_step << "\n";
318 row_step = atof(ptr);
319 cout << " last field = " << ptr << " = " << row_step << "\n";
322 // accuracy code = atod(token)
325 cout << " Accuracy code = " << inum << "\n";
327 cout << " column step = " << col_step <<
328 " row step = " << row_step << "\n";
330 // dimension of arrays to follow (1)
331 token = next_token();
333 // number of profiles
334 dem_num_profiles = cols = next_int();
335 cout << " Expecting " << dem_num_profiles << " profiles\n";
341 // read and parse DEM "B" record
342 void fgDEM::read_b_record( ) {
346 // row / column id of this profile
347 prof_row = next_int();
348 prof_col = next_int();
349 // printf("col id = %d row id = %d\n", prof_col, prof_row);
351 // Number of columns and rows (elevations) in this profile
352 prof_num_rows = rows = next_int();
353 prof_num_cols = next_int();
354 // printf(" profile num rows = %d\n", prof_num_rows);
356 // Ground planimetric coordinates (arc-seconds) of the first
357 // elevation in the profile
358 prof_x1 = next_exp();
359 prof_y1 = next_exp();
360 // printf(" Starting at %.2f %.2f\n", prof_x1, prof_y1);
362 // Elevation of local datum for the profile. Always zero for
363 // 1-degree DEM, the reference is mean sea level.
364 token = next_token();
366 // Minimum and maximum elevations for the profile.
367 token = next_token();
368 token = next_token();
370 // One (usually) dimensional array (prof_num_cols,1) of elevations
371 for ( i = 0; i < prof_num_rows; i++ ) {
372 prof_data = next_int();
373 dem_data[cur_col][i] = (float)prof_data;
379 int fgDEM::parse( ) {
384 if ( !read_a_record() ) {
388 for ( i = 0; i < dem_num_profiles; i++ ) {
389 // printf("Ready to read next b record\n");
393 if ( cur_col % 100 == 0 ) {
394 cout << " loaded " << cur_col << " profiles of data\n";
398 cout << " Done parsing\n";
404 // return the current altitude based on mesh data. We should rewrite
405 // this to interpolate exact values, but for now this is good enough
406 double fgDEM::interpolate_altitude( double lon, double lat ) {
407 // we expect incoming (lon,lat) to be in arcsec for now
409 double xlocal, ylocal, dx, dy, zA, zB, elev;
410 int x1, x2, x3, y1, y2, y3;
414 /* determine if we are in the lower triangle or the upper triangle
422 then calculate our end points
425 xlocal = (lon - originx) / col_step;
426 ylocal = (lat - originy) / row_step;
428 xindex = (int)(xlocal);
429 yindex = (int)(ylocal);
431 // printf("xindex = %d yindex = %d\n", xindex, yindex);
433 if ( xindex + 1 == cols ) {
437 if ( yindex + 1 == rows ) {
441 if ( (xindex < 0) || (xindex + 1 >= cols) ||
442 (yindex < 0) || (yindex + 1 >= rows) ) {
446 dx = xlocal - xindex;
447 dy = ylocal - yindex;
451 // printf(" Lower triangle\n");
455 z1 = dem_data[x1][y1];
459 z2 = dem_data[x2][y2];
463 z3 = dem_data[x3][y3];
465 // printf(" dx = %.2f dy = %.2f\n", dx, dy);
466 // printf(" (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1);
467 // printf(" (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2);
468 // printf(" (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3);
470 zA = dx * (z2 - z1) + z1;
471 zB = dx * (z3 - z1) + z1;
473 // printf(" zA = %.2f zB = %.2f\n", zA, zB);
475 if ( dx > FG_EPSILON ) {
476 elev = dy * (zB - zA) / dx + zA;
482 // printf(" Upper triangle\n");
486 z1 = dem_data[x1][y1];
490 z2 = dem_data[x2][y2];
494 z3 = dem_data[x3][y3];
496 // printf(" dx = %.2f dy = %.2f\n", dx, dy);
497 // printf(" (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1);
498 // printf(" (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2);
499 // printf(" (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3);
501 zA = dy * (z2 - z1) + z1;
502 zB = dy * (z3 - z1) + z1;
504 // printf(" zA = %.2f zB = %.2f\n", zA, zB );
505 // printf(" xB - xA = %.2f\n", col_step * dy / row_step);
507 if ( dy > FG_EPSILON ) {
508 elev = dx * (zB - zA) / dy + zA;
518 // Use least squares to fit a simpler data set to dem data
519 void fgDEM::fit( double error, fgBUCKET *p ) {
520 double x[DEM_SIZE_1], y[DEM_SIZE_1];
521 double m, b, ave_error, max_error;
523 int n, row, start, end;
524 int colmin, colmax, rowmin, rowmax;
526 // FILE *dem, *fit, *fit1;
528 printf("Initializing output mesh structure\n");
531 // determine dimensions
532 colmin = p->x * ( (cols - 1) / 8);
533 colmax = colmin + ( (cols - 1) / 8);
534 rowmin = p->y * ( (rows - 1) / 8);
535 rowmax = rowmin + ( (rows - 1) / 8);
536 printf("Fitting region = %d,%d to %d,%d\n", colmin, rowmin, colmax, rowmax);
538 // include the corners explicitly
539 outputmesh_set_pt(colmin, rowmin, dem_data[colmin][rowmin]);
540 outputmesh_set_pt(colmin, rowmax, dem_data[colmin][rowmax]);
541 outputmesh_set_pt(colmax, rowmax, dem_data[colmax][rowmax]);
542 outputmesh_set_pt(colmax, rowmin, dem_data[colmax][rowmin]);
544 printf("Beginning best fit procedure\n");
546 for ( row = rowmin; row <= rowmax; row++ ) {
547 // fit = fopen("fit.dat", "w");
548 // fit1 = fopen("fit1.dat", "w");
552 // printf(" fitting row = %d\n", row);
554 while ( start < colmax ) {
558 x[(end - start) - 1] = 0.0 + ( start * col_step );
559 y[(end - start) - 1] = dem_data[start][row];
561 while ( (end <= colmax) && good_fit ) {
562 n = (end - start) + 1;
563 // printf("Least square of first %d points\n", n);
564 x[end - start] = 0.0 + ( end * col_step );
565 y[end - start] = dem_data[end][row];
566 least_squares(x, y, n, &m, &b);
567 ave_error = least_squares_error(x, y, n, m, b);
568 max_error = least_squares_max_error(x, y, n, m, b);
571 printf("%d - %d ave error = %.2f max error = %.2f y = %.2f*x + %.2f\n",
572 start, end, ave_error, max_error, m, b);
574 f = fopen("gnuplot.dat", "w");
575 for ( j = 0; j <= end; j++) {
576 fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ),
579 for ( j = start; j <= end; j++) {
580 fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ),
585 printf("Please hit return: "); gets(junk);
588 if ( max_error > error ) {
596 // error exceeded the threshold, back up
597 end -= 2; // back "end" up to the last good enough fit
598 n--; // back "n" up appropriately too
600 // we popped out of the above loop while still within
601 // the error threshold, so we must be at the end of
606 least_squares(x, y, n, &m, &b);
607 ave_error = least_squares_error(x, y, n, m, b);
608 max_error = least_squares_max_error(x, y, n, m, b);
612 printf("%d - %d ave error = %.2f max error = %.2f y = %.2f*x + %.2f\n",
613 start, end, ave_error, max_error, m, b);
616 fprintf(fit1, "%.2f %.2f\n", x[0], m * x[0] + b);
617 fprintf(fit1, "%.2f %.2f\n", x[end-start], m * x[end-start] + b);
620 if ( start > colmin ) {
621 // skip this for the first line segment
623 outputmesh_set_pt(start, row, (lasty + cury) / 2);
624 // fprintf(fit, "%.2f %.2f\n", x[0], (lasty + cury) / 2);
627 lasty = m * x[end-start] + b;
635 dem = fopen("gnuplot.dat", "w");
636 for ( j = 0; j < DEM_SIZE_1; j++) {
637 fprintf(dem, "%.2f %.2f\n", 0.0 + ( j * col_step ),
643 // NOTICE, this is for testing only. This instance of
644 // output_nodes should be removed. It should be called only
645 // once at the end once all the nodes have been generated.
646 // newmesh_output_nodes(&nm, "mesh.node");
647 // printf("Please hit return: "); gets(junk);
650 // outputmesh_output_nodes(fg_root, p);
654 // Initialize output mesh structure
655 void fgDEM::outputmesh_init( void ) {
658 for ( j = 0; j < DEM_SIZE_1; j++ ) {
659 for ( i = 0; i < DEM_SIZE_1; i++ ) {
660 output_data[i][j] = -9999.0;
666 // Get the value of a mesh node
667 double fgDEM::outputmesh_get_pt( int i, int j ) {
668 return ( output_data[i][j] );
672 // Set the value of a mesh node
673 void fgDEM::outputmesh_set_pt( int i, int j, double value ) {
674 // printf("Setting data[%d][%d] = %.2f\n", i, j, value);
675 output_data[i][j] = value;
679 // Write out a node file that can be used by the "triangle" program.
680 // Check for an optional "index.node.ex" file in case there is a .poly
681 // file to go along with this node file. Include these nodes first
682 // since they are referenced by position from the .poly file.
683 void fgDEM::outputmesh_output_nodes( const string& fg_root, fgBUCKET *p ) {
684 double exnodes[MAX_EX_NODES][3];
685 double junk1, junk2, junk3;
686 struct stat stat_buf;
688 char base_path[256], file[256], exfile[256];
695 int colmin, colmax, rowmin, rowmax;
696 int i, j, count, excount, result;
698 // determine dimensions
699 colmin = p->x * ( (cols - 1) / 8);
700 colmax = colmin + ( (cols - 1) / 8);
701 rowmin = p->y * ( (rows - 1) / 8);
702 rowmax = rowmin + ( (rows - 1) / 8);
703 cout << " dumping region = " << colmin << "," << rowmin << " to " <<
704 colmax << "," << rowmax << "\n";
706 // generate the base directory
707 fgBucketGenBasePath(p, base_path);
708 cout << "fg_root = " + fg_root + " Base Path = " + base_path + "\n";
709 dir = fg_root + "/Scenery/" + base_path;
710 cout << "Dir = " + dir + "\n";
712 // stat() directory and create if needed
713 result = stat(dir.c_str(), &stat_buf);
715 cout << "Stat error need to create directory\n";
719 command = "mkdir -p " + dir + "\n";
720 system( command.c_str() );
724 // Cygwin crashes when trying to output to node file
725 // explicitly making directory structure seems OK on Win95
727 extract_path (base_path, tmp_path);
729 dir = fg_root + "/Scenery";
730 if (my_mkdir ( dir.c_str() )) { exit (-1); }
732 dir += fg_root + "/Scenery/" + tmp_path;
733 if (my_mkdir ( dir.c_str() )) { exit (-1); }
735 dir += fg_root + "/Scenery/" + base_path;
736 if (my_mkdir ( dir.c_str() )) { exit (-1); }
741 // assume directory exists
744 // get index and generate output file name
745 index = fgBucketGenIndex(p);
746 sprintf(file, "%s/%ld.node", dir.c_str(), index);
748 // get (optional) extra node file name (in case there is matching
750 strcpy(exfile, file);
751 strcat(exfile, ".ex");
753 // load extra nodes if they exist
755 if ( (fd = fopen(exfile, "r")) != NULL ) {
756 fscanf(fd, "%d %d %d %d", &excount, &junk1, &junk2, &junk3);
758 if ( excount > MAX_EX_NODES - 1 ) {
759 printf("Error, too many 'extra' nodes, increase array size\n");
762 printf(" Expecting %d 'extra' nodes\n", excount);
765 for ( i = 1; i <= excount; i++ ) {
766 fscanf(fd, "%d %lf %lf %lf\n", &junk1,
767 &exnodes[i][0], &exnodes[i][1], &exnodes[i][2]);
768 printf("(extra) %d %.2f %.2f %.2f\n",
769 i, exnodes[i][0], exnodes[i][1], exnodes[i][2]);
774 printf("Creating node file: %s\n", file);
775 fd = fopen(file, "w");
777 // first count regular nodes to generate header
779 for ( j = rowmin; j <= rowmax; j++ ) {
780 for ( i = colmin; i <= colmax; i++ ) {
781 if ( output_data[i][j] > -9000.0 ) {
785 // printf(" count = %d\n", count);
787 fprintf(fd, "%d 2 1 0\n", count + excount);
789 // now write out extra node data
790 for ( i = 1; i <= excount; i++ ) {
791 fprintf(fd, "%d %.2f %.2f %.2f\n",
792 i, exnodes[i][0], exnodes[i][1], exnodes[i][2]);
795 // write out actual node data
797 for ( j = rowmin; j <= rowmax; j++ ) {
798 for ( i = colmin; i <= colmax; i++ ) {
799 if ( output_data[i][j] > -9000.0 ) {
800 fprintf(fd, "%d %.2f %.2f %.2f\n",
802 originx + (double)i * col_step,
803 originy + (double)j * row_step,
807 // printf(" count = %d\n", count);
814 fgDEM::~fgDEM( void ) {
815 // printf("class fgDEM DEstructor called.\n");
822 // Revision 1.15 1998/09/21 20:53:59 curt
823 // minor tweaks to clean a few additional things up after the rewrite.
825 // Revision 1.14 1998/09/19 17:59:45 curt
826 // Use c++ streams (fg_gzifstream). Also converted many character arrays to
829 // Revision 1.13 1998/09/09 16:24:04 curt
830 // Fixed a bug in the handling of exclude files which was causing
831 // a crash by calling fclose() on an invalid file handle.
833 // Revision 1.12 1998/08/24 20:03:31 curt
834 // Eliminated a possible memory overrun error.
835 // Use the proper free() rather than the incorrect delete().
837 // Revision 1.11 1998/07/20 12:46:11 curt
838 // When outputing to a .node file, first check for an optional
839 // "index.node.ex" file in case there is a .poly file to go along with this
840 // node file. Include these nodes first since they are referenced by position
841 // from the .poly file. This is my first pass at adding an area "cutout"
842 // feature to the terrain generation pipeline.
844 // Revision 1.10 1998/07/13 20:58:02 curt
847 // Revision 1.9 1998/07/13 15:29:49 curt
848 // Added #ifdef HAVE_CONFIG_H
850 // Revision 1.8 1998/07/04 00:47:18 curt
851 // typedef'd struct fgBUCKET.
853 // Revision 1.7 1998/06/05 18:14:39 curt
854 // Abort out early when reading the "A" record if it doesn't look like
855 // a proper DEM file.
857 // Revision 1.6 1998/05/02 01:49:21 curt
858 // Fixed a bug where the wrong variable was being initialized.
860 // Revision 1.5 1998/04/25 15:00:32 curt
861 // Changed "r" to "rb" in gzopen() options. This fixes bad behavior in win32.
863 // Revision 1.4 1998/04/22 13:14:46 curt
864 // Fixed a bug in zlib usage.
866 // Revision 1.3 1998/04/18 03:53:05 curt
867 // Added zlib support.
869 // Revision 1.2 1998/04/14 02:43:27 curt
870 // Used "new" to auto-allocate large DEM parsing arrays in class constructor.
872 // Revision 1.1 1998/04/08 22:57:22 curt
873 // Adopted Gnu automake/autoconf system.
875 // Revision 1.3 1998/04/06 21:09:41 curt
876 // Additional win32 support.
877 // Fixed a bad bug in dem file parsing that was causing the output to be
878 // flipped about x = y.
880 // Revision 1.2 1998/03/23 20:35:41 curt
881 // Updated to use FG_EPSILON
883 // Revision 1.1 1998/03/19 02:54:47 curt
884 // Reorganized into a class lib called fgDEM.
886 // Revision 1.1 1998/03/19 01:46:28 curt