From 8509dc1c33f29f69fde08aa46f1b98facf2686cb Mon Sep 17 00:00:00 2001
From: "James.Hester" <jxh@ansto.gov.au>
Date: Sat, 15 Sep 2018 19:39:02 +1000
Subject: [PATCH 01/15] 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 <config.h>
+#endif
+
+#include <simgear/compiler.h>
+
+#include <stdlib.h> // atof()
+#include <iostream>
+
+#ifdef SG_HAVE_STD_INCLUDES
+# include <cerrno>
+#else
+# include <errno.h>
+#endif
+
+#ifdef _MSC_VER
+# include <direct.h>
+#endif
+
+#include <simgear/constants.h>
+#include <simgear/io/lowlevel.hxx>
+#include <simgear/misc/sg_dir.hxx>
+#include <simgear/debug/logstream.hxx>
+
+
+#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 <config.h>
+#endif
+
+#include <simgear/compiler.h>
+
+#include "srtmbase.hxx"
+
+#include <zlib.h>
+
+#include <string>
+
+#include <simgear/bucket/newbucket.hxx>
+#include <simgear/misc/sg_path.hxx>
+
+#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 <config.h>
+#endif
+
+#include <simgear/compiler.h>
+
+#include <cstdlib>
+#include <string>
+#include <iostream>
+
+#include <simgear/bucket/newbucket.hxx>
+#include <simgear/debug/logstream.hxx>
+#include <simgear/misc/sg_path.hxx>
+#include <simgear/misc/sg_dir.hxx>
+
+#include <Include/version.h>
+#include <HGT/dted.hxx>
+
+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] << " <resolution> <dted_file> <work_dir>" << 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" <jxh@ansto.gov.au>
Date: Sat, 6 Oct 2018 18:46:04 +1000
Subject: [PATCH 02/15] 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<count; i++ ) {
+ poly.LoadFromGzFile( fp );
+ if ( poly.Contours()==1 ) { //should always have one contour
+ cliffs_list.push_back(poly.GetContour(0));
+ } else {
+ SG_LOG( SG_GENERAL, SG_WARN, " Found " << poly.Contours() << " contours in " << p.realpath() );
+ }
+ }
+ }
+ }
+}
+
+
void
TGArray::unload( void ) {
if (array_in) {
@@ -316,14 +359,15 @@ void TGArray::remove_voids( ) {
// Return the elevation of the closest non-void grid point to lon, lat
+// Lon, lat in arcsec
double TGArray::closest_nonvoid_elev( double lon, double lat ) const {
double mindist = 99999999999.9;
double minelev = -9999.0;
- SGGeod p0 = SGGeod::fromDeg( lon, lat );
+ SGGeod p0 = SGGeod::fromDeg( lon/3600.0, lat/3600.0 );
for ( int row = 0; row < rows; row++ ) {
for ( int col = 0; col < cols; col++ ) {
- SGGeod p1 = SGGeod::fromDeg( originx + col * col_step, originy + row * row_step );
+ SGGeod p1 = SGGeod::fromDeg( (originx + col * col_step)/3600.0, (originy + row * row_step)/3600.0 );
double dist = SGGeodesy::distanceM( p0, p1 );
double elev = get_array_elev(col, row);
if ( dist < mindist && elev > -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)<SG_EPSILON && fabs(lat1-lat2)<SG_EPSILON) return true;
+ SGGeod pt1 = SGGeod::fromDeg(lon1,lat1);
+ SGGeod pt2 = SGGeod::fromDeg(lon2,lat2);
+ bool same_side = true;
+ for (int i=0;i<cliffs_list.size();i++) {
+ bool check_result = cliffs_list[i].AreSameSide(pt1,pt2);
+ if(!check_result) {
+ SG_LOG(SG_GENERAL, SG_DEBUG, "Cliff " << i <<":" <<pt1 << " and " << pt2 << " on opposite sides");
+ same_side = false;
+ break;
+ }
+ }
+ return same_side;
+}
TGArray::~TGArray( void )
{
@@ -486,3 +694,4 @@ bool TGArray::is_open() const
return false;
}
}
+
diff --git a/src/Lib/Array/array.hxx b/src/Lib/Array/array.hxx
index 654d8ad3..16d623a2 100644
--- a/src/Lib/Array/array.hxx
+++ b/src/Lib/Array/array.hxx
@@ -27,6 +27,11 @@
#include <simgear/bucket/newbucket.hxx>
#include <simgear/math/sg_types.hxx>
#include <simgear/io/iostreams/sgstream.hxx>
+#include <simgear/misc/sg_dir.hxx>
+#include <boost/foreach.hpp>
+
+#include "tg_contour.hxx"
+#include "tg_polygon.hxx"
class TGArray {
@@ -52,6 +57,8 @@ private:
std::vector<SGGeod> corner_list;
std::vector<SGGeod> 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 <simgear/compiler.h>
#include <simgear/bucket/newbucket.hxx>
#include <simgear/misc/sg_path.hxx>
+#include <simgear/debug/logstream.hxx>
#include <sys/types.h>
#include <iostream>
@@ -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<node_list.size()-1;i++) {
+ double nx1 = node_list[i].getLatitudeDeg();
+ 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 denom = xdif*nydif - ydif*nxdif;
+ if (denom != 0) { //Not parallel
+ 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++;
+ }
+ }
+ }
+ 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> 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" <jxh@ansto.gov.au>
Date: Sun, 21 Oct 2018 10:07:23 +1100
Subject: [PATCH 03/15] 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 <cstring>
+#include <iomanip> //for setprecision
#include <simgear/compiler.h>
#include <simgear/io/iostreams/sgstream.hxx>
@@ -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 << "),("<<lon2<<","<<lat2<<")");
default: // all corners
dx = xlocal - xindex;
dy = ylocal - yindex;
diff --git a/src/Lib/terragear/tg_chopper.cxx b/src/Lib/terragear/tg_chopper.cxx
index 8ed5e80a..f9a52e9f 100644
--- a/src/Lib/terragear/tg_chopper.cxx
+++ b/src/Lib/terragear/tg_chopper.cxx
@@ -11,8 +11,7 @@
#include "tg_misc.hxx"
tgPolygon tgChopper::Clip( const tgPolygon& subject,
- const std::string& type,
- SGBucket& b )
+ const std::string& type, SGBucket& b)
{
tgPolygon base, result;
@@ -22,7 +21,7 @@ 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);
if ( result.Contours() > 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<SGG
p0 = subject.GetNode( subject.GetSize() - 1 );
p1 = subject.GetNode( 0 );
- // add start of segment
- result.AddNode( p0 );
-
- // add intermediate points
- AddIntermediateNodes( p0, p1, nodes, result, SG_EPSILON*10, SG_EPSILON*4 );
-
- // maintain original hole flag setting
+ if (!subject.GetOpen()) {
+ // add start of segment
+ result.AddNode( p0 );
+
+ // add intermediate points
+ AddIntermediateNodes( p0, p1, nodes, result, SG_EPSILON*10, SG_EPSILON*4 );
+ }
+ // maintain original hole and openness flag setting
result.SetHole( subject.GetHole() );
+ result.SetOpen( subject.GetOpen() );
return result;
}
@@ -878,15 +884,17 @@ tgContour tgContour::AddColinearNodes( const tgContour& subject, bool preserve3d
p0 = subject.GetNode( subject.GetSize() - 1 );
p1 = subject.GetNode( 0 );
-
- // add start of segment
- result.AddNode( p0 );
-
- // add intermediate points
- AddIntermediateNodes( p0, p1, preserve3d, nodes, result, SG_EPSILON*10, SG_EPSILON*4 );
-
- // maintain original hole flag setting
+
+ if(!subject.GetOpen()) {
+ // add start of segment
+ result.AddNode( p0 );
+
+ // add intermediate points
+ AddIntermediateNodes( p0, p1, preserve3d, nodes, result, SG_EPSILON*10, SG_EPSILON*4 );
+ }
+ // maintain original hole and open flag settings
result.SetHole( subject.GetHole() );
+ result.SetOpen( subject.GetOpen() );
return result;
}
diff --git a/src/Lib/terragear/tg_contour.hxx b/src/Lib/terragear/tg_contour.hxx
index 37e2767e..cb356aaa 100644
--- a/src/Lib/terragear/tg_contour.hxx
+++ b/src/Lib/terragear/tg_contour.hxx
@@ -29,6 +29,7 @@ class tgContour
public:
tgContour() {
hole = false;
+ keep_open = false;
}
void Erase() {
@@ -38,10 +39,19 @@ public:
void SetHole( bool h ) {
hole = h;
}
+
bool GetHole( void ) const {
return hole;
}
+ bool GetOpen(void) const {
+ return keep_open;
+ }
+
+ void SetOpen(bool o) {
+ keep_open = o;
+ }
+
unsigned int GetSize( void ) const {
return node_list.size();
}
@@ -138,6 +148,7 @@ public:
private:
std::vector<SGGeod> node_list;
bool hole;
+ bool keep_open; //If non-closed contour, keep open
};
typedef std::vector <tgContour> 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" <jxh@ansto.gov.au>
Date: Sun, 28 Oct 2018 23:47:23 +1100
Subject: [PATCH 04/15] 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<int> TGArray::collect_bad_points() {
+ //Find and remember all points that are bad because they are
+ //too close to a cliff
+ std::vector<int> bad_points; //local to avoid multi-thread issues
+ for(int horiz=0;horiz<cols;horiz++) {
+ double lon = (originx + col_step*horiz)/3600;
+ for(int vert=0;vert<rows;vert++) {
+ double lat = (originy + row_step*vert)/3600;
+ if(is_near_cliff(lon,lat)) {
+ bad_points.push_back(horiz+vert*cols);
+ }
+ }
+ }
+ return bad_points;
+}
+// Check to see if the specified grid point is bad
+bool TGArray::is_bad_point(const int xgrid, const int ygrid, const std::vector<int> 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<int> 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<int> 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<cliffs_list.size();i++) {
+ double dist = cliffs_list[i].MinDist(pt1);
+ if (dist < mindist) return true;
+ }
+ return false;
+}
+
TGArray::~TGArray( void )
{
if (in_data) {
diff --git a/src/Lib/Array/array.hxx b/src/Lib/Array/array.hxx
index 16d623a2..cfa2b1bb 100644
--- a/src/Lib/Array/array.hxx
+++ b/src/Lib/Array/array.hxx
@@ -60,6 +60,13 @@ private:
// list of cliff contours
tgcontour_list cliffs_list;
void parse_bin();
+
+ // Routines for height rectification
+ void rectify_heights();
+ std::vector<int> collect_bad_points();
+ bool is_bad_point(const int xgrid, const int ygrid, const std::vector<int> bad_points) const;
+ double rectify_point(const int xgrid, const int ygrid, const std::vector<int> 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 <simgear/math/sg_geodesy.hxx>
+#include <simgear/math/SGGeometry.hxx>
#include <simgear/io/lowlevel.hxx>
#include <simgear/debug/logstream.hxx>
@@ -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<node_list.size(); i++) {
a = SGGeod_ToSGVec2d( node_list[i] );
b = SGGeod_ToSGVec2d( node_list[j] );
-
area += (b.x() + a.x()) * (b.y() - a.y());
j=i;
}
@@ -122,7 +122,27 @@ Then the line segments intersect if 0 <= t,u <= 1.
bool isinter = (intersect_ct%2 == 0);
return isinter;
}
-
+
+double tgContour::MinDist(const SGGeod& probe) const
+{
+ SGVec3d probexyz;
+ SGGeodesy::SGGeodToCart(probe,probexyz);
+ double mindist = 100000.0;
+ double dist;
+ if (node_list.size()) {
+ int j = node_list.size() - 1;
+ for (int i=0;i<j;i++) {
+ SGVec3d start,end;
+ SGGeodesy::SGGeodToCart(node_list[i],start);
+ SGGeodesy::SGGeodToCart(node_list[i+1],end);
+ SGLineSegment<double> piece = SGLineSegment<double>(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 <simgear/compiler.h>
#include <simgear/math/sg_types.hxx>
+#include <simgear/math/SGLineSegment.hxx>
+#include <simgear/math/SGGeodesy.hxx>
#include <boost/concept_check.hpp>
#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" <jxh@ansto.gov.au>
Date: Mon, 17 Dec 2018 23:50:28 +1100
Subject: [PATCH 05/15] 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<int> TGArray::collect_bad_points() {
+std::vector<int> 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<int> bad_points; //local to avoid multi-thread issues
@@ -392,7 +433,7 @@ std::vector<int> TGArray::collect_bad_points() {
double lon = (originx + col_step*horiz)/3600;
for(int vert=0;vert<rows;vert++) {
double lat = (originy + row_step*vert)/3600;
- if(is_near_cliff(lon,lat)) {
+ if(is_near_cliff(lon,lat,bad_zone)) {
bad_points.push_back(horiz+vert*cols);
}
}
@@ -412,10 +453,10 @@ bool TGArray::is_bad_point(const int xgrid, const int ygrid, const std::vector<i
//This may collide with other threads, but as they will both be writing
//the correct height, this is harmless.
-void TGArray::rectify_heights() {
+void TGArray::rectify_heights(const double bad_zone) {
double new_ht;
std::vector<int> 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 <pt_cnt; pti++) {
+ totht = totht + get_array_elev(xgrid+corners[pti][0],ygrid+corners[pti][1]);
+ }
+ height = totht/pt_cnt;
+ } else {
+
// We have three points, calculate the height
+ // Set anything very negative to zero
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);
+ if (corner < -9000) corner = 0;
+ if (horiz < -9000) horiz = 0;
+ if (vert < -9000) vert = 0;
+ height = horiz + (vert - corner);
+ }
std::cout << xgrid << "," << ygrid << ": was " << original_height << " , now " << height << std::endl;
return height;
}
@@ -768,13 +834,12 @@ bool TGArray::check_points (const double lon1, const double lat1, const double l
//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 {
+bool TGArray::is_near_cliff(const double lon1, const double lat1, const double bad_zone) const {
if (cliffs_list.size()==0) return false;
SGGeod pt1 = SGGeod::fromDeg(lon1,lat1);
- double mindist = 30; // 1 arcsec
for (int i=0;i<cliffs_list.size();i++) {
double dist = cliffs_list[i].MinDist(pt1);
- if (dist < mindist) return true;
+ if (dist < bad_zone) return true;
}
return false;
}
diff --git a/src/Lib/Array/array.hxx b/src/Lib/Array/array.hxx
index cfa2b1bb..9de84269 100644
--- a/src/Lib/Array/array.hxx
+++ b/src/Lib/Array/array.hxx
@@ -47,6 +47,8 @@ private:
// number of columns and rows
int cols, rows;
+ // Whether or not the input data have been rectified
+ bool rectified;
// Distance between column and row data points (in arc seconds)
double col_step, row_step;
@@ -62,11 +64,10 @@ private:
void parse_bin();
// Routines for height rectification
- void rectify_heights();
- std::vector<int> collect_bad_points();
+ std::vector<int> collect_bad_points(const double bad_zone);
bool is_bad_point(const int xgrid, const int ygrid, const std::vector<int> bad_points) const;
double rectify_point(const int xgrid, const int ygrid, const std::vector<int> 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" <jxh@ansto.gov.au>
Date: Sun, 30 Dec 2018 10:05:35 +1100
Subject: [PATCH 06/15] 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" <jxh@ansto.gov.au>
Date: Sun, 30 Dec 2018 10:07:26 +1100
Subject: [PATCH 07/15] 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" <jxh@ansto.gov.au>
Date: Sun, 30 Dec 2018 10:08:59 +1100
Subject: [PATCH 08/15] 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 <config.h>
+#endif
+
+#include <simgear/compiler.h>
+#include <simgear/bucket/newbucket.hxx>
+#include <simgear/misc/sg_path.hxx>
+#include <simgear/debug/logstream.hxx>
+
+#include <sys/types.h>
+#include <iostream>
+#include <string.h>
+#include <stdlib.h>
+#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=<directory>");
+ SG_LOG(SG_GENERAL, SG_ALERT, " --height-dir=<directory>");
+ SG_LOG(SG_GENERAL, SG_ALERT, "[ --tile-id=<id>]");
+ SG_LOG(SG_GENERAL, SG_ALERT, " --min-lon=<degrees>");
+ SG_LOG(SG_GENERAL, SG_ALERT, " --max-lon=<degrees>");
+ SG_LOG(SG_GENERAL, SG_ALERT, " --min-lat=<degrees>");
+ SG_LOG(SG_GENERAL, SG_ALERT, " --max-lat=<degrees>");
+ SG_LOG(SG_GENERAL, SG_ALERT, "[ --min-dist=<float>]");
+ 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<SGBucket> 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);
+ }
+}
From e2a6d06de13caab4295b25c0bcf98bcc040e27aa Mon Sep 17 00:00:00 2001
From: "James.Hester" <jxh@ansto.gov.au>
Date: Sun, 30 Dec 2018 13:46:42 +1100
Subject: [PATCH 09/15] Added cliff-decode source and make instructions.
---
src/Prep/Cliff/CMakeLists.txt | 14 ++
src/Prep/Cliff/cliff-decode.cxx | 434 ++++++++++++++++++++++++++++++++
2 files changed, 448 insertions(+)
create mode 100644 src/Prep/Cliff/CMakeLists.txt
create mode 100644 src/Prep/Cliff/cliff-decode.cxx
diff --git a/src/Prep/Cliff/CMakeLists.txt b/src/Prep/Cliff/CMakeLists.txt
new file mode 100644
index 00000000..cda27455
--- /dev/null
+++ b/src/Prep/Cliff/CMakeLists.txt
@@ -0,0 +1,14 @@
+include_directories(${GDAL_INCLUDE_DIR})
+add_executable(cliff-decode cliff-decode.cxx)
+
+target_link_libraries(cliff-decode
+ ${GDAL_LIBRARY}
+ terragear
+ ${ZLIB_LIBRARY}
+ ${CMAKE_THREAD_LIBS_INIT}
+ ${SIMGEAR_CORE_LIBRARIES}
+ ${SIMGEAR_CORE_LIBRARY_DEPENDENCIES}
+ ${Boost_LIBRARIES}
+)
+
+install(TARGETS cliff-decode RUNTIME DESTINATION bin)
diff --git a/src/Prep/Cliff/cliff-decode.cxx b/src/Prep/Cliff/cliff-decode.cxx
new file mode 100644
index 00000000..4b38298c
--- /dev/null
+++ b/src/Prep/Cliff/cliff-decode.cxx
@@ -0,0 +1,434 @@
+// cliff-decode.cxx -- process OGR layers and extract lines
+// representing discontinuities, clipping
+// against and sorting
+// them into the relevant tiles.
+//
+// Written by James Hester starting 2018
+//
+// Based on ogr-decode by Ralf Gerlich.
+//
+// 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.
+//
+#include <string>
+#include <map>
+
+#include <boost/thread.hpp>
+#include <ogrsf_frmts.h>
+#include <gdal_priv.h>
+
+#include <simgear/compiler.h>
+#include <simgear/threads/SGThread.hxx>
+#include <simgear/threads/SGQueue.hxx>
+#include <simgear/debug/logstream.hxx>
+#include <simgear/math/sg_geodesy.hxx>
+#include <simgear/misc/sg_path.hxx>
+
+#include <Include/version.h>
+
+#include <terragear/tg_chopper.hxx>
+#include <terragear/tg_shapefile.hxx>
+
+using std::string;
+
+int continue_on_errors=0;
+int num_threads=1;
+bool use_attribute_query=false;
+string attribute_query;
+int start_record=0;
+bool use_spatial_query=false;
+int seperate_segments = 0;
+double spat_min_x, spat_min_y, spat_max_x, spat_max_y;
+bool save_shapefiles=false;
+SGLockedQueue<OGRFeature *> global_workQueue;
+
+class Decoder : public SGThread
+{
+public:
+ Decoder( OGRCoordinateTransformation *poct, tgChopper& c ) : chopper(c) {
+ poCT = poct;
+ }
+
+private:
+ virtual void run();
+
+ void processLineString(OGRLineString* poGeometry);
+
+private:
+ // The transformation for each geometry object
+ OGRCoordinateTransformation *poCT;
+
+ // Store the reults per tile
+ tgChopper& chopper;
+
+};
+
+void Decoder::processLineString(OGRLineString* poGeometry)
+{
+ tgContour line;
+
+ SGGeod p0, p1;
+ double heading, dist, az2;
+ int i, j, numPoints, numSegs;
+ double max_dist;
+
+ numPoints = poGeometry->getNumPoints();
+ if (numPoints < 2) {
+ SG_LOG( SG_GENERAL, SG_WARN, "Skipping line with less than two points" );
+ return;
+ }
+
+ heading = SGGeodesy::courseDeg( p1, p0 );
+
+ // now add the middle points : if they are too far apart, add intermediate nodes
+ for ( i=0;i<numPoints;i++) {
+ p0 = SGGeod::fromDeg( poGeometry->getX(i), poGeometry->getY(i) );
+ line.AddNode( p0 );
+ }
+
+ // Do not close this contour
+ line.SetOpen(true);
+
+ // clip the contour. Do not use usual clipper, as we wish to preserve the
+ // whole contour, even if it goes outside the bounding box, as height
+ // calculations along the boundary might be affected by it.
+ tgPolygon open_poly;
+ open_poly.SetOpen(); //do not try to close this one up
+ open_poly.AddContour(line);
+ chopper.Add( open_poly, "Cliffs" );
+}
+
+void Decoder::run()
+{
+ // as long as we have geometry to parse, do so
+ while (!global_workQueue.empty()) {
+ OGRFeature *poFeature = global_workQueue.pop();
+ SG_LOG( SG_GENERAL, SG_BULK, " remaining features is " << global_workQueue.size() );
+
+ if ( poFeature ) {
+ OGRGeometry *poGeometry = poFeature->GetGeometryRef();
+
+ if (poGeometry==NULL) {
+ SG_LOG( SG_GENERAL, SG_INFO, "Found feature without geometry!" );
+ if (!continue_on_errors) {
+ SG_LOG( SG_GENERAL, SG_ALERT, "Aborting!" );
+ exit( 1 );
+ } else {
+ continue;
+ }
+ }
+
+ OGRwkbGeometryType geoType=wkbFlatten(poGeometry->getGeometryType());
+ if (geoType!=wkbLineString && geoType!=wkbMultiLineString) {
+ SG_LOG( SG_GENERAL, SG_INFO, "Non-line feature" );
+ continue;
+ }
+
+ poGeometry->transform( poCT );
+
+ switch (geoType) {
+ case wkbLineString: {
+ SG_LOG( SG_GENERAL, SG_DEBUG, "LineString feature" );
+ processLineString((OGRLineString*)poGeometry);
+ break;
+ }
+ case wkbMultiLineString: {
+ SG_LOG( SG_GENERAL, SG_DEBUG, "MultiLineString feature" );
+ OGRMultiLineString* multilines=(OGRMultiLineString*)poGeometry;
+ for (int i=0;i<multilines->getNumGeometries();i++) {
+ processLineString((OGRLineString*)(multilines->getGeometryRef(i)));
+ }
+ break;
+ }
+ default:
+ /* Ignore unhandled objects */
+ break;
+ }
+
+ OGRFeature::DestroyFeature( poFeature );
+ }
+ }
+}
+
+// Main Thread
+void processLayer(OGRLayer* poLayer, tgChopper& results )
+{
+ int feature_count=poLayer->GetFeatureCount();
+
+ if (feature_count!=-1 && start_record>0 && start_record>=feature_count) {
+ SG_LOG( SG_GENERAL, SG_ALERT, "Layer has only " << feature_count << " records, but start record is set to " << start_record );
+ exit( 1 );
+ }
+
+ /* determine the indices of the required columns */
+ OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
+ string layername=poFDefn->GetName();
+
+ /* setup a transformation to WGS84 */
+ OGRSpatialReference *oSourceSRS, oTargetSRS;
+ oSourceSRS=poLayer->GetSpatialRef();
+ if (oSourceSRS == NULL) {
+ SG_LOG( SG_GENERAL, SG_ALERT, "Layer " << layername << " has no defined spatial reference system" );
+ exit( 1 );
+ }
+
+ char* srsWkt;
+ oSourceSRS->exportToWkt(&srsWkt);
+ SG_LOG( SG_GENERAL, SG_DEBUG, "Source spatial reference system: " << srsWkt );
+ OGRFree(srsWkt);
+
+ oTargetSRS.SetWellKnownGeogCS( "WGS84" );
+
+ OGRCoordinateTransformation *poCT = OGRCreateCoordinateTransformation(oSourceSRS, &oTargetSRS);
+
+ /* setup attribute and spatial queries */
+ if (use_spatial_query) {
+ double trans_min_x,trans_min_y,trans_max_x,trans_max_y;
+ /* do a simple reprojection of the source SRS */
+ OGRCoordinateTransformation *poCTinverse;
+
+ poCTinverse = OGRCreateCoordinateTransformation(&oTargetSRS, oSourceSRS);
+
+ trans_min_x=spat_min_x;
+ trans_min_y=spat_min_y;
+ trans_max_x=spat_max_x;
+ trans_max_y=spat_max_y;
+
+ poCTinverse->Transform(1,&trans_min_x,&trans_min_y);
+ poCTinverse->Transform(1,&trans_max_x,&trans_max_y);
+
+ poLayer->SetSpatialFilterRect(trans_min_x, trans_min_y,
+ trans_max_x, trans_max_y);
+ }
+
+ if (use_attribute_query) {
+ if (poLayer->SetAttributeFilter(attribute_query.c_str()) != OGRERR_NONE) {
+ SG_LOG( SG_GENERAL, SG_ALERT, "Error in query expression '" << attribute_query << "'" );
+ exit( 1 );
+ }
+ }
+
+ // Generate the work queue for this layer
+ OGRFeature *poFeature;
+ poLayer->SetNextByIndex(start_record);
+ while ( ( poFeature = poLayer->GetNextFeature()) != NULL )
+ {
+ global_workQueue.push( poFeature );
+ }
+
+ // Now process the workqueue with threads
+ std::vector<Decoder *> decoders;
+ for (int i=0; i<num_threads; i++) {
+ Decoder* decoder = new Decoder( poCT, results );
+ decoder->start();
+ decoders.push_back( decoder );
+ }
+
+ // Then wait until they are finished
+ for (unsigned int i=0; i<decoders.size(); i++) {
+ decoders[i]->join();
+ }
+
+ OCTDestroyCoordinateTransformation ( poCT );
+}
+
+void usage(char* progname) {
+ SG_LOG( SG_GENERAL, SG_ALERT, "Usage: " <<
+ progname << " [options...] <work_dir> <datasource> [<layername>...]" );
+ SG_LOG( SG_GENERAL, SG_ALERT, "Options:" );
+ SG_LOG( SG_GENERAL, SG_ALERT, "--log-level priority" );
+ SG_LOG( SG_GENERAL, SG_ALERT, " Width in priority being bulk|debug|info|warn|alert" );
+ SG_LOG( SG_GENERAL, SG_ALERT, "--continue-on-errors" );
+ SG_LOG( SG_GENERAL, SG_ALERT, " Continue even if the file seems fishy" );
+ SG_LOG( SG_GENERAL, SG_ALERT, "--max-segment max_segment_length" );
+ SG_LOG( SG_GENERAL, SG_ALERT, " Maximum segment length in meters" );
+ SG_LOG( SG_GENERAL, SG_ALERT, "--start-record record_no" );
+ SG_LOG( SG_GENERAL, SG_ALERT, " Start processing at the specified record number (first record num=0)" );
+ SG_LOG( SG_GENERAL, SG_ALERT, "--where attrib_query" );
+ SG_LOG( SG_GENERAL, SG_ALERT, " Use an attribute query (like SQL WHERE)" );
+ SG_LOG( SG_GENERAL, SG_ALERT, "--spat xmin ymin xmax ymax" );
+ SG_LOG( SG_GENERAL, SG_ALERT, " spatial query extents" );
+ SG_LOG( SG_GENERAL, SG_ALERT, "--threads" );
+ SG_LOG( SG_GENERAL, SG_ALERT, " Enable multithreading with user specified number of threads" );
+ SG_LOG( SG_GENERAL, SG_ALERT, "--all-threads" );
+ SG_LOG( SG_GENERAL, SG_ALERT, " Enable multithreading with all available cpu cores" );
+ SG_LOG( SG_GENERAL, SG_ALERT, "" );
+ SG_LOG( SG_GENERAL, SG_ALERT, "<work_dir>" );
+ SG_LOG( SG_GENERAL, SG_ALERT, " Directory to put the cliff files in" );
+ SG_LOG( SG_GENERAL, SG_ALERT, "<datasource>" );
+ SG_LOG( SG_GENERAL, SG_ALERT, " The datasource from which to fetch the data" );
+ SG_LOG( SG_GENERAL, SG_ALERT, "<layername>..." );
+ SG_LOG( SG_GENERAL, SG_ALERT, " The layers to process." );
+ SG_LOG( SG_GENERAL, SG_ALERT, " If no layer is given, all layers in the datasource are used" );
+ exit(-1);
+}
+
+void
+setLoggingPriority (const char * p)
+{
+ if (p == 0)
+ return;
+ string priority = p;
+ if (priority == "bulk") {
+ sglog().set_log_priority(SG_BULK);
+ } else if (priority == "debug") {
+ sglog().set_log_priority(SG_DEBUG);
+ } else if (priority.empty() || priority == "info") { // default
+ sglog().set_log_priority(SG_INFO);
+ } else if (priority == "warn") {
+ sglog().set_log_priority(SG_WARN);
+ } else if (priority == "alert") {
+ sglog().set_log_priority(SG_ALERT);
+ } else {
+ SG_LOG(SG_GENERAL, SG_WARN, "Unknown logging priority " << priority);
+ }
+}
+
+int main( int argc, char **argv ) {
+ char* progname=argv[0];
+ string datasource,work_dir;
+
+ sglog().setLogLevels( SG_ALL, SG_INFO );
+
+ while (argc>1) {
+ if (!strcmp(argv[1],"--log-level")) {
+ if (argc<3) {
+ usage(progname);
+ }
+ setLoggingPriority(argv[2]);
+ argv+=2;
+ argc-=2;
+ } else if (!strcmp(argv[1],"--seperate-segments")) {
+ argv++;
+ argc--;
+ seperate_segments=1;
+ } else if (!strcmp(argv[1],"--continue-on-errors")) {
+ argv++;
+ argc--;
+ continue_on_errors=1;
+ } else if (!strcmp(argv[1],"--start-record")) {
+ if (argc<3) {
+ usage(progname);
+ }
+ start_record=atoi(argv[2]);
+ argv+=2;
+ argc-=2;
+ } else if (!strcmp(argv[1],"--where")) {
+ if (argc<3) {
+ usage(progname);
+ }
+ use_attribute_query=true;
+ attribute_query=argv[2];
+ argv+=2;
+ argc-=2;
+ } else if (!strcmp(argv[1],"--spat")) {
+ if (argc<6) {
+ usage(progname);
+ }
+ use_spatial_query=true;
+ spat_min_x=atof(argv[2]);
+ spat_min_y=atof(argv[3]);
+ spat_max_x=atof(argv[4]);
+ spat_max_y=atof(argv[5]);
+ argv+=5;
+ argc-=5;
+ } else if (!strcmp(argv[1],"--threads")) {
+ if (argc<3) {
+ usage(progname);
+ }
+ num_threads=atoi(argv[2]);
+ argv+=2;
+ argc-=2;
+ } else if (!strcmp(argv[1],"--all-threads")) {
+ num_threads=boost::thread::hardware_concurrency();
+ argv+=1;
+ argc-=1;
+ } else if (!strcmp(argv[1],"--debug")) {
+ argv++;
+ argc--;
+ save_shapefiles=true;
+ } else if (!strcmp(argv[1],"--help")) {
+ usage(progname);
+ } else {
+ break;
+ }
+ }
+
+ SG_LOG( SG_GENERAL, SG_ALERT, "\ncliff-decode version " << getTGVersion() );
+
+ if (argc<3) {
+ usage(progname);
+ }
+ work_dir=argv[1];
+ datasource=argv[2];
+
+ SGPath sgp( work_dir );
+ sgp.append( "dummy" );
+ sgp.create_dir( 0755 );
+
+ tgChopper results( work_dir );
+
+ // initialize persistant polygon counter
+ //string counter_file = work_dir + "/poly_counter";
+ //poly_index_init( counter_file );
+
+ // new chop api
+ //std::string counter_file2 = work_dir + "/poly_counter2";
+ //tgPolygon::ChopIdxInit( counter_file );
+
+ SG_LOG( SG_GENERAL, SG_DEBUG, "Opening datasource " << datasource << " for reading." );
+
+ GDALAllRegister();
+ GDALDataset *poDS;
+
+ poDS = (GDALDataset*) GDALOpenEx( datasource.c_str(), GDAL_OF_VECTOR, NULL, NULL, NULL );
+ if( poDS == NULL )
+ {
+ SG_LOG( SG_GENERAL, SG_ALERT, "Failed opening datasource " << datasource );
+ exit( 1 );
+ }
+
+ SG_LOG( SG_GENERAL, SG_ALERT, "Processing datasource " << datasource );
+
+ OGRLayer *poLayer;
+ if (argc>3) {
+ for (int i=3;i<argc;i++) {
+ poLayer = poDS->GetLayerByName( argv[i] );
+
+ if (poLayer == NULL )
+ {
+ SG_LOG( SG_GENERAL, SG_ALERT, "Failed opening layer " << argv[i] << " from datasource " << datasource );
+ exit( 1 );
+ }
+ processLayer(poLayer, results );
+ }
+ } else {
+ for (int i=0;i<poDS->GetLayerCount();i++) {
+ poLayer = poDS->GetLayer(i);
+
+ assert(poLayer != NULL);
+
+ processLayer(poLayer, results );
+ }
+ }
+
+ GDALClose(poDS);
+
+ SG_LOG(SG_GENERAL, SG_ALERT, "Saving to buckets");
+ results.Add_Extension("cliffs");
+ results.Save( save_shapefiles );
+ return 0;
+}
+
+
From 3f9c296bf4f91d9a00ec079b86fe90c25f2ce920 Mon Sep 17 00:00:00 2001
From: "James.Hester" <jxh@ansto.gov.au>
Date: Sun, 30 Dec 2018 16:25:56 +1100
Subject: [PATCH 10/15] Brought cliff-decode more closely in line with current
ogr-decode. Fixed bug introduced during merging.
---
src/Lib/terragear/tg_polygon_clip.cxx | 3 ++-
src/Prep/Cliff/cliff-decode.cxx | 27 ++++++++++++++-------------
2 files changed, 16 insertions(+), 14 deletions(-)
diff --git a/src/Lib/terragear/tg_polygon_clip.cxx b/src/Lib/terragear/tg_polygon_clip.cxx
index 791ca400..e97f9d4e 100644
--- a/src/Lib/terragear/tg_polygon_clip.cxx
+++ b/src/Lib/terragear/tg_polygon_clip.cxx
@@ -164,11 +164,12 @@ tgPolygon tgPolygon::Intersect( const tgPolygon& subject, const tgPolygon& clip)
c.Clear();
c.AddPaths(clipper_subject, ClipperLib::PolyType::Subject, subject.IsClosed());
c.AddPaths(clipper_clip, ClipperLib::PolyType::Clip, true);
+ if(subject.IsClosed()) {
c.Execute(ClipperLib::ClipType::Intersection, clipper_result, ClipperLib::PolyFillType::EvenOdd, ClipperLib::PolyFillType::EvenOdd);
result = tgPolygon::FromClipper( clipper_result );
}
else {
- c.Execute(ClipperLib::ctIntersection, clipper_tree_result, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd);
+ c.Execute(ClipperLib::ClipType::Intersection, clipper_tree_result, ClipperLib::PolyFillType::EvenOdd, ClipperLib::PolyFillType::EvenOdd);
result = tgPolygon::FromClipper( clipper_tree_result );
}
result = tgPolygon::AddColinearNodes( result, all_nodes );
diff --git a/src/Prep/Cliff/cliff-decode.cxx b/src/Prep/Cliff/cliff-decode.cxx
index 4b38298c..c60d31aa 100644
--- a/src/Prep/Cliff/cliff-decode.cxx
+++ b/src/Prep/Cliff/cliff-decode.cxx
@@ -3,9 +3,9 @@
// against and sorting
// them into the relevant tiles.
//
-// Written by James Hester starting 2018
+// Written by James Hester 2018
//
-// Based on ogr-decode by Ralf Gerlich.
+// Largely copied from ogr-decode by Ralf Gerlich.
//
// 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
@@ -21,6 +21,7 @@
// along with this program; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
//
+#include <chrono>
#include <string>
#include <map>
@@ -100,9 +101,9 @@ void Decoder::processLineString(OGRLineString* poGeometry)
// Do not close this contour
line.SetOpen(true);
- // clip the contour. Do not use usual clipper, as we wish to preserve the
- // whole contour, even if it goes outside the bounding box, as height
- // calculations along the boundary might be affected by it.
+ // clip the contour.
+ // TODO: don't bother clipping, as the contour is informational
+ // only, and simply output to all relevant buckets instead.
tgPolygon open_poly;
open_poly.SetOpen(); //do not try to close this one up
open_poly.AddContour(line);
@@ -299,6 +300,8 @@ int main( int argc, char **argv ) {
char* progname=argv[0];
string datasource,work_dir;
+ auto start_time = std::chrono::high_resolution_clock::now();
+
sglog().setLogLevels( SG_ALL, SG_INFO );
while (argc>1) {
@@ -379,14 +382,6 @@ int main( int argc, char **argv ) {
tgChopper results( work_dir );
- // initialize persistant polygon counter
- //string counter_file = work_dir + "/poly_counter";
- //poly_index_init( counter_file );
-
- // new chop api
- //std::string counter_file2 = work_dir + "/poly_counter2";
- //tgPolygon::ChopIdxInit( counter_file );
-
SG_LOG( SG_GENERAL, SG_DEBUG, "Opening datasource " << datasource << " for reading." );
GDALAllRegister();
@@ -428,6 +423,12 @@ int main( int argc, char **argv ) {
SG_LOG(SG_GENERAL, SG_ALERT, "Saving to buckets");
results.Add_Extension("cliffs");
results.Save( save_shapefiles );
+
+ auto finish_time = std::chrono::high_resolution_clock::now();
+ std::chrono::duration<double> elapsed = finish_time - start_time;
+ std::cout << std::endl << "Elapsed time: " << elapsed.count() << " seconds" << std::endl << std::endl;
+
+
return 0;
}
From 3192c91da62ead7e79b11c9a2b5e259627bd196e Mon Sep 17 00:00:00 2001
From: "James.Hester" <jxh@ansto.gov.au>
Date: Mon, 31 Dec 2018 00:11:55 +1100
Subject: [PATCH 11/15] Fixed missing return causing crash in Release mode,
changed debug messages to use standard flightgear tools, improved comments.
---
src/Lib/Array/array.cxx | 10 +++++-----
src/Lib/Array/array.hxx | 2 +-
src/Lib/Array/rectify_height.cxx | 28 +++++++++++++++++++++++++---
3 files changed, 31 insertions(+), 9 deletions(-)
diff --git a/src/Lib/Array/array.cxx b/src/Lib/Array/array.cxx
index 9c159d96..4ec2dbac 100644
--- a/src/Lib/Array/array.cxx
+++ b/src/Lib/Array/array.cxx
@@ -119,7 +119,7 @@ TGArray::close() {
//All polys in the bucket should be contours which we load
//into our contour list.
-bool TGArray::load_cliffs(const string & height_base)
+void TGArray::load_cliffs(const string & height_base)
{
//Get the directory so we can list the children
tgPolygon poly; //actually a contour but whatever...
@@ -467,7 +467,7 @@ void TGArray::rectify_heights(const double bad_zone) {
set_array_elev(xgrid,ygrid,(int) new_ht);
}
}
- std::cout << "Rectified " << rectified.size() << " points " << std::endl;
+ SG_LOG(SG_GENERAL, SG_DEBUG, "Rectified " << rectified.size() << " points ");
if(rectified.size()>0) {
for(auto r : rectified) {
bad_points.erase(std::remove(std::begin(bad_points),std::end(bad_points),r));
@@ -475,7 +475,7 @@ void TGArray::rectify_heights(const double bad_zone) {
rectified.clear();
} else {
if(bad_points.size() > 0) {
- std::cout << "Failed to rectify " << bad_points.size() << " points" << std::endl;
+ SG_LOG(SG_GENERAL, SG_DEBUG, "Failed to rectify " << bad_points.size() << " points");
}
break; // Cant do any more
}
@@ -536,7 +536,7 @@ double TGArray::rectify_point(const int xgrid, const int ygrid, const std::vecto
}
}
- if (pt == pt_cnt) { // perhaps we are in a bay, just take the
+ if (pt == pt_cnt) { // perhaps we have a concave cliff, just take the
// average of the known points
double totht = 0;
for(int pti = 0; pti <pt_cnt; pti++) {
@@ -555,7 +555,7 @@ double TGArray::rectify_point(const int xgrid, const int ygrid, const std::vecto
if (vert < -9000) vert = 0;
height = horiz + (vert - corner);
}
- std::cout << xgrid << "," << ygrid << ": was " << original_height << " , now " << height << std::endl;
+ SG_LOG(SG_GENERAL, SG_DEBUG, xgrid << "," << ygrid << ": was " << original_height << " , now " << height);
return height;
}
diff --git a/src/Lib/Array/array.hxx b/src/Lib/Array/array.hxx
index 9de84269..46fe8854 100644
--- a/src/Lib/Array/array.hxx
+++ b/src/Lib/Array/array.hxx
@@ -81,7 +81,7 @@ public:
bool open ( const std::string& file_base );
// Load contours from polygon files delineating height discontinuities
- bool load_cliffs(const std::string & height_base);
+ void load_cliffs(const std::string & height_base);
// return if array was successfully opened or not
bool is_open() const;
diff --git a/src/Lib/Array/rectify_height.cxx b/src/Lib/Array/rectify_height.cxx
index 06b35323..b6bf7709 100644
--- a/src/Lib/Array/rectify_height.cxx
+++ b/src/Lib/Array/rectify_height.cxx
@@ -1,4 +1,25 @@
-//rectify_height.cxx
+// rectify_height.cxx -- Correct height grid for incorrect heights near
+// discontinuities included in cliff files
+//
+// Written by James Hester 2018
+//
+// 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 <config.h>
#endif
@@ -16,7 +37,7 @@
/* 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. */
+the heights. */
// display usage and exit
static void usage( const std::string name ) {
@@ -38,7 +59,7 @@ int main( int argc, char **argv) {
std::string height_dir = "";
SGGeod min,max;
long tile_id = -1;
- double bad_zone = 30; //distance in m from cliff to rectify
+ double bad_zone = 100; //distance in m from cliff to rectify
sglog().setLogLevels(SG_ALL,SG_INFO);
sglog().set_log_priority(SG_DEBUG);
@@ -122,6 +143,7 @@ int main( int argc, char **argv) {
SG_LOG(SG_GENERAL,SG_DEBUG, "Failed to open array file " << array_path);
continue;
}
+ SG_LOG(SG_GENERAL,SG_DEBUG, "Successfully opened array file");
array.parse(bucket);
array.close();
array.remove_voids();
From e8b87fcf74b2417fbc1f89f06c762102034432c6 Mon Sep 17 00:00:00 2001
From: "James.Hester" <jxh@ansto.gov.au>
Date: Mon, 31 Dec 2018 10:12:49 +1100
Subject: [PATCH 12/15] Fix includes for full terragear tool suite, allow
arr.gz again for terrafit.
---
src/Airports/GenAirports850/CMakeLists.txt | 1 +
src/Prep/DemChop/CMakeLists.txt | 3 +++
src/Prep/TerraFit/terrafit.cc | 3 ++-
3 files changed, 6 insertions(+), 1 deletion(-)
diff --git a/src/Airports/GenAirports850/CMakeLists.txt b/src/Airports/GenAirports850/CMakeLists.txt
index 67722dff..a699042d 100644
--- a/src/Airports/GenAirports850/CMakeLists.txt
+++ b/src/Airports/GenAirports850/CMakeLists.txt
@@ -1,5 +1,6 @@
include_directories(${GDAL_INCLUDE_DIR})
include_directories(${PROJECT_SOURCE_DIR}/src/Lib)
+include_directories(${PROJECT_SOURCE_DIR}/src/Lib/terragear)
add_executable(genapts850
airport.hxx airport.cxx
diff --git a/src/Prep/DemChop/CMakeLists.txt b/src/Prep/DemChop/CMakeLists.txt
index 00d0ea37..1d29e622 100644
--- a/src/Prep/DemChop/CMakeLists.txt
+++ b/src/Prep/DemChop/CMakeLists.txt
@@ -1,3 +1,4 @@
+include_directories(${PROJECT_SOURCE_DIR}/src/Lib/terragear)
add_executable(demchop demchop.cxx)
@@ -49,6 +50,7 @@ endif(TIFF_FOUND)
add_executable(fillvoids fillvoids.cxx)
target_link_libraries(fillvoids
Array
+ terragear
${ZLIB_LIBRARY}
${SIMGEAR_CORE_LIBRARIES}
${SIMGEAR_CORE_LIBRARY_DEPENDENCIES})
@@ -58,6 +60,7 @@ install(TARGETS fillvoids RUNTIME DESTINATION bin)
add_executable(testassem testassem.cxx)
target_link_libraries(testassem
Array
+ 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 d84d671b..43f5cdf2 100644
--- a/src/Prep/TerraFit/terrafit.cc
+++ b/src/Prep/TerraFit/terrafit.cc
@@ -239,7 +239,8 @@ void walk_path(const SGPath& path) {
return;
}
- if ((path.lower_extension() == "arr") || (path.complete_lower_extension() == "arr.rectified.gz")) {
+ if ((path.lower_extension() == "arr") || (path.complete_lower_extension() == "arr.rectified.gz") ||
+ (path.complete_lower_extension() == "arr.gz")) {
SG_LOG(SG_GENERAL, SG_DEBUG, "will queue " << path);
queue_fit_file(path);
} else if (path.isDir()) {
From 5586be95f53b50ece061c5684fc3a52f41d1cde4 Mon Sep 17 00:00:00 2001
From: "James.Hester" <jxh@ansto.gov.au>
Date: Mon, 31 Dec 2018 11:16:16 +1100
Subject: [PATCH 13/15] Add description of cliff generation; remove debugging
printout.
---
README.howto | 54 ++++++++++++++++++++++++++++++++
src/Lib/terragear/tg_chopper.cxx | 10 +-----
2 files changed, 55 insertions(+), 9 deletions(-)
diff --git a/README.howto b/README.howto
index c313d35c..de0f8e16 100644
--- a/README.howto
+++ b/README.howto
@@ -1,5 +1,6 @@
Original version by Alexei Novikov <anovikov@heron.itep.ru> Sep. 3, 1999
Updates by Curtis Olson <http://www.flightgear.org/~curt>
+Cliff description by James Hester 2018
Overview
@@ -120,3 +121,56 @@ tg-construct-client in rude mode. After a while scenery is ready.
That's all folks,
Alexei.
+
+Addendum: Adding cliffs
+=======================
+
+Cliffs represent discontinuities in the height array described above. This
+is problematic because:
+(1) Most height grids will give an average height in the area around the
+grid
+(2) Interpolation between grid points to get the local height will convert
+a cliff into a much more gentle slope.
+
+To add cliffs into the scenery using OpenStreetMap data, the following
+approach can be used:
+
+(1) Extract cliffs from the OSM datafile using your favourite
+tool, for example for all cliffs near Sydney, Australia:
+
+osmosis --read-pbf australia-latest.osm.pbf --bounding-box top=-33 left=150
+bottom=-35 right=152 completeWays=yes --lp --tf accept-ways 'natural=cliff'
+--tf reject-relations --used-node --write-xml file='my_cliffs.osm'
+
+(2) Convert the XML above to Shapefiles:
+
+ogr2ogr -where "OGR_GEOMETRY='LineString'" -f 'ESRI shapefile' -lco SHPT=ARC
+cs_cliffs my_cliffs.osm
+
+(3) Run cliff-decode to place cliff information into the height grid working
+directory created previously (<data>/cs_cliffs in the example below):
+
+cliff-decode <work_base>/SRTM-3 <data>/cs_cliffs
+
+(4) Optionally run rectify_height to fix grid heights. This is recommended
+as otherwise cliffs will have uneven edges and bases as the grid points
+vary in distance from the cliff:
+
+rectify_height --work-dir=<work_base> --height-dir=SRTM-3 --min-lon=151.0
+--max-lon=152.0 --min-lat=-34.0 --max-lat=-33.0 --min-dist=100
+
+The min-dist parameter sets the distance from the cliff beyond which points
+are trustworthy. 100m is the default.
+
+(5) Place the cliffs into the scenery. This is identical to the process for
+rivers and roads and uses the same cliff information from step (2):
+
+ogr-decode --line-width 5 --area-type Rock <work_base>/Cliffs <data>/cs_cliffs"
+
+Note that an area-type of Cliffs will allow custom cliff texturing, if this
+effect is available. This command places a 5m wide strip of rock at the
+position of the cliffs. One side of this strip will have a calculated elevation
+at the top of the cliff, and the other will be at the bottom.
+
+(6) Run tg-construct as usual, including 'Rock' or 'Cliffs' polygon types to
+make sure that cliffs appear.
\ No newline at end of file
diff --git a/src/Lib/terragear/tg_chopper.cxx b/src/Lib/terragear/tg_chopper.cxx
index 25cc71e4..a57fe34d 100644
--- a/src/Lib/terragear/tg_chopper.cxx
+++ b/src/Lib/terragear/tg_chopper.cxx
@@ -24,15 +24,7 @@ tgPolygon tgChopper::Clip( const tgPolygon& subject,
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 );
From dfe81ce9fa0e68ddb025983f8f85123135bfca7a Mon Sep 17 00:00:00 2001
From: "James.Hester" <jxh@ansto.gov.au>
Date: Sat, 12 Jan 2019 09:29:39 +1100
Subject: [PATCH 14/15] Changed exit(1) to EXIT_FAILURE; use delegation in
constructors; check for endianness only once per invocation.
---
src/Lib/HGT/dted.cxx | 14 ++++++--------
1 file changed, 6 insertions(+), 8 deletions(-)
diff --git a/src/Lib/HGT/dted.cxx b/src/Lib/HGT/dted.cxx
index 05d7de53..56700c08 100644
--- a/src/Lib/HGT/dted.cxx
+++ b/src/Lib/HGT/dted.cxx
@@ -64,11 +64,8 @@ TGDted::TGDted( int _res )
}
-TGDted::TGDted( int _res, const SGPath &file )
+TGDted::TGDted( int _res, const SGPath &file ):TGDted(_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::open( file );
}
@@ -109,7 +106,7 @@ TGDted::open ( const SGPath &f ) {
cout << "Proceeding with " << file_name.str() << endl;
} else {
SG_LOG(SG_GENERAL, SG_ALERT, "Failed to issue system call " << command );
- exit(1);
+ return EXIT_FAILURE;
}
}
@@ -168,6 +165,7 @@ TGDted::close () {
bool
TGDted::load( ) {
int size;
+ bool little_endian = sgIsLittleEndian();
if ( dted_resolution == 1 ) {
cols = rows = size = 3601;
col_step = row_step = 1;
@@ -179,7 +177,7 @@ TGDted::load( ) {
cout << " are supported!" << endl;
return false;
}
- if (sgIsLittleEndian()) {
+ if (little_endian) {
cout << "Little Endian: swapping input values" << endl;
}
@@ -200,7 +198,7 @@ TGDted::load( ) {
gzread(fd,&dummy,3); //block count
gzread(fd,&longct,2); //Longitude count
gzread(fd,&latct,2); //Latitude count
- if ( sgIsLittleEndian() ) {
+ if ( little_endian ) {
sgEndianSwap(&longct);
}
// cout << "Longitude count " << longct << endl;
@@ -209,7 +207,7 @@ TGDted::load( ) {
if ( gzread ( fd, var, 2 ) != sizeof(short) ) {
return false;
}
- if ( sgIsLittleEndian() ) {
+ if ( little_endian ) {
sgEndianSwap( (unsigned short int*)var);
}
}
From 3421ce2018611e93d434afba0845ee277fdb297d Mon Sep 17 00:00:00 2001
From: "James.Hester" <jxh@ansto.gov.au>
Date: Sun, 13 Jan 2019 09:59:00 +1100
Subject: [PATCH 15/15] Numerous small fixes following code review:
(1) Whitespace use matched better to surrounding code
(2) While(1) in TGArray::rectify_heights replaced with do/while loop
(3) Incorrect interpolation for 2 points fixed in altitude_from_grid
(4) exit(1) replaced with return EXIT_FAILURE in cliff-decode.cxx
---
src/Lib/Array/array.cxx | 595 ++++++++++++++------------
src/Lib/terragear/tg_chopper.cxx | 11 +-
src/Lib/terragear/tg_contour.cxx | 138 +++---
src/Lib/terragear/tg_polygon_clip.cxx | 8 +-
src/Prep/Cliff/cliff-decode.cxx | 4 +-
5 files changed, 404 insertions(+), 352 deletions(-)
diff --git a/src/Lib/Array/array.cxx b/src/Lib/Array/array.cxx
index 4ec2dbac..340eb275 100644
--- a/src/Lib/Array/array.cxx
+++ b/src/Lib/Array/array.cxx
@@ -67,14 +67,14 @@ bool TGArray::open( const string& file_base ) {
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;
- }
+ // 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);
@@ -121,34 +121,34 @@ TGArray::close() {
void 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.lower_extension();
- if (lext == "cliffs") {
- 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<count; i++ ) {
- poly.LoadFromGzFile( fp );
- if ( poly.Contours()==1 ) { //should always have one contour
- cliffs_list.push_back(poly.GetContour(0));
- } else {
- SG_LOG( SG_GENERAL, SG_WARN, " Found " << poly.Contours() << " contours in " << p.realpath() );
+ //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);
+ for (const SGPath& p: files) {
+ if (p.file_base() != b.file_base()) {
+ continue;
+ }
+
+ string lext = p.lower_extension();
+ if (lext == "cliffs") {
+ 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<count; i++ ) {
+ poly.LoadFromGzFile( fp );
+ if ( poly.Contours()==1 ) { //should always have one contour
+ cliffs_list.push_back(poly.GetContour(0));
+ } else {
+ SG_LOG( SG_GENERAL, SG_WARN, " Found " << poly.Contours() << " contours in " << p.realpath() );
+ }
+ }
}
- }
}
- }
}
@@ -408,7 +408,7 @@ double TGArray::closest_nonvoid_elev( double lon, double lat ) const {
for ( int row = 0; row < rows; row++ ) {
for ( int col = 0; col < cols; col++ ) {
- SGGeod p1 = SGGeod::fromDeg( (originx + col * col_step)/3600.0, (originy + row * row_step)/3600.0 );
+ SGGeod p1 = SGGeod::fromDeg( (originx + col * col_step)/3600.0, (originy + row * row_step)/3600.0 );
double dist = SGGeodesy::distanceM( p0, p1 );
double elev = get_array_elev(col, row);
if ( dist < mindist && elev > -9000 ) {
@@ -425,61 +425,68 @@ double TGArray::closest_nonvoid_elev( double lon, double lat ) const {
}
}
+//Find and remember all points that are bad because they are
+//too close to a cliff
std::vector<int> 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<int> bad_points; //local to avoid multi-thread issues
- for(int horiz=0;horiz<cols;horiz++) {
- double lon = (originx + col_step*horiz)/3600;
- for(int vert=0;vert<rows;vert++) {
- double lat = (originy + row_step*vert)/3600;
- if(is_near_cliff(lon,lat,bad_zone)) {
- bad_points.push_back(horiz+vert*cols);
+
+ std::vector<int> bad_points; //local to avoid multi-thread issues
+
+ for( int horiz=0;horiz<cols;horiz++ ) {
+ double lon = (originx + col_step*horiz)/3600;
+ for( int vert=0;vert<rows;vert++ ) {
+ double lat = (originy + row_step*vert)/3600;
+ if( is_near_cliff(lon,lat,bad_zone) ) {
+ bad_points.push_back(horiz+vert*cols);
+ }
}
}
- }
- return bad_points;
+
+ return bad_points;
}
// Check to see if the specified grid point is bad
bool TGArray::is_bad_point(const int xgrid, const int ygrid, const std::vector<int> 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;
- }
+ 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(const double bad_zone) {
- double new_ht;
- std::vector<int> rectified,bad_points;
- bad_points = collect_bad_points(bad_zone);
- 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);
- }
- }
- SG_LOG(SG_GENERAL, SG_DEBUG, "Rectified " << rectified.size() << " points ");
- 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 {
- if(bad_points.size() > 0) {
+void TGArray::rectify_heights( const double bad_zone ) {
+ double new_ht;
+ std::vector<int> rectified,bad_points;
+ int total_rectified;
+ bad_points = collect_bad_points( bad_zone );
+
+ do {
+ 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 );
+ }
+ }
+ total_rectified = rectified.size();
+ SG_LOG(SG_GENERAL, SG_DEBUG, "Rectified " << total_rectified << " points ");
+
+ if( total_rectified > 0 ) {
+ for( auto r : rectified ) {
+ bad_points.erase( std::remove( std::begin(bad_points), std::end(bad_points),r) );
+ }
+ rectified.clear();
+ }
+ } while ( total_rectified > 0 );
+
+ if( bad_points.size() > 0 ) {
SG_LOG(SG_GENERAL, SG_DEBUG, "Failed to rectify " << bad_points.size() << " points");
}
- break; // Cant do any more
- }
- }
+
}
/* If we have cliffs, it is possible that a grid point will be too close
@@ -496,67 +503,82 @@ through the three known points.
* * *
TODO: Handle points on the boundaries. */
+
double TGArray::rectify_point(const int xgrid, const int ygrid, const std::vector<int> 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) {
- 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
- corners[pt_cnt][0] = horiz;
- corners[pt_cnt][1] = vert;
- pt_cnt++;
+ //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 ) {
+ 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
+
+ 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;
- 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)) {
- 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;
+ } // 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;
+ 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 ) ) {
+
+ 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 have a concave cliff, just take the
- // average of the known points
- double totht = 0;
- for(int pti = 0; pti <pt_cnt; pti++) {
- totht = totht + get_array_elev(xgrid+corners[pti][0],ygrid+corners[pti][1]);
- }
- height = totht/pt_cnt;
- } else {
+ if (pt == pt_cnt) {
+ // perhaps we have a concave cliff, just take the
+ // average of the known points
+ double totht = 0;
+ for( int pti = 0; pti <pt_cnt; pti++ ) {
+ totht = totht + get_array_elev( xgrid+corners[pti][0],ygrid+corners[pti][1] );
+ }
+
+ height = totht/pt_cnt;
+
+ } else {
- // We have three points, calculate the height
- // Set anything very negative to zero
- 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);
- if (corner < -9000) corner = 0;
- if (horiz < -9000) horiz = 0;
- if (vert < -9000) vert = 0;
- height = horiz + (vert - corner);
- }
- SG_LOG(SG_GENERAL, SG_DEBUG, xgrid << "," << ygrid << ": was " << original_height << " , now " << height);
- return height;
+ // We have three points, calculate the height
+ // Set anything very negative to zero
+ 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 );
+ if ( corner < -9000 ) corner = 0;
+ if ( horiz < -9000 ) horiz = 0;
+ if ( vert < -9000 ) vert = 0;
+ height = horiz + ( vert - corner );
+ }
+
+ SG_LOG(SG_GENERAL, SG_DEBUG, xgrid << "," << ygrid << ": was " << original_height << " , now " << height);
+
+ return height;
}
// return the current altitude based on grid data.
@@ -578,7 +600,7 @@ double TGArray::altitude_from_grid( double lon, double lat ) const {
------
then calculate our end points
- */
+ */
// Store in degrees for later
double londeg = lon/3600;
@@ -589,8 +611,6 @@ double TGArray::altitude_from_grid( double lon, double lat ) const {
xindex = (int)(xlocal);
yindex = (int)(ylocal);
- // printf("xindex = %d yindex = %d\n", xindex, yindex);
-
if ( xindex + 1 == cols ) {
xindex--;
}
@@ -601,7 +621,9 @@ double TGArray::altitude_from_grid( double lon, double lat ) const {
if ( (xindex < 0) || (xindex + 1 >= cols) ||
(yindex < 0) || (yindex + 1 >= rows) ) {
+
SG_LOG(SG_GENERAL, SG_DEBUG, "WARNING: Attempt to interpolate value outside of array!!!" );
+
return -9999;
}
@@ -613,95 +635,101 @@ double TGArray::altitude_from_grid( double lon, double lat ) const {
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
+ 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) {
+
+ switch ( missing ) {
case 0: //SW corner missing
- x1 = corners[0][0];
- y1 = corners[0][1];
+ x1 = corners[0][0];
+ y1 = corners[0][1];
- x2 = corners[1][0];
- y2 = corners[1][1];
+ x2 = corners[1][0];
+ y2 = corners[1][1];
- x3 = corners[2][0];
- y3 = corners[2][1];
+ x3 = corners[2][0];
+ y3 = corners[2][1];
- dy = 1 - dy;
- break;
+ dy = 1 - dy;
+ break;
case 1: //NW corner missing
- x1 = corners[0][0];
- y1 = corners[0][1];
+ x1 = corners[0][0];
+ y1 = corners[0][1];
- x2 = corners[2][0];
- y2 = corners[2][1];
+ x2 = corners[2][0];
+ y2 = corners[2][1];
- x3 = corners[1][0];
- y3 = corners[1][1];
+ x3 = corners[1][0];
+ y3 = corners[1][1];
- break;
+ break;
case 2: //NE corner missing
- x1 = corners[2][0];
- y1 = corners[2][1];
+ x1 = corners[2][0];
+ y1 = corners[2][1];
- x2 = corners[0][0];
- y2 = corners[0][1];
+ x2 = corners[0][0];
+ y2 = corners[0][1];
- x3 = corners[1][0];
- y3 = corners[1][1];
+ x3 = corners[1][0];
+ y3 = corners[1][1];
- dx = 1 - dx; //x1 is SE corner
- break;
+ dx = 1 - dx; //x1 is SE corner
+ break;
case 3: //SE corner missing
- x1 = corners[2][0];
- y1 = corners[2][1];
+ x1 = corners[2][0];
+ y1 = corners[2][1];
- x2 = corners[1][0];
- y2 = corners[1][1];
+ x2 = corners[1][0];
+ y2 = corners[1][1];
- x3 = corners[0][0];
- y3 = corners[0][1];
+ x3 = corners[0][0];
+ y3 = corners[0][1];
- dx = 1 - dx; //x1 is NE corner
- dy = 1 - dy;
- break;
+ 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);
@@ -709,13 +737,14 @@ double TGArray::altitude_from_grid( double lon, double lat ) const {
zB = dx * (z3 - z1) + z1;
if ( dx > SG_EPSILON ) {
- elev = dy * (zB - zA) / dx + zA;
+ elev = dy * (zB - zA) / dx + zA;
} else {
- elev = zA;
+ elev = zA;
}
break;
- case 2: //project onto line connecting two points
+
+ case 2: //project onto line connecting two points
x1 = corners[0][0];
y1 = corners[0][1];
z1 = get_array_elev(x1,y1);
@@ -729,119 +758,135 @@ double TGArray::altitude_from_grid( double lon, double lat ) const {
dx = xlocal - x1;
dy = ylocal - y1;
if (x1==x2) {
- elev = z1+dy*(z2-z1);
+ elev = z1+dy*(z2-z1)/(y2-y1);
}
else if (y1==y2) {
- elev = z1+dx*(z2-z1);
+ elev = z1+dx*(z2-z1)/(x2-x1);
}
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);
+ 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
+
+ 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 " << std::setprecision(10) << londeg << "," << latdeg );
SG_LOG(SG_GENERAL, SG_WARN, "Grid points ("<< std::setprecision(9) << lon1 << "," << lat1 << "),("<<lon2<<","<<lat2<<")");
- default: // all corners
+
+ default: // all corners
dx = xlocal - xindex;
dy = ylocal - yindex;
if ( dx > dy ) {
- // lower triangle
+ // lower triangle
- x1 = xindex;
- y1 = yindex;
- z1 = get_array_elev(x1, y1);
+ x1 = xindex;
+ y1 = yindex;
+ z1 = get_array_elev(x1, y1);
- x2 = xindex + 1;
- y2 = yindex;
- z2 = get_array_elev(x2, y2);
+ x2 = xindex + 1;
+ y2 = yindex;
+ 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 ) {
- // don't interpolate off a void
- return closest_nonvoid_elev( lon, lat );
- }
+ 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 ) {
- elev = dy * (zB - zA) / dx + zA;
- } else {
- elev = zA;
- }
+ if ( dx > SG_EPSILON ) {
+ elev = dy * (zB - zA) / dx + zA;
+ } else {
+ elev = zA;
+ }
} else {
- // upper triangle
+ // 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 ) {
- // don't interpolate off a void
- return closest_nonvoid_elev( lon, lat );
- }
+ 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 ) {
- elev = dx * (zB - zA) / dy + zA;
- } else {
- elev = zA;
- }
+ if ( dy > SG_EPSILON ) {
+ elev = dx * (zB - zA) / dy + zA;
+ } 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)<SG_EPSILON && fabs(lat1-lat2)<SG_EPSILON) return true;
- SGGeod pt1 = SGGeod::fromDeg(lon1,lat1);
- SGGeod pt2 = SGGeod::fromDeg(lon2,lat2);
- bool same_side = true;
- for (int i=0;i<cliffs_list.size();i++) {
- bool check_result = cliffs_list[i].AreSameSide(pt1,pt2);
- if(!check_result) {
- SG_LOG(SG_GENERAL, SG_DEBUG, "Cliff " << i <<":" <<pt1 << " and " << pt2 << " on opposite sides");
- same_side = false;
- break;
+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)<SG_EPSILON && fabs(lat1-lat2)<SG_EPSILON ) return true;
+
+ SGGeod pt1 = SGGeod::fromDeg( lon1,lat1 );
+ SGGeod pt2 = SGGeod::fromDeg( lon2,lat2 );
+ bool same_side = true;
+
+ for ( int i=0;i<cliffs_list.size();i++ ) {
+ bool check_result = cliffs_list[i].AreSameSide( pt1,pt2 );
+
+ if(!check_result) {
+
+ SG_LOG(SG_GENERAL, SG_DEBUG, "Cliff " << i <<":" <<pt1 << " and " << pt2 << " on opposite sides");
+
+ same_side = false;
+ break;
+ }
}
- }
- return same_side;
+
+ 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 double bad_zone) const {
- if (cliffs_list.size()==0) return false;
- SGGeod pt1 = SGGeod::fromDeg(lon1,lat1);
- for (int i=0;i<cliffs_list.size();i++) {
- double dist = cliffs_list[i].MinDist(pt1);
- if (dist < bad_zone) return true;
- }
- return false;
+bool TGArray::is_near_cliff( const double lon1, const double lat1, const double bad_zone ) const {
+
+ if (cliffs_list.size()==0) return false;
+
+ SGGeod pt1 = SGGeod::fromDeg(lon1,lat1);
+
+ for ( int i=0;i<cliffs_list.size();i++ ) {
+ double dist = cliffs_list[i].MinDist(pt1);
+ if (dist < bad_zone) return true;
+ }
+
+ return false;
}
TGArray::~TGArray( void )
@@ -875,10 +920,10 @@ void TGArray::set_array_elev( int col, int row, int val )
bool TGArray::is_open() const
{
- if ( array_in != NULL ) {
- return true;
- } else {
- return false;
- }
+ if ( array_in != NULL ) {
+ return true;
+ } else {
+ return false;
+ }
}
diff --git a/src/Lib/terragear/tg_chopper.cxx b/src/Lib/terragear/tg_chopper.cxx
index a57fe34d..0206373b 100644
--- a/src/Lib/terragear/tg_chopper.cxx
+++ b/src/Lib/terragear/tg_chopper.cxx
@@ -11,7 +11,7 @@
#include "tg_misc.hxx"
tgPolygon tgChopper::Clip( const tgPolygon& subject,
- const std::string& type, SGBucket& b)
+ const std::string& type, SGBucket& b )
{
tgPolygon base, result;
@@ -21,8 +21,7 @@ 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);
- // Debug: See if numbers of nodes have changed
+ result = tgPolygon::Intersect( subject, base );
if ( result.Contours() > 0 ) {
if ( subject.GetPreserve3D() ) {
@@ -68,7 +67,7 @@ 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 )
@@ -95,7 +94,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
{
@@ -122,7 +121,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() ) {
diff --git a/src/Lib/terragear/tg_contour.cxx b/src/Lib/terragear/tg_contour.cxx
index 2d5ac262..b8497b4c 100644
--- a/src/Lib/terragear/tg_contour.cxx
+++ b/src/Lib/terragear/tg_contour.cxx
@@ -61,7 +61,7 @@ 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<node_list.size(); i++) {
a = SGGeod_ToSGVec2d( node_list[i] );
@@ -77,80 +77,88 @@ double tgContour::GetArea( void ) const
// 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 = x2-x1;
- double ydif = y2-y1;
- /*We describe a line parametrically:
+ //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 = x2-x1;
+ double ydif = y2-y1;
+
+ /*We describe a line parametrically:
- x1 (x2-x1)
- L = + t
- y1 (y2-y1)
+ 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.
+ 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).
+ 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
+ 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)
+ 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;
- if (node_list.size()) {
- int j = node_list.size() - 1;
- for (int i=0;i<node_list.size()-1;i++) {
- double nx1 = node_list[i].getLatitudeDeg();
- double ny1 = node_list[i].getLongitudeDeg();
- double nx2 = node_list[i+1].getLatitudeDeg();
- double ny2 = node_list[i+1].getLongitudeDeg();
- double nydif = ny2-ny1;
- double nxdif = nx2-nx1;
- double denom = xdif*nydif - ydif*nxdif;
- if (denom != 0) { //Not parallel
- 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 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++;
- }
+ 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;
+ if (node_list.size()) {
+ int j = node_list.size() - 1;
+ for (int i=0;i<node_list.size()-1;i++) {
+ double nx1 = node_list[i].getLatitudeDeg();
+ double ny1 = node_list[i].getLongitudeDeg();
+ double nx2 = node_list[i+1].getLatitudeDeg();
+ double ny2 = node_list[i+1].getLongitudeDeg();
+ double nydif = ny2-ny1;
+ double nxdif = nx2-nx1;
+ double denom = xdif*nydif - ydif*nxdif;
+
+ if (denom != 0) { //Not parallel
+ 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 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++;
+ }
+ }
}
- }
- bool isinter = (intersect_ct%2 == 0);
- return isinter;
+
+ bool isinter = (intersect_ct%2 == 0);
+ return isinter;
}
-double tgContour::MinDist(const SGGeod& probe) const
-{
- SGVec3d probexyz;
- SGGeodesy::SGGeodToCart(probe,probexyz);
- double mindist = 100000.0;
- double dist;
- if (node_list.size()) {
- int j = node_list.size() - 1;
- for (int i=0;i<j;i++) {
- SGVec3d start,end;
- SGGeodesy::SGGeodToCart(node_list[i],start);
- SGGeodesy::SGGeodToCart(node_list[i+1],end);
- SGLineSegment<double> piece = SGLineSegment<double>(start,end);
- dist = distSqr(piece,probexyz);
- if (dist < mindist) mindist = dist;
+double tgContour::MinDist(const SGGeod& probe) const {
+ SGVec3d probexyz;
+ SGGeodesy::SGGeodToCart( probe,probexyz );
+ double mindist = 100000.0;
+ double dist;
+
+ if ( node_list.size() ) {
+
+ int j = node_list.size() - 1;
+
+ for (int i=0;i<j;i++) {
+ SGVec3d start,end;
+ SGGeodesy::SGGeodToCart( node_list[i],start );
+ SGGeodesy::SGGeodToCart( node_list[i+1],end );
+ SGLineSegment<double> piece = SGLineSegment<double>(start,end);
+ dist = distSqr( piece,probexyz );
+ if (dist < mindist) mindist = dist;
+ }
}
- }
- return sqrt(mindist);
+
+ return sqrt(mindist);
}
bool tgContour::IsInside( const tgContour& inside, const tgContour& outside )
diff --git a/src/Lib/terragear/tg_polygon_clip.cxx b/src/Lib/terragear/tg_polygon_clip.cxx
index e97f9d4e..1c3de84d 100644
--- a/src/Lib/terragear/tg_polygon_clip.cxx
+++ b/src/Lib/terragear/tg_polygon_clip.cxx
@@ -137,7 +137,7 @@ tgPolygon tgPolygon::Diff( const tgPolygon& subject, tgPolygon& clip )
}
//Intersect will keep open paths open
-tgPolygon tgPolygon::Intersect( const tgPolygon& subject, const tgPolygon& clip)
+tgPolygon tgPolygon::Intersect( const tgPolygon& subject, const tgPolygon& clip )
{
tgPolygon result;
UniqueSGGeodSet all_nodes;
@@ -162,14 +162,14 @@ tgPolygon tgPolygon::Intersect( const tgPolygon& subject, const tgPolygon& clip)
ClipperLib::Clipper c;
c.Clear();
- c.AddPaths(clipper_subject, ClipperLib::PolyType::Subject, subject.IsClosed());
+ c.AddPaths(clipper_subject, ClipperLib::PolyType::Subject, subject.IsClosed() );
c.AddPaths(clipper_clip, ClipperLib::PolyType::Clip, true);
if(subject.IsClosed()) {
- c.Execute(ClipperLib::ClipType::Intersection, clipper_result, ClipperLib::PolyFillType::EvenOdd, ClipperLib::PolyFillType::EvenOdd);
+ c.Execute(ClipperLib::ClipType::Intersection, clipper_result, ClipperLib::PolyFillType::EvenOdd, ClipperLib::PolyFillType::EvenOdd );
result = tgPolygon::FromClipper( clipper_result );
}
else {
- c.Execute(ClipperLib::ClipType::Intersection, clipper_tree_result, ClipperLib::PolyFillType::EvenOdd, ClipperLib::PolyFillType::EvenOdd);
+ c.Execute(ClipperLib::ClipType::Intersection, clipper_tree_result, ClipperLib::PolyFillType::EvenOdd, ClipperLib::PolyFillType::EvenOdd );
result = tgPolygon::FromClipper( clipper_tree_result );
}
result = tgPolygon::AddColinearNodes( result, all_nodes );
diff --git a/src/Prep/Cliff/cliff-decode.cxx b/src/Prep/Cliff/cliff-decode.cxx
index c60d31aa..5d398263 100644
--- a/src/Prep/Cliff/cliff-decode.cxx
+++ b/src/Prep/Cliff/cliff-decode.cxx
@@ -391,7 +391,7 @@ int main( int argc, char **argv ) {
if( poDS == NULL )
{
SG_LOG( SG_GENERAL, SG_ALERT, "Failed opening datasource " << datasource );
- exit( 1 );
+ return EXIT_FAILURE;
}
SG_LOG( SG_GENERAL, SG_ALERT, "Processing datasource " << datasource );
@@ -404,7 +404,7 @@ int main( int argc, char **argv ) {
if (poLayer == NULL )
{
SG_LOG( SG_GENERAL, SG_ALERT, "Failed opening layer " << argv[i] << " from datasource " << datasource );
- exit( 1 );
+ return EXIT_FAILURE;
}
processLayer(poLayer, results );
}