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

@@ -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
+Cliff                   cliff   # A cliff face
 Glacier			other	# Solid ice/snow
 PackIce			other	# Water with ice packs
 PolarIce		other

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
+#endif // _PRIORITIES_HXX

src/BuildTiles/Main/tgconstruct_elevation.cxx

@@ -222,4 +222,4 @@ void TGConstruct::CalcElevations( void )
             }
         }
     }
-}
+}

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

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

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

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;

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

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

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
+#endif // _TGCONTOUR_HXX

src/Prep/CMakeLists.txt

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