From 8509dc1c33f29f69fde08aa46f1b98facf2686cb Mon Sep 17 00:00:00 2001 From: "James.Hester" Date: Sat, 15 Sep 2018 19:39:02 +1000 Subject: [PATCH 1/8] Added processing of DTED files. --- src/Lib/HGT/CMakeLists.txt | 1 + src/Lib/HGT/dted.cxx | 234 ++++++++++++++++++++++++++++++++ src/Lib/HGT/dted.hxx | 84 ++++++++++++ src/Prep/DemChop/CMakeLists.txt | 11 ++ src/Prep/DemChop/dtedchop.cxx | 110 +++++++++++++++ 5 files changed, 440 insertions(+) create mode 100644 src/Lib/HGT/dted.cxx create mode 100644 src/Lib/HGT/dted.hxx create mode 100644 src/Prep/DemChop/dtedchop.cxx diff --git a/src/Lib/HGT/CMakeLists.txt b/src/Lib/HGT/CMakeLists.txt index b120f7fa..9e1d239f 100644 --- a/src/Lib/HGT/CMakeLists.txt +++ b/src/Lib/HGT/CMakeLists.txt @@ -1,4 +1,5 @@ add_library(HGT STATIC + dted.cxx dted.hxx hgt.cxx hgt.hxx srtmbase.cxx srtmbase.hxx ) diff --git a/src/Lib/HGT/dted.cxx b/src/Lib/HGT/dted.cxx new file mode 100644 index 00000000..05d7de53 --- /dev/null +++ b/src/Lib/HGT/dted.cxx @@ -0,0 +1,234 @@ +// dted.cxx -- SRTM "dted" data management class +// +// Written by James Hester based on hgt code of +// Curtis Olson, started February 2003. +// +// Copyright (C) 2003 Curtis L. Olson - http://www.flightgear.org/~curt +// Copyright (C) 2018 James Hester +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. +// +// $Id: hgt.cxx,v 1.7 2005-12-19 16:06:45 curt Exp $ + + +#ifdef HAVE_CONFIG_H +# include +#endif + +#include + +#include // atof() +#include + +#ifdef SG_HAVE_STD_INCLUDES +# include +#else +# include +#endif + +#ifdef _MSC_VER +# include +#endif + +#include +#include +#include +#include + + +#include "dted.hxx" + +using std::cout; +using std::endl; +using std::string; + + +TGDted::TGDted( int _res ) +{ + dted_resolution = _res; + + data = new short int[MAX_DTED_SIZE][MAX_DTED_SIZE]; + output_data = new short int[MAX_DTED_SIZE][MAX_DTED_SIZE]; +} + + +TGDted::TGDted( int _res, const SGPath &file ) +{ + dted_resolution = _res; + data = new short int[MAX_DTED_SIZE][MAX_DTED_SIZE]; + output_data = new short int[MAX_DTED_SIZE][MAX_DTED_SIZE]; + + TGDted::open( file ); +} + + +// open an DTED file +bool +TGDted::open ( const SGPath &f ) { + SGPath file_name = f; + + // open input file (or read from stdin) + if ( file_name.str() == "-" ) { + cout << "Loading DTED data file: stdin" << endl; + if ( (fd = gzdopen(0, "r")) == NULL ) { // 0 == STDIN_FILENO + cout << "ERROR: opening stdin" << endl; + return false; + } + } else { + if ( file_name.extension() == "zip" ) { + // extract the .zip file to /tmp and point the file name + // to the extracted file + tmp_dir = simgear::Dir::tempDir("dted"); + + cout << "Extracting " << file_name.str() << " to " << tmp_dir.path().str() << endl; + string command = "unzip -d \"" + tmp_dir.path().str() + "\" " + file_name.base(); + if ( system( command.c_str() ) != -1 ) + { + simgear::PathList files = tmp_dir.children(simgear::Dir::TYPE_FILE | simgear::Dir::NO_DOT_OR_DOTDOT); + for (const SGPath& file : files) { + string ext = file.lower_extension(); + if ( ext == "dted" ) { + file_name = file; + break; + } + } + + remove_tmp_file = true; + cout << "Proceeding with " << file_name.str() << endl; + } else { + SG_LOG(SG_GENERAL, SG_ALERT, "Failed to issue system call " << command ); + exit(1); + } + } + + cout << "Loading DTED data file: " << file_name.str() << endl; + if ( (fd = gzopen( file_name.c_str(), "rb" )) == NULL ) { + SGPath file_name_gz = file_name; + file_name_gz.append( ".gz" ); + if ( (fd = gzopen( file_name_gz.c_str(), "rb" )) == NULL ) { + cout << "ERROR: opening " << file_name.str() << " or " + << file_name_gz.str() << " for reading!" << endl; + return false; + } + } + } + + // Determine originx/originy from file contents + // User Header Label + char header[3]; + char location[7]; + int degrees,minutes,seconds; + char hemisphere; + // Check header + gzread(fd,header,3); + if (strncmp(header,"UHL",3) != 0) { + cout << "UHL User Header Label not found" << endl; + return false; + } + gzread(fd,header,1); //dummy + gzread(fd,location,8);//longitude + sscanf(location,"%3d%2d%2d%c",°rees,&minutes,&seconds,&hemisphere); + originx = degrees *3600 + minutes*60 + seconds; + if(hemisphere == 'W') { + originx = - originx; + } + gzread(fd,location,8);//latitude + sscanf(location,"%3d%2d%2d%c",°rees,&minutes,&seconds,&hemisphere); + originy = degrees *3600 + minutes*60 + seconds; + if(hemisphere == 'S') { + originy = - originy; + } + cout << " Origin = " << originx << ", " << originy << endl; + + return true; +} + + +// close an DTED file +bool +TGDted::close () { + gzclose(fd); + return true; +} + + +// load an DTED file +bool +TGDted::load( ) { + int size; + if ( dted_resolution == 1 ) { + cols = rows = size = 3601; + col_step = row_step = 1; + } else if ( dted_resolution == 3 ) { + cols = rows = size = 1201; + col_step = row_step = 3; + } else { + cout << "Unknown DTED resolution, only 1 and 3 arcsec formats" << endl; + cout << " are supported!" << endl; + return false; + } + if (sgIsLittleEndian()) { + cout << "Little Endian: swapping input values" << endl; + } + + //Skip headers + gzseek(fd,3428,SEEK_SET); + unsigned short int latct, longct; + short int *var; + int dummy; //to read column header + for ( int col = 0; col < size; ++col ) { + dummy = 0; // zero out all bytes + longct = 0; + latct = 0; + gzread(fd,&dummy,1); //sentinel + if(dummy != 170) { + cout << "Failed to find sentinel at col " << col << endl; + return false; + } + gzread(fd,&dummy,3); //block count + gzread(fd,&longct,2); //Longitude count + gzread(fd,&latct,2); //Latitude count + if ( sgIsLittleEndian() ) { + sgEndianSwap(&longct); + } + // cout << "Longitude count " << longct << endl; + for ( int row = 0; row < size; ++row ) { + var = &data[col][row]; + if ( gzread ( fd, var, 2 ) != sizeof(short) ) { + return false; + } + if ( sgIsLittleEndian() ) { + sgEndianSwap( (unsigned short int*)var); + } + } + gzread(fd,&dummy,4); //Checksum + // Check values are right + if (col == 0) { + cout << data[col][0] << endl; + cout << data[col][1] << endl; + cout << data[col][2] << endl; + } + } + + return true; +} + + + +TGDted::~TGDted() { + // printf("class TGSrtmBase DEstructor called.\n"); + delete [] data; + delete [] output_data; +} diff --git a/src/Lib/HGT/dted.hxx b/src/Lib/HGT/dted.hxx new file mode 100644 index 00000000..fbd2d426 --- /dev/null +++ b/src/Lib/HGT/dted.hxx @@ -0,0 +1,84 @@ +// dted.hxx -- SRTM "dted" data management class +// +// Written by James Hester based on hgt.hxx, which was +// written by Curtis Olson who started February 2003. +// +// Copyright (C) 2003 Curtis L. Olson - http://www.flightgear.org/~curt +// Copyright (c) 2018 James Hester +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. +// + + +#ifndef _DTED_HXX +#define _DTED_HXX + +#ifdef HAVE_CONFIG_H +# include +#endif + +#include + +#include "srtmbase.hxx" + +#include + +#include + +#include +#include + +#define MAX_DTED_SIZE 3601 + + +class TGDted : public TGSrtmBase { + +private: + + // file pointer for input + gzFile fd; + + int dted_resolution; + + // pointers to the actual grid data allocated here + short int (*data)[MAX_DTED_SIZE]; + short int (*output_data)[MAX_DTED_SIZE]; + +public: + + // Constructor, _res must be either "1" for the 1arcsec data or + // "3" for the 3arcsec data. + TGDted( int _res ); + TGDted( int _res, const SGPath &file ); + + // Destructor + ~TGDted(); + + // open an DTED file (use "-" if input is coming from stdin) + bool open ( const SGPath &file ); + + // close an DTED file + bool close(); + + // load an dted file + bool load(); + + virtual short height( int x, int y ) const { return data[x][y]; } +}; + + +#endif // _DTED_HXX + + diff --git a/src/Prep/DemChop/CMakeLists.txt b/src/Prep/DemChop/CMakeLists.txt index 2f2a041b..00d0ea37 100644 --- a/src/Prep/DemChop/CMakeLists.txt +++ b/src/Prep/DemChop/CMakeLists.txt @@ -19,6 +19,17 @@ target_link_libraries(hgtchop install(TARGETS hgtchop RUNTIME DESTINATION bin) +add_executable(dtedchop dtedchop.cxx) + +target_link_libraries(dtedchop + HGT + ${ZLIB_LIBRARY} + ${SIMGEAR_CORE_LIBRARIES} + ${SIMGEAR_CORE_LIBRARY_DEPENDENCIES}) + +install(TARGETS dtedchop RUNTIME DESTINATION bin) + + if(TIFF_FOUND) if(MSVC AND CMAKE_CL_64) set( SRTMCHOP_LIBRARIES ${JPEG_LIBRARY} ) diff --git a/src/Prep/DemChop/dtedchop.cxx b/src/Prep/DemChop/dtedchop.cxx new file mode 100644 index 00000000..e629c494 --- /dev/null +++ b/src/Prep/DemChop/dtedchop.cxx @@ -0,0 +1,110 @@ +// dtedchop.cxx -- chop up a dted file into it's corresponding pieces and stuff +// them into the workspace directory +// +// Adapted by James Hester from hgtchop +// Written by Curtis Olson, started March 1999. +// +// Copyright (C) 1997 Curtis L. Olson - http://www.flightgear.org/~curt +// Copyright (C) 2018 James Hester +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. +// + +#ifdef HAVE_CONFIG_H +# include +#endif + +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +using std::cout; +using std::endl; +using std::string; + + +int main(int argc, char **argv) { + sglog().setLogLevels( SG_ALL, SG_WARN ); + SG_LOG( SG_GENERAL, SG_ALERT, "dtedchop version " << getTGVersion() << "\n" ); + + if ( argc != 4 ) { + cout << "Usage " << argv[0] << " " << endl; + cout << endl; + cout << "\tresolution must be either 1 or 3 (1-arc-sec or 3-arc-sec)" << endl; + + return EXIT_FAILURE; + } + + int resolution = std::stoi(string(argv[1])); + string dted_name = string(argv[2]); + string work_dir = string(argv[3]); + + // determine if file is 1arc-sec or 3arc-sec variety + if ( resolution != 1 && resolution != 3 ) { + cout << "ERROR: resolution must be 1 or 3." << endl; + return EXIT_FAILURE; + } + + SGPath sgp( work_dir ); + simgear::Dir workDir(sgp); + workDir.create(0755); + + TGDted dted(resolution, dted_name); + dted.load(); + dted.close(); + + SGGeod min = SGGeod::fromDeg( dted.get_originx() / 3600.0 + SG_HALF_BUCKET_SPAN, + dted.get_originy() / 3600.0 + SG_HALF_BUCKET_SPAN ); + SGGeod max = SGGeod::fromDeg( (dted.get_originx() + dted.get_cols() * dted.get_col_step()) / 3600.0 - SG_HALF_BUCKET_SPAN, + (dted.get_originy() + dted.get_rows() * dted.get_row_step()) / 3600.0 - SG_HALF_BUCKET_SPAN ); + SGBucket b_min( min ); + SGBucket b_max( max ); + + if ( b_min == b_max ) { + dted.write_area( work_dir, b_min ); + } else { + SGBucket b_cur; + + int dx, dy; + sgBucketDiff(b_min, b_max, &dx, &dy); + + cout << "DTED file spans tile boundaries (ok)" << endl; + cout << " dx = " << dx << " dy = " << dy << endl; + + if ( (dx > 20) || (dy > 20) ) { + cout << "somethings really wrong!!!!" << endl; + return EXIT_FAILURE; + } + + for ( int j = 0; j <= dy; ++j ) { + for ( int i = 0; i <= dx; ++i ) { + b_cur = b_min.sibling(i, j); + dted.write_area( work_dir, b_cur ); + } + } + } + + return EXIT_SUCCESS; +} From 005fd6886a83f6449a2eac681194a521bd84751f Mon Sep 17 00:00:00 2001 From: "James.Hester" Date: Sat, 6 Oct 2018 18:46:04 +1000 Subject: [PATCH 2/8] Added cliff calculation: 1. Added tg_contour includes to cmake files 2. Added calculations based on cliff contours to array.cxx --- src/BuildTiles/CMakeLists.txt | 2 + src/BuildTiles/Main/CMakeLists.txt | 1 + src/BuildTiles/Main/default_priorities.txt | 1 + src/BuildTiles/Main/priorities.hxx | 10 +- src/BuildTiles/Main/tgconstruct_elevation.cxx | 2 +- src/Lib/Array/CMakeLists.txt | 4 + src/Lib/Array/array.cxx | 295 +++++++++++++++--- src/Lib/Array/array.hxx | 12 + src/Lib/Array/testarray.cxx | 4 + src/Lib/terragear/CMakeLists.txt | 3 +- src/Lib/terragear/tg_contour.cxx | 45 +++ src/Lib/terragear/tg_contour.hxx | 5 +- src/Prep/CMakeLists.txt | 1 + 13 files changed, 338 insertions(+), 47 deletions(-) diff --git a/src/BuildTiles/CMakeLists.txt b/src/BuildTiles/CMakeLists.txt index f071e5a5..3442ec32 100644 --- a/src/BuildTiles/CMakeLists.txt +++ b/src/BuildTiles/CMakeLists.txt @@ -1,5 +1,7 @@ include_directories(${PROJECT_SOURCE_DIR}/src/Lib) include_directories(${PROJECT_SOURCE_DIR}/src/BuildTiles) +include_directories(${PROJECT_SOURCE_DIR}/src/Lib/Array) +include_directories(${PROJECT_SOURCE_DIR}/src/Lib/terragear) #add_subdirectory(Parallel) add_subdirectory(Main) diff --git a/src/BuildTiles/Main/CMakeLists.txt b/src/BuildTiles/Main/CMakeLists.txt index fefcb625..699b9c26 100644 --- a/src/BuildTiles/Main/CMakeLists.txt +++ b/src/BuildTiles/Main/CMakeLists.txt @@ -1,4 +1,5 @@ include_directories(${GDAL_INCLUDE_DIR}) +include_directories(${PROJECT_SOURCE_DIR}/src/Lib/terragear) add_executable(tg-construct tgconstruct.hxx diff --git a/src/BuildTiles/Main/default_priorities.txt b/src/BuildTiles/Main/default_priorities.txt index bf718752..65142186 100644 --- a/src/BuildTiles/Main/default_priorities.txt +++ b/src/BuildTiles/Main/default_priorities.txt @@ -30,6 +30,7 @@ River stream IntermittentStream stream Watercourse stream Canal stream +Cliff cliff # A cliff face Glacier other # Solid ice/snow PackIce other # Water with ice packs PolarIce other diff --git a/src/BuildTiles/Main/priorities.hxx b/src/BuildTiles/Main/priorities.hxx index cd0d2cb6..f1247efa 100644 --- a/src/BuildTiles/Main/priorities.hxx +++ b/src/BuildTiles/Main/priorities.hxx @@ -140,6 +140,14 @@ public: } } + bool is_cliff_area( unsigned int p ) const { + if (area_defs[p].GetCategory() == "cliff" ) { + return true; + } else { + return false; + } + } + std::string const& get_area_name( unsigned int p ) const { return area_defs[p].GetName(); } @@ -172,4 +180,4 @@ private: unsigned int sliver_area_priority; }; -#endif // _PRIORITIES_HXX \ No newline at end of file +#endif // _PRIORITIES_HXX diff --git a/src/BuildTiles/Main/tgconstruct_elevation.cxx b/src/BuildTiles/Main/tgconstruct_elevation.cxx index 99be7bff..1c1c2773 100644 --- a/src/BuildTiles/Main/tgconstruct_elevation.cxx +++ b/src/BuildTiles/Main/tgconstruct_elevation.cxx @@ -222,4 +222,4 @@ void TGConstruct::CalcElevations( void ) } } } -} \ No newline at end of file +} diff --git a/src/Lib/Array/CMakeLists.txt b/src/Lib/Array/CMakeLists.txt index 869a607f..989241aa 100644 --- a/src/Lib/Array/CMakeLists.txt +++ b/src/Lib/Array/CMakeLists.txt @@ -3,6 +3,9 @@ add_library(Array STATIC array.cxx array.hxx ) + +include_directories(${PROJECT_SOURCE_DIR}/src/Lib/terragear) + install( TARGETS Array RUNTIME DESTINATION bin LIBRARY DESTINATION lib @@ -14,6 +17,7 @@ add_executable(test_array testarray.cxx) target_link_libraries(test_array Array + terragear ${ZLIB_LIBRARY} ${SIMGEAR_CORE_LIBRARIES} ${SIMGEAR_CORE_LIBRARY_DEPENDENCIES}) diff --git a/src/Lib/Array/array.cxx b/src/Lib/Array/array.cxx index ee8e31b8..faf35fef 100644 --- a/src/Lib/Array/array.cxx +++ b/src/Lib/Array/array.cxx @@ -56,6 +56,7 @@ TGArray::TGArray( const string &file ): // open an Array file (and fitted file if it exists) +// Also open cliffs file if it exists bool TGArray::open( const string& file_base ) { // open array data file string array_name = file_base + ".arr.gz"; @@ -79,7 +80,8 @@ bool TGArray::open( const string& file_base ) { } else { SG_LOG(SG_GENERAL, SG_DEBUG, " Opening fitted data file: " << fitted_name ); } - + // open any cliffs data file + load_cliffs(file_base); return (array_in != NULL) ? true : false; } @@ -101,6 +103,47 @@ TGArray::close() { return true; } +//This code adapted from tgconstruct::LoadLandclassPolys +//All polys in the bucket should be contours which we load +//into our contour list. + +bool TGArray::load_cliffs(const string & height_base) +{ + //Get the directory so we can list the children + tgPolygon poly; //actually a contour but whatever... + int total_contours_read = 0; + SGPath b(height_base); + simgear::Dir d(b.dir()); + simgear::PathList files = d.children(simgear::Dir::TYPE_FILE); + BOOST_FOREACH(const SGPath& p, files) { + if (p.file_base() != b.file_base()) { + continue; + } + + string lext = p.complete_lower_extension(); + if ((lext == "arr") || (lext == "arr.gz") || (lext == "btg.gz") || + (lext == "fit") || (lext == "fit.gz") || (lext == "ind")) + { + // skipped! + } else { + gzFile fp = gzopen( p.c_str(), "rb" ); + unsigned int count; + sgReadUInt( fp, &count ); + SG_LOG( SG_GENERAL, SG_DEBUG, " Load " << count << " contours from " << p.realpath() ); + + for ( unsigned int i=0; i -9000 ) { @@ -362,6 +406,9 @@ double TGArray::altitude_from_grid( double lon, double lat ) const { then calculate our end points */ + // Store in degrees for later + double londeg = lon/3600; + double latdeg = lat/3600; xlocal = (lon - originx) / col_step; ylocal = (lat - originy) / row_step; @@ -384,70 +431,231 @@ double TGArray::altitude_from_grid( double lon, double lat ) const { return -9999; } - dx = xlocal - xindex; - dy = ylocal - yindex; + // Now check if we are on the same side of any cliffs - if ( dx > dy ) { - // lower triangle + // Collect lat,long at corners of area + // remember the missing corner if three found + // Go around the rectangle clockwise from SW corner + int corners[4][2]; + int ccnt = 0; + int missing = -1; //the missing point when 3 available + double lon1 = (originx+(xindex*col_step))/3600; + double lat1 = (originy+(yindex*row_step))/3600; + double lon2 = lon1 + col_step/3600; + double lat2 = lat1 + row_step/3600; + if (check_points(lon1,lat1,londeg,latdeg)) { + corners[ccnt][0] = xindex; + corners[ccnt][1] = yindex; + ccnt++; + } else missing = 0; + if (check_points(lon1,lat2,londeg,latdeg)) { + corners[ccnt][0] = xindex; + corners[ccnt][1] = yindex+1; + ccnt++; + } else missing = 1; + if (check_points(lon2,lat2,londeg,latdeg)) { + corners[ccnt][0] = xindex+1; + corners[ccnt][1] = yindex+1; + ccnt++; + } else missing = 2; + if (check_points(lon2,lat1,londeg,latdeg)) { + corners[ccnt][0] = xindex+1; + corners[ccnt][1] = yindex; + ccnt++; + } else missing = 3; + + switch (ccnt) { + case 3: //3 points are corners of a rectangle + // choose the points so that x2 is the right angle + // and x1-x2 is the x arm of the triangle + // dx,dy are the (positive) distances from the x1 corner + SG_LOG(SG_GENERAL, SG_DEBUG, "3 points, missing #" << missing); + dx = xlocal -xindex; + dy = ylocal -yindex; + switch (missing) { + case 0: //SW corner missing + x1 = corners[0][0]; + y1 = corners[0][1]; - x1 = xindex; - y1 = yindex; - z1 = get_array_elev(x1, y1); + x2 = corners[1][0]; + y2 = corners[1][1]; - x2 = xindex + 1; - y2 = yindex; - z2 = get_array_elev(x2, y2); + x3 = corners[2][0]; + y3 = corners[2][1]; - x3 = xindex + 1; - y3 = yindex + 1; - z3 = get_array_elev(x3, y3); + dy = 1 - dy; + break; + case 1: //NW corner missing + x1 = corners[0][0]; + y1 = corners[0][1]; - if ( z1 < -9000 || z2 < -9000 || z3 < -9000 ) { + x2 = corners[2][0]; + y2 = corners[2][1]; + + x3 = corners[1][0]; + y3 = corners[1][1]; + + break; + case 2: //NE corner missing + x1 = corners[2][0]; + y1 = corners[2][1]; + + x2 = corners[0][0]; + y2 = corners[0][1]; + + x3 = corners[1][0]; + y3 = corners[1][1]; + + dx = 1 - dx; //x1 is SE corner + break; + + case 3: //SE corner missing + x1 = corners[2][0]; + y1 = corners[2][1]; + + x2 = corners[1][0]; + y2 = corners[1][1]; + + x3 = corners[0][0]; + y3 = corners[0][1]; + + dx = 1 - dx; //x1 is NE corner + dy = 1 - dy; + break; + + } + // Now do the calcs on the triangle + // We interpolate on height along x1-x2 and + // x1 - x3. Then interpolate between these + // two points along y. + z1 = get_array_elev(x1,y1); + z2 = get_array_elev(x2,y2); + z3 = get_array_elev(x3,y3); + zA = dx * (z2 - z1) + z1; + zB = dx * (z3 - z1) + z1; + + if ( dx > SG_EPSILON ) { + elev = dy * (zB - zA) / dx + zA; + } else { + elev = zA; + } + + break; + case 2: //project onto line connecting two points + x1 = corners[0][0]; + y1 = corners[0][1]; + z1 = get_array_elev(x1,y1); + + x2 = corners[1][0]; + y2 = corners[1][1]; + z2 = get_array_elev(x2,y2); + + //two points are either a side of the rectangle, or + //else the diagonal + dx = xlocal - x1; + dy = ylocal - y1; + if (x1==x2) { + elev = z1+dy*(z2-z1); + } + else if (y1==y2) { + elev = z1+dx*(z2-z1); + } + else { //diagonal: project onto 45 degree line + int comp1 = x2-x1; + int comp2 = y2-y1; + double dotprod = (dx*comp1 + dy*comp2)/sqrt(2); + double projlen = sqrt(dx*dx+dy*dy)*dotprod; + elev = (z2-z1)*projlen/sqrt(2); + } + break; + case 1: //only one point found + elev = get_array_elev(corners[0][0],corners[0][1]); + break; + case 0: // all points on wrong side, fall through to normal calc + SG_LOG(SG_GENERAL, SG_WARN, "All elevation grid points on wrong side of cliff for " << londeg << "," << latdeg ); + default: // all corners + dx = xlocal - xindex; + dy = ylocal - yindex; + + if ( dx > dy ) { + // lower triangle + + x1 = xindex; + y1 = yindex; + z1 = get_array_elev(x1, y1); + + x2 = xindex + 1; + y2 = yindex; + z2 = get_array_elev(x2, y2); + + x3 = xindex + 1; + y3 = yindex + 1; + z3 = get_array_elev(x3, y3); + + if ( z1 < -9000 || z2 < -9000 || z3 < -9000 ) { // don't interpolate off a void return closest_nonvoid_elev( lon, lat ); - } + } - zA = dx * (z2 - z1) + z1; - zB = dx * (z3 - z1) + z1; + zA = dx * (z2 - z1) + z1; + zB = dx * (z3 - z1) + z1; - if ( dx > SG_EPSILON ) { + if ( dx > SG_EPSILON ) { elev = dy * (zB - zA) / dx + zA; - } else { + } else { elev = zA; - } - } else { - // upper triangle + } + } else { + // upper triangle - x1 = xindex; - y1 = yindex; - z1 = get_array_elev(x1, y1); + x1 = xindex; + y1 = yindex; + z1 = get_array_elev(x1, y1); - x2 = xindex; - y2 = yindex + 1; - z2 = get_array_elev(x2, y2); + x2 = xindex; + y2 = yindex + 1; + z2 = get_array_elev(x2, y2); - x3 = xindex + 1; - y3 = yindex + 1; - z3 = get_array_elev(x3, y3); + x3 = xindex + 1; + y3 = yindex + 1; + z3 = get_array_elev(x3, y3); - if ( z1 < -9000 || z2 < -9000 || z3 < -9000 ) { + if ( z1 < -9000 || z2 < -9000 || z3 < -9000 ) { // don't interpolate off a void return closest_nonvoid_elev( lon, lat ); - } + } - zA = dy * (z2 - z1) + z1; - zB = dy * (z3 - z1) + z1; + zA = dy * (z2 - z1) + z1; + zB = dy * (z3 - z1) + z1; - if ( dy > SG_EPSILON ) { + if ( dy > SG_EPSILON ) { elev = dx * (zB - zA) / dy + zA; - } else { + } else { elev = zA; - } - } - + } + } + } return elev; } +// Check that two points are on the same side of all cliff contours +// Could speed up by checking bounding box first +bool TGArray::check_points (const double lon1, const double lat1, const double lon2, const double lat2) const { + if (cliffs_list.size()==0) return true; + if (fabs(lon1-lon2) #include #include +#include +#include + +#include "tg_contour.hxx" +#include "tg_polygon.hxx" class TGArray { @@ -52,6 +57,8 @@ private: std::vector corner_list; std::vector fitted_list; + // list of cliff contours + tgcontour_list cliffs_list; void parse_bin(); public: @@ -65,6 +72,9 @@ public: // open an Array file (use "-" if input is coming from stdin) bool open ( const std::string& file_base ); + // Load contours from polygon files delineating height discontinuities + bool load_cliffs(const std::string & height_base); + // return if array was successfully opened or not bool is_open() const; @@ -103,6 +113,8 @@ public: int get_array_elev( int col, int row ) const; void set_array_elev( int col, int row, int val ); + // Check whether or not two points are on the same side of contour + bool check_points (const double a,const double b, const double c, const double d) const; // reset Array to initial state - ready to load another elevation file void unload( void ); }; diff --git a/src/Lib/Array/testarray.cxx b/src/Lib/Array/testarray.cxx index a0e6a268..d544ef26 100644 --- a/src/Lib/Array/testarray.cxx +++ b/src/Lib/Array/testarray.cxx @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -42,6 +43,9 @@ int main( int argc, char **argv ) { double lon, lat; check_for_help(argc, argv); + + sglog().setLogLevels(SG_ALL,SG_INFO); + sglog().set_log_priority(SG_DEBUG); if ( argc != 4 ) { cout << "ERROR: Needs 3 arguments!" << endl; diff --git a/src/Lib/terragear/CMakeLists.txt b/src/Lib/terragear/CMakeLists.txt index f9c9a257..e5ac6b02 100644 --- a/src/Lib/terragear/CMakeLists.txt +++ b/src/Lib/terragear/CMakeLists.txt @@ -1,4 +1,5 @@ include_directories(${GDAL_INCLUDE_DIR}) +include_directories(${PROJECT_SOURCE_DIR}/src/Lib/terragear) include( ${CGAL_USE_FILE} ) @@ -32,4 +33,4 @@ add_library(terragear STATIC tg_unique_vec2f.hxx tg_unique_vec3d.hxx tg_unique_vec3f.hxx -) \ No newline at end of file +) diff --git a/src/Lib/terragear/tg_contour.cxx b/src/Lib/terragear/tg_contour.cxx index f5d22974..2937a5fa 100644 --- a/src/Lib/terragear/tg_contour.cxx +++ b/src/Lib/terragear/tg_contour.cxx @@ -74,6 +74,51 @@ double tgContour::GetArea( void ) const return fabs(area * 0.5); } +// Check that the two supplied points are on the same side of the contour +bool tgContour::AreSameSide( const SGGeod& firstpt, const SGGeod& secondpt) const +{ + //Find equation of line segment joining the points + double x1 = firstpt.getLatitudeDeg(); + double x2 = secondpt.getLatitudeDeg(); + double y1 = firstpt.getLongitudeDeg(); + double y2 = secondpt.getLongitudeDeg(); + //Store differences for later + double xdif = x1-x2; + double ydif = y1-y2; + /*We describe a line parametrically: + + x1 (x2-x1) + L = + t + y1 (y2-y1) + +with u the parametric coefficient for the second line. +Then the line segments intersect if 0 <= t,u <= 1. + + */ + //Now cycle over all nodes and count how many times we intersect + int intersect_ct = 0; + if (node_list.size()) { + int j = node_list.size() - 1; + for (int i=0;i= 0.0 && t <= 1.0 && u >= 0.0 && u <= 1.0) intersect_ct++; + } + } + } + bool isinter = (intersect_ct%2 == 0); + return isinter; +} + bool tgContour::IsInside( const tgContour& inside, const tgContour& outside ) { // first contour is inside second if the intersection of first with second is == first diff --git a/src/Lib/terragear/tg_contour.hxx b/src/Lib/terragear/tg_contour.hxx index def48cf4..37e2767e 100644 --- a/src/Lib/terragear/tg_contour.hxx +++ b/src/Lib/terragear/tg_contour.hxx @@ -101,6 +101,9 @@ public: } + // Return true if the two points are on the same side of the contour + bool AreSameSide(const SGGeod& firstpt, const SGGeod& secondpt) const; + static tgContour Snap( const tgContour& subject, double snap ); static tgContour RemoveDups( const tgContour& subject ); static tgContour SplitLongEdges( const tgContour& subject, double dist ); @@ -141,4 +144,4 @@ typedef std::vector tgcontour_list; typedef tgcontour_list::iterator tgcontour_list_iterator; typedef tgcontour_list::const_iterator const_tgcontour_list_iterator; -#endif // _TGCONTOUR_HXX \ No newline at end of file +#endif // _TGCONTOUR_HXX diff --git a/src/Prep/CMakeLists.txt b/src/Prep/CMakeLists.txt index f70dbfa1..d2e471ba 100644 --- a/src/Prep/CMakeLists.txt +++ b/src/Prep/CMakeLists.txt @@ -5,3 +5,4 @@ add_subdirectory(DemChop) add_subdirectory(Terra) add_subdirectory(TerraFit) add_subdirectory(OGRDecode) +add_subdirectory(Cliff) From d36116d9657dbf0c0d40a0763388eeb0fde21d6e Mon Sep 17 00:00:00 2001 From: "James.Hester" Date: Sun, 21 Oct 2018 10:07:23 +1100 Subject: [PATCH 3/8] Cliffs now appear more or less correctly. Added cliffs to default_priorities Ignore cliff files when constructing terrain Use cliff files when calculating elevations Set up chopper to allow file extension to be specified Specify cliff file name when chopping Allow chopper to chop lines as well as polygons. --- src/BuildTiles/Main/default_priorities.txt | 4 ++- src/BuildTiles/Main/tgconstruct_poly.cxx | 4 ++- src/Lib/Array/array.cxx | 12 +++---- src/Lib/terragear/tg_chopper.cxx | 14 ++++---- src/Lib/terragear/tg_chopper.hxx | 5 +++ src/Lib/terragear/tg_contour.cxx | 40 +++++++++++++--------- src/Lib/terragear/tg_contour.hxx | 11 ++++++ src/Lib/terragear/tg_polygon.hxx | 15 ++++++++ src/Lib/terragear/tg_polygon_clip.cxx | 30 ++++++++++++---- 9 files changed, 97 insertions(+), 38 deletions(-) diff --git a/src/BuildTiles/Main/default_priorities.txt b/src/BuildTiles/Main/default_priorities.txt index 65142186..779098c6 100644 --- a/src/BuildTiles/Main/default_priorities.txt +++ b/src/BuildTiles/Main/default_priorities.txt @@ -30,7 +30,7 @@ River stream IntermittentStream stream Watercourse stream Canal stream -Cliff cliff # A cliff face +Cliffs cliff # A cliff face Glacier other # Solid ice/snow PackIce other # Water with ice packs PolarIce other @@ -86,7 +86,9 @@ ShrubGrassCover other # Mixed Shrubland/Grassland SavannaCover other # Savanna Orchard other # Orchard DeciduousForest other # Deciduous Forest +DeciduousBroadCover other # Deciduous Forest EvergreenForest other # Evergreen Forest +EvergreenBroadCover other # Evergreen Forest MixedForest other # Mixed Forest RainForest other # Rain Forest BarrenCover other # Barren or Sparsely Vegetated diff --git a/src/BuildTiles/Main/tgconstruct_poly.cxx b/src/BuildTiles/Main/tgconstruct_poly.cxx index 90d52634..5ebec27a 100644 --- a/src/BuildTiles/Main/tgconstruct_poly.cxx +++ b/src/BuildTiles/Main/tgconstruct_poly.cxx @@ -67,8 +67,10 @@ int TGConstruct::LoadLandclassPolys( void ) { } string lext = p.complete_lower_extension(); + string lastext = p.extension(); if ((lext == "arr") || (lext == "arr.gz") || (lext == "btg.gz") || - (lext == "fit") || (lext == "fit.gz") || (lext == "ind")) + (lext == "fit") || (lext == "fit.gz") || (lext == "ind") || + (lastext == "cliffs")) { // skipped! } else { diff --git a/src/Lib/Array/array.cxx b/src/Lib/Array/array.cxx index faf35fef..5138beb8 100644 --- a/src/Lib/Array/array.cxx +++ b/src/Lib/Array/array.cxx @@ -24,6 +24,7 @@ #endif #include +#include //for setprecision #include #include @@ -120,12 +121,8 @@ bool TGArray::load_cliffs(const string & height_base) continue; } - string lext = p.complete_lower_extension(); - if ((lext == "arr") || (lext == "arr.gz") || (lext == "btg.gz") || - (lext == "fit") || (lext == "fit.gz") || (lext == "ind")) - { - // skipped! - } else { + string lext = p.lower_extension(); + if (lext == "cliffs") { gzFile fp = gzopen( p.c_str(), "rb" ); unsigned int count; sgReadUInt( fp, &count ); @@ -572,7 +569,8 @@ double TGArray::altitude_from_grid( double lon, double lat ) const { elev = get_array_elev(corners[0][0],corners[0][1]); break; case 0: // all points on wrong side, fall through to normal calc - SG_LOG(SG_GENERAL, SG_WARN, "All elevation grid points on wrong side of cliff for " << londeg << "," << latdeg ); + SG_LOG(SG_GENERAL, SG_WARN, "All elevation grid points on wrong side of cliff for " << std::setprecision(10) << londeg << "," << latdeg ); + SG_LOG(SG_GENERAL, SG_WARN, "Grid points ("<< std::setprecision(9) << lon1 << "," << lat1 << "),("< 0 ) { if ( subject.GetPreserve3D() ) { result.InheritElevations( subject ); @@ -63,7 +62,8 @@ void tgChopper::ClipRow( const tgPolygon& subject, const double& center_lat, con } } -void tgChopper::Add( const tgPolygon& subject, const std::string& type ) + +void tgChopper::Add( const tgPolygon& subject, const std::string& type) { // bail out immediately if polygon is empty if ( subject.Contours() == 0 ) @@ -90,7 +90,7 @@ void tgChopper::Add( const tgPolygon& subject, const std::string& type ) // We just have a single row - no need to intersect first SG_LOG( SG_GENERAL, SG_DEBUG, " UN_CLIPPED row - center lat is " << b_min.get_center_lat() ); - ClipRow( subject, b_min.get_center_lat(), type ); + ClipRow( subject, b_min.get_center_lat(), type); } else { @@ -117,7 +117,7 @@ void tgChopper::Add( const tgPolygon& subject, const std::string& type ) clip_row.AddNode( 0, SGGeod::fromDeg( 180.0, clip_top) ); clip_row.AddNode( 0, SGGeod::fromDeg(-180.0, clip_top) ); - clipped = tgPolygon::Intersect( subject, clip_row ); + clipped = tgPolygon::Intersect( subject, clip_row); if ( clipped.TotalNodes() > 0 ) { if ( subject.GetPreserve3D() ) { @@ -205,7 +205,7 @@ void tgChopper::Save( bool DebugShapefiles ) long int poly_index = GenerateIndex( path ); sprintf( poly_ext, "%ld", poly_index ); - polyfile = polyfile + "." + poly_ext; + polyfile = polyfile + "." + poly_ext + "." + extra_extension; gzFile fp; if ( (fp = gzopen( polyfile.c_str(), "wb9" )) == NULL ) { diff --git a/src/Lib/terragear/tg_chopper.hxx b/src/Lib/terragear/tg_chopper.hxx index 47546192..e01f0215 100644 --- a/src/Lib/terragear/tg_chopper.hxx +++ b/src/Lib/terragear/tg_chopper.hxx @@ -11,10 +11,14 @@ class tgChopper public: tgChopper( const std::string& path ) { root_path = path; + extra_extension = ""; } void Add( const tgPolygon& poly, const std::string& type ); void Save( bool DebugShapes ); + void Add_Extension( const std::string& extension) { + extra_extension = extension; + } private: long int GenerateIndex( std::string path ); @@ -25,4 +29,5 @@ private: std::string root_path; bucket_polys_map bp_map; SGMutex lock; + std::string extra_extension; //add at end of file name }; diff --git a/src/Lib/terragear/tg_contour.cxx b/src/Lib/terragear/tg_contour.cxx index 2937a5fa..661279de 100644 --- a/src/Lib/terragear/tg_contour.cxx +++ b/src/Lib/terragear/tg_contour.cxx @@ -111,7 +111,11 @@ Then the line segments intersect if 0 <= t,u <= 1. double crossx = x1-nx1; double crossy = y1-ny1; double t = (crossx*nydif - crossy*nxdif)/denom; double u = -1*(xdif*crossy - ydif*crossx)/denom; - if (t >= 0.0 && t <= 1.0 && u >= 0.0 && u <= 1.0) intersect_ct++; + // We consider that an intersection at the edge of the line have + // not crossed + // over, that is, they lie on the same side, so we do not + // include equality in the comparisons + if (t > 0.0 && t < 1.0 && u > 0.0 && u < 1.0) intersect_ct++; } } } @@ -848,14 +852,16 @@ tgContour tgContour::AddColinearNodes( const tgContour& subject, std::vector node_list; bool hole; + bool keep_open; //If non-closed contour, keep open }; typedef std::vector tgcontour_list; diff --git a/src/Lib/terragear/tg_polygon.hxx b/src/Lib/terragear/tg_polygon.hxx index 0329d089..1cf7cef7 100644 --- a/src/Lib/terragear/tg_polygon.hxx +++ b/src/Lib/terragear/tg_polygon.hxx @@ -243,6 +243,7 @@ class tgPolygon public: tgPolygon() { preserve3d = false; + closed = true; } ~tgPolygon() { contours.clear(); @@ -351,6 +352,18 @@ public: flag = f; } + void SetOpen( void ) { //Do not try to close the contours + closed = false; + } + + void SetClosed (void) { + closed = true; + } + + bool IsClosed( void ) const { + return closed; + } + bool GetPreserve3D( void ) const { return preserve3d; } @@ -436,6 +449,7 @@ public: // Conversions static ClipperLib::Paths ToClipper( const tgPolygon& subject ); static tgPolygon FromClipper( const ClipperLib::Paths& subject ); + static tgPolygon FromClipper( const ClipperLib::PolyTree& subject ); static void ToClipperFile( const tgPolygon& subject, const std::string& path, const std::string& filename ); // T-Junctions and segment search @@ -460,6 +474,7 @@ private: bool preserve3d; unsigned int id; // unique polygon id for debug tgTexParams tp; + bool closed; // if we treat it as a closed shape }; #endif // _POLYGON_HXX diff --git a/src/Lib/terragear/tg_polygon_clip.cxx b/src/Lib/terragear/tg_polygon_clip.cxx index 359f9ef5..ec1c8436 100644 --- a/src/Lib/terragear/tg_polygon_clip.cxx +++ b/src/Lib/terragear/tg_polygon_clip.cxx @@ -136,7 +136,8 @@ tgPolygon tgPolygon::Diff( const tgPolygon& subject, tgPolygon& clip ) return result; } -tgPolygon tgPolygon::Intersect( const tgPolygon& subject, const tgPolygon& clip ) +//Intersect will keep open paths open +tgPolygon tgPolygon::Intersect( const tgPolygon& subject, const tgPolygon& clip) { tgPolygon result; UniqueSGGeodSet all_nodes; @@ -156,17 +157,23 @@ tgPolygon tgPolygon::Intersect( const tgPolygon& subject, const tgPolygon& clip ClipperLib::Paths clipper_subject = tgPolygon::ToClipper( subject ); ClipperLib::Paths clipper_clip = tgPolygon::ToClipper( clip ); - ClipperLib::Paths clipper_result; + ClipperLib::Paths clipper_result; // for closed polygons + ClipperLib::PolyTree clipper_tree_result; //for open polygons ClipperLib::Clipper c; c.Clear(); - c.AddPaths(clipper_subject, ClipperLib::ptSubject, true); + c.AddPaths(clipper_subject, ClipperLib::ptSubject, subject.IsClosed()); c.AddPaths(clipper_clip, ClipperLib::ptClip, true); - c.Execute(ClipperLib::ctIntersection, clipper_result, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd); + if(subject.IsClosed()) { + c.Execute(ClipperLib::ctIntersection, clipper_result, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd); - result = tgPolygon::FromClipper( clipper_result ); + result = tgPolygon::FromClipper( clipper_result ); + } + else { + c.Execute(ClipperLib::ctIntersection, clipper_tree_result, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd); + result = tgPolygon::FromClipper( clipper_tree_result ); + } result = tgPolygon::AddColinearNodes( result, all_nodes ); - result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); result.SetId( subject.GetId() ); @@ -199,6 +206,17 @@ tgPolygon tgPolygon::FromClipper( const ClipperLib::Paths& subject ) return result; } +// A polytree is returned when an open polygon has been clipped +tgPolygon tgPolygon::FromClipper( const ClipperLib::PolyTree& subject ) +{ + ClipperLib::Paths path_result; + tgPolygon result; + ClipperLib::PolyTreeToPaths(subject,path_result); + result = FromClipper(path_result); + return result; +} + + ClipperLib::Path tgTriangle::ToClipper( const tgTriangle& subject ) { ClipperLib::Path contour; From 34b3bd89ffda21473ab62fd451d604abc21ea10c Mon Sep 17 00:00:00 2001 From: "James.Hester" Date: Sun, 28 Oct 2018 23:47:23 +1100 Subject: [PATCH 4/8] Height adjustments based on cliff positions works, but is slow. --- src/Lib/Array/array.cxx | 124 +++++++++++++++++++++++++++++++ src/Lib/Array/array.hxx | 7 ++ src/Lib/terragear/tg_contour.cxx | 27 ++++++- src/Lib/terragear/tg_contour.hxx | 5 +- 4 files changed, 159 insertions(+), 4 deletions(-) diff --git a/src/Lib/Array/array.cxx b/src/Lib/Array/array.cxx index 5138beb8..d5b1f1a5 100644 --- a/src/Lib/Array/array.cxx +++ b/src/Lib/Array/array.cxx @@ -206,6 +206,9 @@ TGArray::parse( SGBucket& b ) { } } + // Fix heights near cliffs + if(cliffs_list.size()>0) rectify_heights(); + return true; } @@ -381,7 +384,115 @@ double TGArray::closest_nonvoid_elev( double lon, double lat ) const { } } +std::vector TGArray::collect_bad_points() { + //Find and remember all points that are bad because they are + //too close to a cliff + std::vector bad_points; //local to avoid multi-thread issues + for(int horiz=0;horiz bad_points) const { + int grididx; + grididx = xgrid+ygrid*cols; + auto result = std::find(std::begin(bad_points),std::end(bad_points),grididx); + if (result != std::end(bad_points)) return true; + return false; + } + + +//This may collide with other threads, but as they will both be writing +//the correct height, this is harmless. +void TGArray::rectify_heights() { + double new_ht; + std::vector rectified,bad_points; + bad_points = collect_bad_points(); + while(1) { + for (auto pt : bad_points) { + int ygrid = pt/cols; + int xgrid = pt - ygrid*cols; + new_ht = rectify_point(xgrid,ygrid,bad_points); + if (new_ht > -9999) { + rectified.push_back(pt); + set_array_elev(xgrid,ygrid,(int) new_ht); + } + } + std::cout << "Rectified " << rectified.size() << " points " << std::endl; + if(rectified.size()>0) { + for(auto r : rectified) { + bad_points.erase(std::remove(std::begin(bad_points),std::end(bad_points),r)); + } + rectified.clear(); + } else { + break; + } + } +} + +/* If we have cliffs, it is possible that a grid point will be too close +to the cliff. In this case, the SRTM-3 data appears to average the height +in the region of the point, which makes the height unreliable. This routine +searches for three neighbouring points that are reliable, and form a rectangle +with the target point, and calculates the height from the plane passing +through the three known points. + +* * * + +* x * + +* * * + +TODO: Handle points on the boundaries. */ +double TGArray::rectify_point(const int xgrid, const int ygrid, const std::vector bad_points) const { + //xgrid: grid units horizontally + //ygrid: grid units vertically + //Loop over corner points, if no points available, give up + int corners[4][2]; //possible corners + int final_pts[3][2]; // rectangle corners + int pt_cnt = 0; + double centre_long, centre_lat; + double cliff_error = col_step; //Assume row step, col step the same + int original_height = get_array_elev(xgrid,ygrid); + centre_long = (originx + col_step*xgrid)/3600; + centre_lat = (originy + row_step*ygrid)/3600; + for (int horiz = -1; horiz <= 1; horiz+=2) { + double test_long = centre_long + (col_step*horiz)/3600; + for (int vert = -1; vert <= 1; vert+=2) { + double test_lat = centre_lat + (row_step*vert)/3600; + if (!is_bad_point(xgrid+horiz,ygrid+vert,bad_points) && //can trust height + check_points(test_long,test_lat,centre_long,centre_lat)) { //same side + corners[pt_cnt][0] = horiz; + corners[pt_cnt][1] = vert; + pt_cnt++; + } + } + } // end of search for corners + if (pt_cnt == 0) return -9999; // no corners found + // Find two points that form a rectangle with a corner + int pt; + for (pt = 0; pt < pt_cnt; pt++) { + if (!is_bad_point(xgrid+corners[pt][0],ygrid,bad_points) && + !is_bad_point(xgrid, ygrid+corners[pt][1],bad_points)) break; + } + if (pt == pt_cnt) return -9999; // not found + // We have three points, calculate the height + double corner = get_array_elev(xgrid+corners[pt][0],ygrid+corners[pt][1]); + double horiz = get_array_elev(xgrid,ygrid+corners[pt][1]); + double vert = get_array_elev(xgrid+corners[pt][0],ygrid); + double height = horiz + (vert - corner); + std::cout << xgrid << "," << ygrid << ": was " << original_height << " , now " << height << std::endl; + return height; +} + // return the current altitude based on grid data. // TODO: We should rewrite this to interpolate exact values, but for now this is good enough double TGArray::altitude_from_grid( double lon, double lat ) const { @@ -655,6 +766,19 @@ bool TGArray::check_points (const double lon1, const double lat1, const double l return same_side; } +//Check that a point is more than given distance from any cliff +//Could speed up by checking bounding box +bool TGArray::is_near_cliff(const double lon1, const double lat1) const { + if (cliffs_list.size()==0) return false; + SGGeod pt1 = SGGeod::fromDeg(lon1,lat1); + double mindist = 30; // 1 arcsec + for (int i=0;i collect_bad_points(); + bool is_bad_point(const int xgrid, const int ygrid, const std::vector bad_points) const; + double rectify_point(const int xgrid, const int ygrid, const std::vector bad_points) const; + bool is_near_cliff(const double lon1,const double lon2) const; public: // Constructor diff --git a/src/Lib/terragear/tg_contour.cxx b/src/Lib/terragear/tg_contour.cxx index 661279de..b4aea891 100644 --- a/src/Lib/terragear/tg_contour.cxx +++ b/src/Lib/terragear/tg_contour.cxx @@ -1,4 +1,5 @@ #include +#include #include #include @@ -60,12 +61,11 @@ double tgContour::GetArea( void ) const SGVec2d a, b; unsigned int i, j; - if ( node_list.size() ) { + if (node_list.size() ) { j = node_list.size() - 1; for (i=0; i piece = SGLineSegment(start,end); + dist = distSqr(piece,probexyz); + if (dist < mindist) mindist = dist; + } + } + return sqrt(mindist); +} + bool tgContour::IsInside( const tgContour& inside, const tgContour& outside ) { // first contour is inside second if the intersection of first with second is == first @@ -473,6 +493,7 @@ tgContour tgContour::FromClipper( const ClipperLib::Path& subject ) return result; } + tgRectangle tgContour::GetBoundingBox( void ) const { SGGeod min, max; diff --git a/src/Lib/terragear/tg_contour.hxx b/src/Lib/terragear/tg_contour.hxx index cb356aaa..4848705c 100644 --- a/src/Lib/terragear/tg_contour.hxx +++ b/src/Lib/terragear/tg_contour.hxx @@ -7,6 +7,8 @@ #include #include +#include +#include #include #include "tg_unique_geod.hxx" @@ -113,7 +115,8 @@ public: // Return true if the two points are on the same side of the contour bool AreSameSide(const SGGeod& firstpt, const SGGeod& secondpt) const; - + // Return minimum distance of point from contour + double MinDist(const SGGeod& probe) const; static tgContour Snap( const tgContour& subject, double snap ); static tgContour RemoveDups( const tgContour& subject ); static tgContour SplitLongEdges( const tgContour& subject, double dist ); From a7af59b61f47fa6d8c57a64c1e244a735b64e4c3 Mon Sep 17 00:00:00 2001 From: "James.Hester" Date: Mon, 17 Dec 2018 23:50:28 +1100 Subject: [PATCH 5/8] Fixed height rectification along boundaries and proper clipping of open contours. --- src/BuildTiles/Main/tgconstruct_poly.cxx | 2 +- src/Lib/Array/CMakeLists.txt | 11 ++- src/Lib/Array/array.cxx | 107 ++++++++++++++++++----- src/Lib/Array/array.hxx | 13 ++- src/Lib/terragear/tg_chopper.cxx | 15 +++- src/Lib/terragear/tg_contour.cxx | 20 +++-- src/Lib/terragear/tg_polygon_clip.cxx | 6 +- 7 files changed, 141 insertions(+), 33 deletions(-) diff --git a/src/BuildTiles/Main/tgconstruct_poly.cxx b/src/BuildTiles/Main/tgconstruct_poly.cxx index 5ebec27a..bfb6e76d 100644 --- a/src/BuildTiles/Main/tgconstruct_poly.cxx +++ b/src/BuildTiles/Main/tgconstruct_poly.cxx @@ -70,7 +70,7 @@ int TGConstruct::LoadLandclassPolys( void ) { string lastext = p.extension(); if ((lext == "arr") || (lext == "arr.gz") || (lext == "btg.gz") || (lext == "fit") || (lext == "fit.gz") || (lext == "ind") || - (lastext == "cliffs")) + (lastext == "cliffs") || (lext == "arr.rectified.gz") || (lext == "arr.new.gz")) { // skipped! } else { diff --git a/src/Lib/Array/CMakeLists.txt b/src/Lib/Array/CMakeLists.txt index 989241aa..2d4f870e 100644 --- a/src/Lib/Array/CMakeLists.txt +++ b/src/Lib/Array/CMakeLists.txt @@ -15,6 +15,8 @@ install( FILES array.hxx DESTINATION include/tg ) add_executable(test_array testarray.cxx) +add_executable(rectify_height rectify_height.cxx) + target_link_libraries(test_array Array terragear @@ -22,7 +24,14 @@ target_link_libraries(test_array ${SIMGEAR_CORE_LIBRARIES} ${SIMGEAR_CORE_LIBRARY_DEPENDENCIES}) -install(TARGETS test_array RUNTIME DESTINATION bin) +target_link_libraries(rectify_height + Array + terragear + ${ZLIB_LIBRARY} + ${SIMGEAR_CORE_LIBRARIES} + ${SIMGEAR_CORE_LIBRARY_DEPENDENCIES}) + +install(TARGETS test_array rectify_height RUNTIME DESTINATION bin) if (MSVC) set_target_properties( test_array PROPERTIES DEBUG_POSTFIX d ) endif () diff --git a/src/Lib/Array/array.cxx b/src/Lib/Array/array.cxx index d5b1f1a5..9c159d96 100644 --- a/src/Lib/Array/array.cxx +++ b/src/Lib/Array/array.cxx @@ -57,16 +57,27 @@ TGArray::TGArray( const string &file ): // open an Array file (and fitted file if it exists) -// Also open cliffs file if it exists +// Also open cliffs file if it exists. By default a +// rectified file is searched for, if it doesn't exist +// we load the unrectified version. bool TGArray::open( const string& file_base ) { // open array data file - string array_name = file_base + ".arr.gz"; - + string array_name = file_base + ".arr.rectified.gz"; + rectified = true; + array_in = gzopen( array_name.c_str(), "rb" ); if (array_in == NULL) { + // try unrectified + array_name = file_base + ".arr.gz"; + array_in = gzopen(array_name.c_str(), "rb"); + if (array_in == NULL) { return false; - } + } else { + rectified = false; + } + } + SG_LOG(SG_GENERAL,SG_DEBUG,"Loaded height array " << array_name); // open fitted data file string fitted_name = file_base + ".fit.gz"; fitted_in = new sg_gzifstream( fitted_name ); @@ -206,10 +217,7 @@ TGArray::parse( SGBucket& b ) { } } - // Fix heights near cliffs - if(cliffs_list.size()>0) rectify_heights(); - - return true; + return true; } void TGArray::parse_bin() @@ -240,6 +248,40 @@ void TGArray::parse_bin() sgReadShort(array_in, cols * rows, in_data); } +// Write out an array. If rectified is true, the heights have been adjusted +// for discontinuities. +void TGArray::write_bin(const string root_dir, bool rectified, SGBucket& b) { + // generate output file name + string base = b.gen_base_path(); + string path = root_dir + "/" + base; + string extension = ".arr.new.gz"; + if (rectified) extension = ".arr.rectified.gz"; + SGPath sgp( path ); + sgp.append( "dummy" ); + sgp.create_dir( 0755 ); + + string array_file = path + "/" + b.gen_index_str() + extension; + SG_LOG(SG_GENERAL, SG_DEBUG, "array_file = " << array_file ); + + // write the file + gzFile fp; + if ( (fp = gzopen( array_file.c_str(), "wb9" )) == NULL ) { + SG_LOG(SG_GENERAL, SG_ALERT, "ERROR: cannot open " << array_file << " for writing!" ); + return; + } + + int32_t header = 0x54474152; //'TGAR' + sgWriteLong(fp,header); + sgWriteInt(fp,originx); + sgWriteInt(fp,originy); + sgWriteInt(fp,cols); + sgWriteInt(fp,col_step); + sgWriteInt(fp,rows); + sgWriteInt(fp,row_step); + sgWriteShort(fp, rows*cols, in_data); + gzclose(fp); +} + // write an Array file bool TGArray::write( const string root_dir, SGBucket& b ) { // generate output file name @@ -273,7 +315,6 @@ bool TGArray::write( const string root_dir, SGBucket& b ) { return true; } - // do our best to remove voids by picking data from the nearest neighbor. void TGArray::remove_voids( ) { // need two passes to ensure that all voids are removed (unless entire @@ -384,7 +425,7 @@ double TGArray::closest_nonvoid_elev( double lon, double lat ) const { } } -std::vector TGArray::collect_bad_points() { +std::vector TGArray::collect_bad_points(const double bad_zone) { //Find and remember all points that are bad because they are //too close to a cliff std::vector bad_points; //local to avoid multi-thread issues @@ -392,7 +433,7 @@ std::vector TGArray::collect_bad_points() { double lon = (originx + col_step*horiz)/3600; for(int vert=0;vert rectified,bad_points; - bad_points = collect_bad_points(); + bad_points = collect_bad_points(bad_zone); while(1) { for (auto pt : bad_points) { int ygrid = pt/cols; @@ -433,7 +474,10 @@ void TGArray::rectify_heights() { } rectified.clear(); } else { - break; + if(bad_points.size() > 0) { + std::cout << "Failed to rectify " << bad_points.size() << " points" << std::endl; + } + break; // Cant do any more } } } @@ -465,8 +509,10 @@ double TGArray::rectify_point(const int xgrid, const int ygrid, const std::vecto centre_long = (originx + col_step*xgrid)/3600; centre_lat = (originy + row_step*ygrid)/3600; for (int horiz = -1; horiz <= 1; horiz+=2) { + if (xgrid + horiz >= cols || xgrid + horiz < 0) continue; //edge of bucket double test_long = centre_long + (col_step*horiz)/3600; for (int vert = -1; vert <= 1; vert+=2) { + if (ygrid + vert >= rows || ygrid + vert < 0) continue; //edge of bucket double test_lat = centre_lat + (row_step*vert)/3600; if (!is_bad_point(xgrid+horiz,ygrid+vert,bad_points) && //can trust height check_points(test_long,test_lat,centre_long,centre_lat)) { //same side @@ -479,16 +525,36 @@ double TGArray::rectify_point(const int xgrid, const int ygrid, const std::vecto if (pt_cnt == 0) return -9999; // no corners found // Find two points that form a rectangle with a corner int pt; + double height = 0; for (pt = 0; pt < pt_cnt; pt++) { if (!is_bad_point(xgrid+corners[pt][0],ygrid,bad_points) && - !is_bad_point(xgrid, ygrid+corners[pt][1],bad_points)) break; - } - if (pt == pt_cnt) return -9999; // not found + !is_bad_point(xgrid, ygrid+corners[pt][1],bad_points)) { + double test_horiz = centre_long + corners[pt][0]*col_step/3600; + double test_vert = centre_lat + corners[pt][1]*row_step/3600; + if (check_points(test_horiz,centre_lat,centre_long,centre_lat) && + check_points(centre_long,test_vert,centre_long,centre_lat)) break; + } + } + + if (pt == pt_cnt) { // perhaps we are in a bay, just take the + // average of the known points + double totht = 0; + for(int pti = 0; pti collect_bad_points(); + std::vector collect_bad_points(const double bad_zone); bool is_bad_point(const int xgrid, const int ygrid, const std::vector bad_points) const; double rectify_point(const int xgrid, const int ygrid, const std::vector bad_points) const; - bool is_near_cliff(const double lon1,const double lon2) const; + bool is_near_cliff(const double lon1,const double lon2, const double bad_zone) const; public: // Constructor @@ -94,10 +95,16 @@ public: // write an Array file bool write( const std::string root_dir, SGBucket& b ); + // write an Array file in binary format. If ht_rect is true, + // the file will have extension 'arr.rectified.gz' + void write_bin(const std::string root_dir, bool ht_rect, SGBucket& b); + // do our best to remove voids by picking data from the nearest // neighbor. void remove_voids(); + void rectify_heights(const double bad_zone); + // Return the elevation of the closest non-void grid point to lon, lat double closest_nonvoid_elev( double lon, double lat ) const; diff --git a/src/Lib/terragear/tg_chopper.cxx b/src/Lib/terragear/tg_chopper.cxx index f9a52e9f..25cc71e4 100644 --- a/src/Lib/terragear/tg_chopper.cxx +++ b/src/Lib/terragear/tg_chopper.cxx @@ -21,8 +21,18 @@ tgPolygon tgChopper::Clip( const tgPolygon& subject, base.AddNode( 0, b.get_corner( SG_BUCKET_NE ) ); base.AddNode( 0, b.get_corner( SG_BUCKET_NW ) ); - result = tgPolygon::Intersect( subject, base); + result = tgPolygon::Intersect( subject, base); + // Debug: See if numbers of nodes have changed if ( result.Contours() > 0 ) { + /*if (result.ContourSize(0) != subject.ContourSize(0)) { + tgRectangle rr = subject.GetBoundingBox(); + SG_LOG(SG_GENERAL,SG_INFO,"---- Bucket ID " << b.gen_index() << " ------- " ); + SG_LOG(SG_GENERAL,SG_INFO,"A contour has changed size"); + SG_LOG(SG_GENERAL,SG_INFO,"Bounding box " << rr.getMin() << " , " << rr.getMax() ); + SG_LOG(SG_GENERAL,SG_INFO,"Old size " << subject.ContourSize(0) << " New size " << result.ContourSize(0)); + SG_LOG(SG_GENERAL,SG_INFO,"Old: First node " << subject.GetNode(0,0) << " New: First node " << result.GetNode(0,0)); + SG_LOG(SG_GENERAL,SG_INFO,"Clipping rectangle: " << base.GetContour(0)); + }*/ if ( subject.GetPreserve3D() ) { result.InheritElevations( subject ); result.SetPreserve3D( true ); @@ -33,6 +43,9 @@ tgPolygon tgChopper::Clip( const tgPolygon& subject, result.SetTexMethod( TG_TEX_BY_GEODE, b.get_center_lat() ); } result.SetFlag(type); + if (!subject.IsClosed()) { + result.SetOpen(); + } lock.lock(); bp_map[b.gen_index()].push_back( result ); diff --git a/src/Lib/terragear/tg_contour.cxx b/src/Lib/terragear/tg_contour.cxx index b4aea891..160f92f3 100644 --- a/src/Lib/terragear/tg_contour.cxx +++ b/src/Lib/terragear/tg_contour.cxx @@ -83,8 +83,8 @@ bool tgContour::AreSameSide( const SGGeod& firstpt, const SGGeod& secondpt) cons double y1 = firstpt.getLongitudeDeg(); double y2 = secondpt.getLongitudeDeg(); //Store differences for later - double xdif = x1-x2; - double ydif = y1-y2; + double xdif = x2-x1; + double ydif = y2-y1; /*We describe a line parametrically: x1 (x2-x1) @@ -94,6 +94,16 @@ bool tgContour::AreSameSide( const SGGeod& firstpt, const SGGeod& secondpt) cons with u the parametric coefficient for the second line. Then the line segments intersect if 0 <= t,u <= 1. +To determine t and u we use the approach of Goldman ("Graphics +Gems" as described in Stack Overflow question 563198). + +if r x s = r_x * s_y - r_y * s_x, then + +t = (q - p) x s / (r x s) +and +u = (q - p) x r / (r x s) + +for line 1 = p + t r, line 2 = q + u s */ //Now cycle over all nodes and count how many times we intersect int intersect_ct = 0; @@ -104,11 +114,11 @@ Then the line segments intersect if 0 <= t,u <= 1. double ny1 = node_list[i].getLongitudeDeg(); double nx2 = node_list[i+1].getLatitudeDeg(); double ny2 = node_list[i+1].getLongitudeDeg(); - double nydif = ny1-ny2; - double nxdif = nx1-nx2; + double nydif = ny2-ny1; + double nxdif = nx2-nx1; double denom = xdif*nydif - ydif*nxdif; if (denom != 0) { //Not parallel - double crossx = x1-nx1; double crossy = y1-ny1; + double crossx = nx1-x1; double crossy = ny1-y1; double t = (crossx*nydif - crossy*nxdif)/denom; double u = -1*(xdif*crossy - ydif*crossx)/denom; // We consider that an intersection at the edge of the line have diff --git a/src/Lib/terragear/tg_polygon_clip.cxx b/src/Lib/terragear/tg_polygon_clip.cxx index ec1c8436..c73b7436 100644 --- a/src/Lib/terragear/tg_polygon_clip.cxx +++ b/src/Lib/terragear/tg_polygon_clip.cxx @@ -178,7 +178,11 @@ tgPolygon tgPolygon::Intersect( const tgPolygon& subject, const tgPolygon& clip) result.SetTexParams( subject.GetTexParams() ); result.SetId( subject.GetId() ); result.SetPreserve3D( subject.GetPreserve3D() ); - + if(subject.IsClosed()) { + result.SetClosed(); + } else { + result.SetOpen(); + } return result; } From 62359aede25420ee456a55f67b5506952f9873c2 Mon Sep 17 00:00:00 2001 From: "James.Hester" Date: Sun, 30 Dec 2018 10:05:35 +1100 Subject: [PATCH 6/8] Adjust limits for line intersection calculation to catch points slightly past the end of a line segment. This is to catch situations where the cliff line has been chopped at a bucket edge and the two candidate points are on the edge. --- src/Lib/terragear/tg_contour.cxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Lib/terragear/tg_contour.cxx b/src/Lib/terragear/tg_contour.cxx index 160f92f3..9cd4e7de 100644 --- a/src/Lib/terragear/tg_contour.cxx +++ b/src/Lib/terragear/tg_contour.cxx @@ -121,11 +121,11 @@ for line 1 = p + t r, line 2 = q + u s double crossx = nx1-x1; double crossy = ny1-y1; double t = (crossx*nydif - crossy*nxdif)/denom; double u = -1*(xdif*crossy - ydif*crossx)/denom; - // We consider that an intersection at the edge of the line have - // not crossed - // over, that is, they lie on the same side, so we do not - // include equality in the comparisons - if (t > 0.0 && t < 1.0 && u > 0.0 && u < 1.0) intersect_ct++; + // We consider that an intersection at the edge of the line has + // crossed + // over, that is, they lie on opposite sides. This way we capture + // places where the chopper has clipped a cliff on the tile edge + if (t > -0.0001 && t < 1.0001 && u > -0.0001 && u < 1.0001) intersect_ct++; } } } From 990417c00094dc5aa74a26d87820288aed8b2a1a Mon Sep 17 00:00:00 2001 From: "James.Hester" Date: Sun, 30 Dec 2018 10:07:26 +1100 Subject: [PATCH 7/8] Teach terrafit about rectified height files. --- src/Prep/TerraFit/CMakeLists.txt | 3 +++ src/Prep/TerraFit/terrafit.cc | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/Prep/TerraFit/CMakeLists.txt b/src/Prep/TerraFit/CMakeLists.txt index 1337d3b0..bcd7e752 100644 --- a/src/Prep/TerraFit/CMakeLists.txt +++ b/src/Prep/TerraFit/CMakeLists.txt @@ -1,8 +1,11 @@ add_executable(terrafit terrafit.cc) +include_directories(${PROJECT_SOURCE_DIR}/src/Lib/terragear) + target_link_libraries(terrafit Array Terra + terragear ${ZLIB_LIBRARY} ${SIMGEAR_CORE_LIBRARIES} ${SIMGEAR_CORE_LIBRARY_DEPENDENCIES}) diff --git a/src/Prep/TerraFit/terrafit.cc b/src/Prep/TerraFit/terrafit.cc index 4743133c..6ecce7d5 100644 --- a/src/Prep/TerraFit/terrafit.cc +++ b/src/Prep/TerraFit/terrafit.cc @@ -239,7 +239,7 @@ void walk_path(const SGPath& path) { return; } - if ((path.lower_extension() == "arr") || (path.complete_lower_extension() == "arr.gz")) { + if ((path.lower_extension() == "arr") || (path.complete_lower_extension() == "arr.rectified.gz")) { SG_LOG(SG_GENERAL, SG_DEBUG, "will queue " << path); queue_fit_file(path); } else if (path.isDir()) { From 09d61732cd6dd58527759e851f73db3a0f388775 Mon Sep 17 00:00:00 2001 From: "James.Hester" Date: Sun, 30 Dec 2018 10:08:59 +1100 Subject: [PATCH 8/8] Added a small tool to rectify heights using cliff contours placed in the height array directory. --- src/Lib/Array/rectify_height.cxx | 131 +++++++++++++++++++++++++++++++ 1 file changed, 131 insertions(+) create mode 100644 src/Lib/Array/rectify_height.cxx diff --git a/src/Lib/Array/rectify_height.cxx b/src/Lib/Array/rectify_height.cxx new file mode 100644 index 00000000..06b35323 --- /dev/null +++ b/src/Lib/Array/rectify_height.cxx @@ -0,0 +1,131 @@ +//rectify_height.cxx +#ifdef HAVE_CONFIG_H +# include +#endif + +#include +#include +#include +#include + +#include +#include +#include +#include +#include "array.hxx" + +/* This program will find all height files in the directory provided, +and if they have an associated cliff file, will adjust and then output +the heights. Single threaded to avoid reading/writing the same file. */ + +// display usage and exit +static void usage( const std::string name ) { + SG_LOG(SG_GENERAL, SG_ALERT, "Usage: " << name); + SG_LOG(SG_GENERAL, SG_ALERT, " --work-dir="); + SG_LOG(SG_GENERAL, SG_ALERT, " --height-dir="); + SG_LOG(SG_GENERAL, SG_ALERT, "[ --tile-id=]"); + SG_LOG(SG_GENERAL, SG_ALERT, " --min-lon="); + SG_LOG(SG_GENERAL, SG_ALERT, " --max-lon="); + SG_LOG(SG_GENERAL, SG_ALERT, " --min-lat="); + SG_LOG(SG_GENERAL, SG_ALERT, " --max-lat="); + SG_LOG(SG_GENERAL, SG_ALERT, "[ --min-dist=]"); + exit(-1); +} + +int main( int argc, char **argv) { + std::string output_dir = "."; + std::string work_dir = "."; + std::string height_dir = ""; + SGGeod min,max; + long tile_id = -1; + double bad_zone = 30; //distance in m from cliff to rectify + + sglog().setLogLevels(SG_ALL,SG_INFO); + sglog().set_log_priority(SG_DEBUG); + + // + // Parse the command-line arguments. + // + int arg_pos; + for (arg_pos = 1; arg_pos < argc; arg_pos++) { + std::string arg = argv[arg_pos]; + + if (arg.find("--work-dir=") == 0) { + work_dir = arg.substr(11); + } else if (arg.find("--height-dir=") == 0) { + height_dir = arg.substr(13); + } else if (arg.find("--tile-id=") == 0) { + tile_id = atol(arg.substr(10).c_str()); + } else if ( arg.find("--min-lon=") == 0 ) { + min.setLongitudeDeg(atof( arg.substr(10).c_str() )); + } else if ( arg.find("--max-lon=") == 0 ) { + max.setLongitudeDeg(atof( arg.substr(10).c_str() )); + } else if ( arg.find("--min-lat=") == 0 ) { + min.setLatitudeDeg(atof( arg.substr(10).c_str() )); + } else if ( arg.find("--max-lat=") == 0 ) { + max.setLatitudeDeg(atof( arg.substr(10).c_str() )); + } else if ( arg.find("--min-dist=") == 0) { + bad_zone = atof(arg.substr(11).c_str()); + } else if (arg.find("--") == 0) { + usage(argv[0]); + } else { + break; + } + } + + SG_LOG(SG_GENERAL, SG_ALERT, "Working directory is " << work_dir); + SG_LOG(SG_GENERAL, SG_ALERT, "Heights are in " << height_dir); + SG_LOG(SG_GENERAL, SG_ALERT, "Rectification zone within " << bad_zone << "m of cliffs"); + + if ( tile_id > 0 ) { + SG_LOG(SG_GENERAL, SG_ALERT, "Tile id is " << tile_id); + } else { + if (min.isValid() && max.isValid() && (min != max)) + { + SG_LOG(SG_GENERAL, SG_ALERT, "Longitude = " << min.getLongitudeDeg() << ':' << max.getLongitudeDeg()); + SG_LOG(SG_GENERAL, SG_ALERT, "Latitude = " << min.getLatitudeDeg() << ':' << max.getLatitudeDeg()); + } else + { + SG_LOG(SG_GENERAL, SG_ALERT, "Lon/Lat unset or wrong"); + exit(1); + } + } + + // Now generate the buckets + std::vector bucketList; + if (tile_id == -1) { + // build all the tiles in an area + SG_LOG(SG_GENERAL, SG_ALERT, "Fixing heights within given bounding box"); + + SGBucket b_min( min ); + SGBucket b_max( max ); + + if ( b_min == b_max ) { + bucketList.push_back( b_min ); + } else { + SG_LOG(SG_GENERAL, SG_ALERT, " adjustment area spans tile boundaries"); + sgGetBuckets( min, max, bucketList ); + } + } else { + // adjust the specified tile + SG_LOG(SG_GENERAL, SG_ALERT, "Adjusting tile " << tile_id); + bucketList.push_back( SGBucket( tile_id ) ); + } + + // Finally, loop over the buckets adjusting heights + for (unsigned int i=0; i < bucketList.size(); i++) { + SGBucket bucket = bucketList[i]; + TGArray array; + std::string base_path = bucket.gen_base_path(); + std::string array_path = work_dir + "/" + height_dir + "/" + base_path + "/" + bucket.gen_index_str(); + if (!array.open(array_path)) { + SG_LOG(SG_GENERAL,SG_DEBUG, "Failed to open array file " << array_path); + continue; + } + array.parse(bucket); + array.close(); + array.remove_voids(); + array.rectify_heights(bad_zone); + array.write_bin(work_dir + "/" + height_dir,true,bucket); + } +}