Merge branch 'cliffs' of ../../stable/terragear into next

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)

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

src/BuildTiles/Main/default_priorities.txt

@@ -30,6 +30,7 @@ River			stream
 IntermittentStream	stream
 Watercourse		stream
 Canal			stream
+Cliffs                  cliff   # A cliff face
 Glacier			other	# Solid ice/snow
 PackIce			other	# Water with ice packs
 PolarIce		other
@@ -85,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

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();
     }

src/BuildTiles/Main/tgconstruct_poly.cxx

@@ -65,8 +65,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") || (lext == "arr.rectified.gz") || (lext == "arr.new.gz"))
             {
                 // skipped!
             } else {

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
@@ -12,13 +15,23 @@ 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
 	${ZLIB_LIBRARY}
 	${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 ()

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>
@@ -56,15 +57,27 @@ TGArray::TGArray( const string &file ):
 
 
 // open an Array file (and fitted 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 );
@@ -79,7 +92,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 +115,43 @@ 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.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() );
+        }
+      }
+    }
+  }
+}
+
+  
 void
 TGArray::unload( void ) {
     if (array_in) {
@@ -197,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
@@ -230,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
@@ -316,14 +400,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 ) {
@@ -340,6 +425,139 @@ double TGArray::closest_nonvoid_elev( double lon, double lat ) const {
     }
 }
 
+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);
+        }
+    }
+  }
+  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(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);
+      }
+  }
+  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 {
+    if(bad_points.size() > 0) {
+    std::cout << "Failed to rectify " << bad_points.size() << " points" << std::endl;
+    }
+    break;   // Cant do any more
+  }
+  }
+}
+
+/* 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) {
+    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;
+        }
+    }
+  
+  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);
+  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;
+}
   
 // return the current altitude based on grid data.
 // TODO: We should rewrite this to interpolate exact values, but for now this is good enough
@@ -362,6 +580,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,6 +605,150 @@ double TGArray::altitude_from_grid( double lon, double lat ) const {
 	return -9999;
     }
 
+    // Now check if we are on the same side of any cliffs
+
+    // 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];
+
+          x2 = corners[1][0];
+          y2 = corners[1][1];
+
+          x3 = corners[2][0];
+          y3 = corners[2][1];
+
+          dy = 1 - dy;
+          break;
+        case 1:                 //NW corner missing
+          x1 = corners[0][0];  
+          y1 = corners[0][1];
+
+          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 " << 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;
 
@@ -444,10 +809,40 @@ double TGArray::altitude_from_grid( double lon, double lat ) const {
 	    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;
+}
+
+//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;
+}
       
 TGArray::~TGArray( void )
 {
@@ -486,3 +881,4 @@ bool TGArray::is_open() const
       return false;
   }
 }
+

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 {
 
@@ -42,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;
 
@@ -52,7 +59,15 @@ private:
     std::vector<SGGeod> corner_list;
     std::vector<SGGeod> fitted_list;
 
+  // list of cliff contours
+  tgcontour_list cliffs_list;
     void parse_bin();
+
+  // Routines for height rectification
+  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 double bad_zone) const;
 public:
 
     // Constructor
@@ -65,6 +80,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;
 
@@ -77,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;
 
@@ -103,6 +127,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 );
 };

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);
+  }
+}

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>
@@ -43,6 +44,9 @@ int main( int argc, char **argv ) {
 
     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;
         give_help( argv[0] );

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
 )

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",&degrees,&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",&degrees,&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;
+}

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
+
+

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} )
 

src/Lib/terragear/tg_chopper.cxx

@@ -11,8 +11,7 @@
 #include "tg_misc.hxx"
 
 tgPolygon tgChopper::Clip( const tgPolygon& subject,
-                      const std::string& type,
+                      const std::string& type, SGBucket& b)
-                      SGBucket& b )
 {
     tgPolygon base, result;
 
@@ -23,7 +22,17 @@ tgPolygon tgChopper::Clip( const tgPolygon& subject,
     base.AddNode( 0, b.get_corner( SG_BUCKET_NW ) );
 
     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 );
@@ -34,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 );
@@ -63,6 +75,7 @@ void tgChopper::ClipRow( const tgPolygon& subject, const double& center_lat, con
     }
 }
 
+
 void tgChopper::Add( const tgPolygon& subject, const std::string& type)
 {
     // bail out immediately if polygon is empty
@@ -205,7 +218,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 ) {

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
 };

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>
 
@@ -65,7 +66,6 @@ double tgContour::GetArea( void ) const
         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;
         }
@@ -74,6 +74,85 @@ 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 = x2-x1;
+  double ydif = y2-y1;
+  /*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.
+
+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;
+  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;
+}
+
+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
@@ -424,6 +503,7 @@ tgContour tgContour::FromClipper( const ClipperLib::Path& subject )
     return result;
 }
 
+
 tgRectangle tgContour::GetBoundingBox( void ) const
 {
     SGGeod min, max;
@@ -803,14 +883,16 @@ tgContour tgContour::AddColinearNodes( const tgContour& subject, std::vector<SGG
     p0 = subject.GetNode( subject.GetSize() - 1 );
     p1 = subject.GetNode( 0 );
 
+    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 flag setting
+    // maintain original hole and openness flag setting
     result.SetHole( subject.GetHole() );
+    result.SetOpen( subject.GetOpen() );
 
     return result;
 }
@@ -834,14 +916,16 @@ tgContour tgContour::AddColinearNodes( const tgContour& subject, bool preserve3d
     p0 = subject.GetNode( subject.GetSize() - 1 );
     p1 = subject.GetNode( 0 );
 
+    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 flag setting
+    // maintain original hole and open flag settings
     result.SetHole( subject.GetHole() );
+    result.SetOpen( subject.GetOpen() );
     
     return result;
 }

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"
@@ -29,6 +31,7 @@ class tgContour
 public:
     tgContour() {
         hole = false;
+        keep_open = false;
     }
 
     void Erase() {
@@ -38,10 +41,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();
     }
@@ -101,6 +113,10 @@ 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 );
@@ -135,6 +151,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;

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

src/Lib/terragear/tg_polygon_clip.cxx

@@ -136,6 +136,7 @@ tgPolygon tgPolygon::Diff( const tgPolygon& subject, tgPolygon& clip )
     return result;
 }
 
+//Intersect will keep open paths open
 tgPolygon tgPolygon::Intersect( const tgPolygon& subject, const tgPolygon& clip)
 {
     tgPolygon result;
@@ -156,22 +157,30 @@ 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::PolyType::Subject, true);
+    c.AddPaths(clipper_subject, ClipperLib::PolyType::Subject, subject.IsClosed());
     c.AddPaths(clipper_clip, ClipperLib::PolyType::Clip, true);
     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);
+      result = tgPolygon::FromClipper( clipper_tree_result );
+    }
     result = tgPolygon::AddColinearNodes( result, all_nodes );
-
     result.SetMaterial( subject.GetMaterial() );
     result.SetTexParams( subject.GetTexParams() );
     result.SetId( subject.GetId() );
     result.SetPreserve3D( subject.GetPreserve3D() );
-    
+    if(subject.IsClosed()) {
+      result.SetClosed();
+    } else {
+      result.SetOpen();
+    }
     return result;
 }
 
@@ -199,6 +208,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;

src/Prep/CMakeLists.txt

@@ -5,3 +5,4 @@ add_subdirectory(DemChop)
 add_subdirectory(Terra)
 add_subdirectory(TerraFit)
 add_subdirectory(OGRDecode)
+add_subdirectory(Cliff)

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} )

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;
+}

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})

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()) {