GEOSSnap: Preserve Z values from input

capi/geos_c.h.in

@@ -5505,6 +5505,9 @@ extern GEOSGeometry GEOS_DLL *GEOSGeom_transformXYZ(
 * Snaps the vertices and segments of the first geometry to vertices of the
 * second geometry within the given tolerance.
 *
+* Z values from the input will be preserved.
+* M values from the input will be discarded.
+*
 * Where possible, this operation tries to avoid creating invalid geometries;
 * however, it does not guarantee that output geometries will be valid.  It is
 * the responsibility of the caller to check for and handle invalid geometries.

include/geos/operation/overlay/snap/GeometrySnapper.h

@@ -138,7 +138,7 @@ class GEOS_DLL GeometrySnapper {
     const geom::Geometry& srcGeom;
 
     /// Extract target (unique) coordinates
-    std::unique_ptr<geom::Coordinate::ConstVect> extractTargetCoordinates(
+    static std::unique_ptr<geom::Coordinate::ConstVect> extractTargetCoordinates(
         const geom::Geometry& g);
 
     // Declare type as noncopyable

src/operation/overlay/snap/GeometrySnapper.cpp

@@ -51,6 +51,8 @@ class SnapTransformer: public geos::geom::util::GeometryTransformer {
     double snapTol;
 
     const Coordinate::ConstVect& snapPts;
+    bool constructZ;
+    bool constructM;
 
     CoordinateSequence::Ptr
     snapLine(
@@ -59,7 +61,7 @@ class SnapTransformer: public geos::geom::util::GeometryTransformer {
         using std::unique_ptr;
 
         assert(srcPts);
-        auto coords = detail::make_unique<CoordinateSequence>();
+        auto coords = detail::make_unique<CoordinateSequence>(0, constructZ, constructM);
         coords->add(*srcPts);
         LineStringSnapper snapper(*coords, snapTol);
         return snapper.snapTo(snapPts);
@@ -68,10 +70,13 @@ class SnapTransformer: public geos::geom::util::GeometryTransformer {
 public:
 
     SnapTransformer(double nSnapTol,
-                    const Coordinate::ConstVect& nSnapPts)
+                    const Coordinate::ConstVect& nSnapPts,
+                    bool p_constructZ, bool p_constructM)
         :
         snapTol(nSnapTol),
-        snapPts(nSnapPts)
+        snapPts(nSnapPts),
+        constructZ(p_constructZ),
+        constructM(p_constructM)
     {
     }
 
@@ -112,7 +117,7 @@ GeometrySnapper::snapTo(const geom::Geometry& g, double snapTolerance)
 
     // Apply a SnapTransformer to source geometry
     // (we need a pointer for dynamic polymorphism)
-    std::unique_ptr<GeometryTransformer> snapTrans(new SnapTransformer(snapTolerance, *snapPts));
+    std::unique_ptr<GeometryTransformer> snapTrans(new SnapTransformer(snapTolerance, *snapPts, srcGeom.hasZ(), false));
     return snapTrans->transform(&srcGeom);
 }
 
@@ -129,7 +134,7 @@ GeometrySnapper::snapToSelf(double snapTolerance, bool cleanResult)
 
     // Apply a SnapTransformer to source geometry
     // (we need a pointer for dynamic polymorphism)
-    std::unique_ptr<GeometryTransformer> snapTrans(new SnapTransformer(snapTolerance, *snapPts));
+    std::unique_ptr<GeometryTransformer> snapTrans(new SnapTransformer(snapTolerance, *snapPts, srcGeom.hasZ(), false));
 
     std::unique_ptr<Geometry> result = snapTrans->transform(&srcGeom);
 

src/operation/overlay/snap/LineStringSnapper.cpp

@@ -153,7 +153,8 @@ LineStringSnapper::snapVertices(geom::CoordinateList& srcCoords,
 #if GEOS_DEBUG
         std::cerr << " Vertex to be snapped found, snapping" << std::endl;
 #endif
-        *vertpos = snapPt;
+        vertpos->x = snapPt.x;
+        vertpos->y = snapPt.y;
 
         // keep final closing point in synch (rings only)
         if(vertpos == srcCoords.begin() && isClosed) {
@@ -162,7 +163,8 @@ LineStringSnapper::snapVertices(geom::CoordinateList& srcCoords,
 #if GEOS_DEBUG
             std::cerr << " Snapped vertex was first in a closed line, also snapping last" << std::endl;
 #endif
-            *vertpos = snapPt;
+            vertpos->x = snapPt.x;
+            vertpos->y = snapPt.y;
         }
 
 #if GEOS_DEBUG

tests/unit/capi/GEOSSnapTest.cpp

@@ -220,5 +220,26 @@ void object::test<11>()
     ensure("curved geometry not supported", result_ == nullptr);
 }
 
+template<>
+template<>
+void object::test<12>
+()
+{
+    set_test_name("POLYGON ZM snapped to Point");
+
+    geom1_ = fromWKT("POLYGON ZM ((0 0 4 3, 10 0 2 8, 10 10 6 1, 0 1 -3 4, 0 0 4 3))");
+    geom2_ = fromWKT("POINT (0.5 0)");
+    ensure(geom1_);
+    ensure(geom2_);
+
+    result_ = GEOSSnap(geom1_, geom2_, 1);
+    ensure(result_);
+
+    expected_ = fromWKT("POLYGON Z ((0.5 0 4, 10 0 2, 10 10 6, 0 1 -3, 0.5 0 4))");
+    ensure(expected_);
+
+    ensure_geometry_equals_identical(result_, expected_);
+}
+
 } // namespace tut