From 005fd6886a83f6449a2eac681194a521bd84751f Mon Sep 17 00:00:00 2001 From: "James.Hester" Date: Sat, 6 Oct 2018 18:46:04 +1000 Subject: [PATCH] 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)