Added support for including height information in the intermediate polygon

format.  The big trick was that the polygon clipper is completely 2d.  So
I needed to add code to preserve the elevations in the clipped output and
fill in plausible elevations for any new points created as a result of the
clipping.

src/BuildTiles/Clipper/Makefile.am

@@ -11,6 +11,7 @@ testclipper_LDADD = \
 	$(top_builddir)/src/BuildTiles/Osgb36/libOsgb36.a \
 	$(top_builddir)/src/BuildTiles/Triangulate/libTriangulate.a \
 	$(top_builddir)/src/Lib/Polygon/libPolygon.a \
+	$(top_builddir)/src/Lib/Geometry/libGeometry.a \
 	$(top_builddir)/src/Lib/landcover/liblandcover.a \
 	$(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \
 	-lsgbucket -lsgdebug -lsgmisc -lsgxml -lgenpolyclip -lz 

src/Lib/Polygon/polygon.cxx

@@ -36,6 +36,7 @@ extern "C" {
 
 SG_USING_STD(endl);
 
+#include <Geometry/trinodes.hxx>
 #include <poly2tri/interface.h>
 
 #include "polygon.hxx"
@@ -254,6 +255,95 @@ void make_gpc_poly( const TGPolygon& in, gpc_polygon *out ) {
 }
 
 
+// Set the elevations of points in dest based on the elevations of
+// points in source.  For points in dest that are not in source,
+// propogate the value from the nearest matching point.
+static TGPolygon preserve_elevations( const TGPolygon &source,
+                                      const TGPolygon &dest )
+{
+    TGTriNodes nodes;
+    nodes.clear();
+
+    int i, j;
+
+    // build a list of points from the source polygon
+
+    for ( i = 0; i < source.contours(); ++i ) {
+        for ( j = 0; j < source.contour_size(i); ++j ) {
+            Point3D p = source.get_pt( i, j );
+            nodes.unique_add( p );
+        }
+    }
+
+    // traverse the dest polygon and build a mirror image but with
+    // elevations from the source polygon
+
+    TGPolygon result;
+    result.erase();
+
+    for ( i = 0; i < dest.contours(); ++i ) {
+        for ( j = 0; j < dest.contour_size(i); ++j ) {
+            Point3D p = dest.get_pt( i, j );
+            int index = nodes.find( p );
+            if ( index >= 0 ) {
+                Point3D ref = nodes.get_node( index );
+                p.setz( ref.z() );
+            } else {
+                p.setz( -9999.0 );
+            }
+            result.add_node( i, p );
+        }
+    }
+
+    // now post process result to catch any nodes that weren't updated
+    // (because the clipping process may have added points which
+    // weren't in the original.)
+
+    double last = -9999.0;
+    for ( i = 0; i < result.contours(); ++i ) {
+        // go front ways
+        last = -9999.0;
+        for ( j = 0; j < result.contour_size(i); ++j ) {
+            Point3D p = result.get_pt( i, j );
+            if ( p.z() > -9000 ) {
+                last = p.z();
+            } else {
+               if ( last > -9000 ) {
+                   p.setz( last );
+                   result.set_pt( i, j, p );
+               }
+            }
+        }
+
+        // go back ways
+        last = -9999.0;
+        for ( j = result.contour_size(i) - 1; j >= 0; --j ) {
+            Point3D p = result.get_pt( i, j );
+            if ( p.z() > -9000 ) {
+                last = p.z();
+            } else {
+               if ( last > -9000 ) {
+                   p.setz( last );
+                   result.set_pt( i, j, p );
+               }
+            }
+        }
+
+        // finally drop ten and punt on any points that are still
+        // elevation-less and set their elevations to zero.
+        for ( j = 0; j < result.contour_size(i); ++j ) {
+            Point3D p = result.get_pt( i, j );
+            if ( p.z() < -9000 ) {
+                p.setz( 0.0 );
+                result.set_pt( i, j, p );
+            }
+        }
+    }
+
+    return result;
+}
+
+
 // Set operation type
 typedef enum {
     POLY_DIFF,			// Difference
@@ -332,6 +422,8 @@ TGPolygon polygon_clip( clip_op poly_op, const TGPolygon& subject,
     gpc_free_polygon( gpc_clip );
     gpc_free_polygon( gpc_result );
 
+    result = preserve_elevations( clip, result );
+
     return result;
 }
 

src/Lib/Polygon/split-bin.cxx

@@ -40,7 +40,8 @@
 
 
 static void clip_and_write_poly( string root, long int p_index, AreaType area, 
-				 SGBucket b, const TGPolygon& shape ) {
+				 SGBucket b, const TGPolygon& shape,
+                                 bool preserve3d ) {
     Point3D c, min, max, p;
     c = Point3D( b.get_center_lon(), b.get_center_lat(), 0 );
     double span = sg_bucket_span( c.y() );
@@ -112,6 +113,11 @@ static void clip_and_write_poly( string root, long int p_index, AreaType area,
 	    throw sg_exception("unknown area type in clip_and_write_poly()!");
 	
 	FILE *rfp= fopen( polyfile.c_str(), "w" );
+        if ( preserve3d ) {
+            fprintf( rfp, "#3D\n" );
+        } else {
+            fprintf( rfp, "#2D\n" );
+        }
 	fprintf( rfp, "%s\n", poly_type.c_str() );
 
 	fprintf( rfp, "%d\n", result.contours() );
@@ -120,7 +126,11 @@ static void clip_and_write_poly( string root, long int p_index, AreaType area,
 	    fprintf( rfp, "%d\n", result.get_hole_flag(i) );
 	    for ( int j = 0; j < result.contour_size(i); ++j ) {
 		p = result.get_pt( i, j );
-		fprintf( rfp, "%.15f  %.15f\n", p.x(), p.y() );
+                if ( preserve3d ) {
+                    fprintf( rfp, "%.15f  %.15f %.15f\n", p.x(), p.y(), p.z() );
+                } else {
+                    fprintf( rfp, "%.15f  %.15f\n", p.x(), p.y() );
+                }
 	    }
 	}
 	fclose( rfp );
@@ -172,7 +182,7 @@ void tgSplitPolygon( const string& path, AreaType area,
 
     if ( b_min == b_max ) {
 	// shape entirely contained in a single bucket, write and bail
-	clip_and_write_poly( path, index, area, b_min, shape );
+	clip_and_write_poly( path, index, area, b_min, shape, preserve3d );
 	return;
     }
 
@@ -194,7 +204,8 @@ void tgSplitPolygon( const string& path, AreaType area,
 	for ( j = 0; j <= 1; ++j ) {
 	    for ( i = 0; i <= dx; ++i ) {
 	        b_cur = sgBucketOffset(min.x(), min.y(), i, j);
-	        clip_and_write_poly( path, index, area, b_cur, shape );
+	        clip_and_write_poly( path, index, area, b_cur, shape,
+                                     preserve3d );
 	    }
 	}
 	return;

src/Prep/E00Lines/Makefile.am

@@ -3,8 +3,8 @@ bin_PROGRAMS = e00lines
 e00lines_SOURCES = main.cxx
 
 e00lines_LDADD = \
-	$(top_builddir)/src/Lib/Geometry/libGeometry.a \
 	$(top_builddir)/src/Lib/Polygon/libPolygon.a \
+	$(top_builddir)/src/Lib/Geometry/libGeometry.a \
 	$(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \
 	$(top_builddir)/src/Lib/TriangleJRS/libTriangleJRS.a \
 	$(top_builddir)/src/Lib/e00/libe00.a \

src/Prep/GSHHS/Makefile.am

@@ -7,6 +7,7 @@ gshhs_SOURCES = \
 
 gshhs_LDADD = \
 	$(top_builddir)/src/Lib/Polygon/libPolygon.a \
+	$(top_builddir)/src/Lib/Geometry/libGeometry.a \
 	$(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \
 	-lsgdebug -lsgbucket -lsgmisc -lsgxml -lgenpolyclip -lz 
 
@@ -15,6 +16,7 @@ debug_SOURCES = \
 
 debug_LDADD = \
 	$(top_builddir)/src/Lib/Polygon/libPolygon.a \
+	$(top_builddir)/src/Lib/Geometry/libGeometry.a \
 	$(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \
 	-lsgdebug -lsgbucket -lsgmisc -lsgxml -lgenpolyclip -lz
 

src/Prep/MergerClipper/Makefile.am

@@ -9,6 +9,7 @@ testmerger_SOURCES = testmerger.cxx
 testmerger_LDADD = \
 	$(top_builddir)/src/Prep/MergerClipper/libMergerClipper.a \
 	$(top_builddir)/src/Lib/Polygon/libPolygon.a \
+	$(top_builddir)/src/Lib/Geometry/libGeometry.a \
 	$(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \
 	-lsgdebug -lsgbucket -lsgmisc -lsgxml -lgenpolyclip -lz
 

src/Prep/ShapeFile/Makefile.am

@@ -6,12 +6,14 @@ noaa_decode_SOURCES = noaa-decode.cxx
 
 shape_decode_LDADD = \
 	$(top_builddir)/src/Lib/Polygon/libPolygon.a \
+	$(top_builddir)/src/Lib/Geometry/libGeometry.a \
 	$(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \
 	$(top_builddir)/src/Lib/shapelib/libshape.a \
 	-lsgdebug -lsgbucket -lsgmisc -lsgxml -lgenpolyclip -lz
 
 noaa_decode_LDADD = \
 	$(top_builddir)/src/Lib/Polygon/libPolygon.a \
+	$(top_builddir)/src/Lib/Geometry/libGeometry.a \
 	$(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \
 	$(top_builddir)/src/Lib/shapelib/libshape.a \
 	-lsgdebug -lsgbucket -lsgmisc -lsgxml -lgenpolyclip -lz

src/Prep/TGVPF/Makefile.am

@@ -3,8 +3,8 @@ bin_PROGRAMS = tgvpf
 tgvpf_SOURCES = tgvpf.cxx
 
 tgvpf_LDADD = \
-	$(top_builddir)/src/Lib/Geometry/libGeometry.a \
 	$(top_builddir)/src/Lib/Polygon/libPolygon.a \
+	$(top_builddir)/src/Lib/Geometry/libGeometry.a \
 	$(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \
 	$(top_builddir)/src/Lib/TriangleJRS/libTriangleJRS.a \
 	$(top_builddir)/src/Lib/vpf/libvpf.a \

src/Prep/UserDef/Makefile.am

@@ -3,8 +3,8 @@ bin_PROGRAMS = tguserdef
 tguserdef_SOURCES = tguserdef.cxx
 
 tguserdef_LDADD = \
-	$(top_builddir)/src/Lib/Geometry/libGeometry.a \
 	$(top_builddir)/src/Lib/Polygon/libPolygon.a \
+	$(top_builddir)/src/Lib/Geometry/libGeometry.a \
 	$(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \
 	$(top_builddir)/src/Lib/TriangleJRS/libTriangleJRS.a \
 	-lsgbucket -lsgmisc -lsgmath -lsgio -lsgxml -lsgdebug -lgenpolyclip -lz