geo/geomfn: implement ST_Snap

docs/generated/sql/functions.md

@@ -2217,6 +2217,10 @@ The paths themselves are given in the direction of the first geometry.</p>
 </span></td></tr>
 <tr><td><a name="st_simplifypreservetopology"></a><code>st_simplifypreservetopology(geometry: geometry, tolerance: <a href="float.html">float</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Simplifies the given geometry using the Douglas-Peucker algorithm, avoiding the creation of invalid geometries.</p>
 </span></td></tr>
+<tr><td><a name="st_snap"></a><code>st_snap(input: geometry, target: geometry, tolerance: <a href="float.html">float</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Snaps the vertices and segments of input geometry the target geometry’s vertices.
+Tolerance is used to control where snapping is performed. The result geometry is the input geometry with the vertices snapped.
+If no snapping occurs then the input geometry is returned unchanged.</p>
+</span></td></tr>
 <tr><td><a name="st_snaptogrid"></a><code>st_snaptogrid(geometry: geometry, origin_x: <a href="float.html">float</a>, origin_y: <a href="float.html">float</a>, size_x: <a href="float.html">float</a>, size_y: <a href="float.html">float</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Snap a geometry to a grid of with X coordinates snapped to size_x and Y coordinates snapped to size_y based on an origin of (origin_x, origin_y).</p>
 </span></td></tr>
 <tr><td><a name="st_snaptogrid"></a><code>st_snaptogrid(geometry: geometry, size: <a href="float.html">float</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Snap a geometry to a grid of the given size.</p>

pkg/geo/geomfn/BUILD.bazel

@@ -26,6 +26,7 @@ go_library(
         "reverse.go",
         "segmentize.go",
         "shift_longitude.go",
+        "snap.go",
         "snap_to_grid.go",
         "subdivide.go",
         "swap_ordinates.go",
@@ -80,6 +81,7 @@ go_test(
         "reverse_test.go",
         "segmentize_test.go",
         "shift_longitude_test.go",
+        "snap_test.go",
         "snap_to_grid_test.go",
         "subdivide_test.go",
         "swap_ordinates_test.go",

pkg/geo/geomfn/snap.go

@@ -0,0 +1,27 @@
+// Copyright 2021 The Cockroach Authors.
+//
+// Use of this software is governed by the Business Source License
+// included in the file licenses/BSL.txt.
+//
+// As of the Change Date specified in that file, in accordance with
+// the Business Source License, use of this software will be governed
+// by the Apache License, Version 2.0, included in the file
+// licenses/APL.txt.
+
+package geomfn
+
+import (
+	"github.com/cockroachdb/cockroach/pkg/geo"
+	"github.com/cockroachdb/cockroach/pkg/geo/geos"
+)
+
+// Snap returns the input geometry with the vertices snapped to the target
+// geometry. Tolerance is used to control where snapping is performed.
+// If no snapping occurs then the input geometry is returned unchanged.
+func Snap(input, target geo.Geometry, tolerance float64) (geo.Geometry, error) {
+	snappedEWKB, err := geos.Snap(input.EWKB(), target.EWKB(), tolerance)
+	if err != nil {
+		return geo.Geometry{}, err
+	}
+	return geo.ParseGeometryFromEWKB(snappedEWKB)
+}

pkg/geo/geomfn/snap_test.go

@@ -0,0 +1,87 @@
+// Copyright 2021 The Cockroach Authors.
+//
+// Use of this software is governed by the Business Source License
+// included in the file licenses/BSL.txt.
+//
+// As of the Change Date specified in that file, in accordance with
+// the Business Source License, use of this software will be governed
+// by the Apache License, Version 2.0, included in the file
+// licenses/APL.txt.
+
+package geomfn
+
+import (
+	"testing"
+
+	"github.com/cockroachdb/cockroach/pkg/geo"
+	"github.com/stretchr/testify/require"
+	"github.com/twpayne/go-geom"
+)
+
+func TestSnap(t *testing.T) {
+	testCases := []struct {
+		desc      string
+		input     geom.T
+		target    geom.T
+		tolerance float64
+		expected  geom.T
+	}{
+		{
+			desc:      "Test snap on two linestrings with tolerance of 0",
+			input:     geom.NewLineStringFlat(geom.XY, []float64{20, 25, 30, 35}),
+			target:    geom.NewLineStringFlat(geom.XY, []float64{10, 15, 20, 25}),
+			tolerance: 0,
+			expected:  geom.NewLineStringFlat(geom.XY, []float64{20, 25, 30, 35}),
+		},
+		{
+			desc:      "Test snapping a polygon on a linestring with tolerance of 150",
+			input:     geom.NewLineStringFlat(geom.XY, []float64{10, 55, 78, 84, 100, 200}),
+			target:    geom.NewPolygonFlat(geom.XY, []float64{26, 125, 26, 200, 126, 200, 126, 125, 26, 125}, []int{10}),
+			tolerance: 150,
+			expected:  geom.NewLineStringFlat(geom.XY, []float64{10, 55, 26, 125, 26, 200, 126, 200, 126, 125}),
+		},
+		{
+			desc:      "Test snapping a linestring on a polygon with tolerance of 150",
+			target:    geom.NewLineStringFlat(geom.XY, []float64{10, 55, 78, 84, 100, 200}),
+			input:     geom.NewPolygonFlat(geom.XY, []float64{26, 125, 26, 200, 126, 200, 126, 125, 26, 125}, []int{10}),
+			tolerance: 150,
+			expected:  geom.NewPolygonFlat(geom.XY, []float64{10, 55, 26, 200, 100, 200, 78, 84, 10, 55}, []int{10}),
+		},
+		{
+			desc:      "Test snapping a linestring on a multipolygon with tolerance of 200",
+			input:     geom.NewLineStringFlat(geom.XY, []float64{5, 107, 54, 84, 101, 100}),
+			target:    geom.NewMultiPolygonFlat(geom.XY, []float64{1, 1, 2, 2, 3, 3, 1, 1, 4, 4, 5, 5, 6, 6, 4, 4}, [][]int{{8}, {16}}),
+			tolerance: 200,
+			expected:  geom.NewLineStringFlat(geom.XY, []float64{5, 107, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 101, 100}),
+		},
+		{
+			desc:      "Test snapping a multipolygon on a linestring with tolerance of 200",
+			input:     geom.NewMultiPolygonFlat(geom.XY, []float64{1, 1, 2, 2, 3, 3, 1, 1, 4, 4, 5, 5, 6, 6, 4, 4}, [][]int{{8}, {16}}),
+			target:    geom.NewLineStringFlat(geom.XY, []float64{5, 107, 54, 84, 101, 100}),
+			tolerance: 200,
+			expected:  geom.NewMultiPolygonFlat(geom.XY, []float64{1, 1, 2, 2, 5, 107, 54, 84, 101, 100, 1, 1, 4, 4, 5, 5, 5, 107, 54, 84, 101, 100, 4, 4}, [][]int{{8}, {16}}),
+		},
+	}
+
+	for _, tc := range testCases {
+		t.Run(tc.desc, func(t *testing.T) {
+			input, err := geo.MakeGeometryFromGeomT(tc.input)
+			require.NoError(t, err)
+			target, err := geo.MakeGeometryFromGeomT(tc.target)
+			require.NoError(t, err)
+
+			actual, err := Snap(input, target, tc.tolerance)
+			require.NoError(t, err)
+
+			// Compare FlatCoords and assert they are within epsilon.
+			// This is because they exact matches may encounter rounding issues.
+			actualGeomT, err := actual.AsGeomT()
+
+			require.NoError(t, err)
+			require.Equal(t, tc.expected.SRID(), actualGeomT.SRID())
+			require.Equal(t, tc.expected.Layout(), actualGeomT.Layout())
+			require.IsType(t, tc.expected, actualGeomT)
+			require.InEpsilonSlice(t, tc.expected.FlatCoords(), actualGeomT.FlatCoords(), 0.00001)
+		})
+	}
+}

pkg/geo/geos/geos.cc

@@ -174,6 +174,8 @@ typedef char (*CR_GEOS_EqualsExact_r)(CR_GEOS_Handle, CR_GEOS_Geometry,
 
 typedef CR_GEOS_Geometry (*CR_GEOS_MinimumRotatedRectangle_r)(CR_GEOS_Handle, CR_GEOS_Geometry);
 
+typedef CR_GEOS_Geometry (*CR_GEOS_Snap_r)(CR_GEOS_Handle, CR_GEOS_Geometry, CR_GEOS_Geometry, double);
+
 std::string ToString(CR_GEOS_Slice slice) { return std::string(slice.data, slice.len); }
 
 }  // namespace
@@ -283,6 +285,8 @@ struct CR_GEOS {
 
   CR_GEOS_Node_r GEOSNode_r;
 
+  CR_GEOS_Snap_r GEOSSnap_r;
+
   CR_GEOS(dlhandle geoscHandle, dlhandle geosHandle)
       : geoscHandle(geoscHandle), geosHandle(geosHandle) {}
 
@@ -385,6 +389,7 @@ struct CR_GEOS {
     INIT(GEOSWKBWriter_write_r);
     INIT(GEOSClipByRect_r);
     INIT(GEOSNode_r);
+    INIT(GEOSSnap_r);
     return nullptr;
 
 #undef INIT
@@ -1523,3 +1528,27 @@ CR_GEOS_Status CR_GEOS_MinimumRotatedRectangle(CR_GEOS* lib, CR_GEOS_Slice g, CR
   lib->GEOS_finish_r(handle);
   return toGEOSString(error.data(), error.length());
 }
+
+CR_GEOS_Status CR_GEOS_Snap(CR_GEOS* lib, CR_GEOS_Slice input, CR_GEOS_Slice target, double tolerance, CR_GEOS_String* snappedEWKB) {
+  std::string error;
+  auto handle = initHandleWithErrorBuffer(lib, &error);
+  auto gGeomInput = CR_GEOS_GeometryFromSlice(lib, handle, input);
+  auto gGeomTarget = CR_GEOS_GeometryFromSlice(lib, handle, target);
+  *snappedEWKB = {.data = NULL, .len = 0};
+  if (gGeomInput != nullptr && gGeomTarget != nullptr) {
+    auto r = lib->GEOSSnap_r(handle, gGeomInput, gGeomTarget, tolerance);
+    if (r != NULL) {
+      auto srid = lib->GEOSGetSRID_r(handle, r);
+      CR_GEOS_writeGeomToEWKB(lib, handle, r, snappedEWKB, srid);
+      lib->GEOSGeom_destroy_r(handle, r);
+    }
+  }
+  if (gGeomInput != nullptr) {
+    lib->GEOSGeom_destroy_r(handle, gGeomInput);
+  }
+  if (gGeomTarget != nullptr) {
+    lib->GEOSGeom_destroy_r(handle, gGeomTarget);
+  }
+  lib->GEOS_finish_r(handle);
+  return toGEOSString(error.data(), error.length());
+}

pkg/geo/geos/geos.go

@@ -1073,3 +1073,20 @@ func MinimumRotatedRectangle(ewkb geopb.EWKB) (geopb.EWKB, error) {
 	}
 	return cStringToSafeGoBytes(cEWKB), nil
 }
+
+// Snap returns the input EWKB with the vertices snapped to the target
+// EWKB. Tolerance is used to control where snapping is performed.
+// If no snapping occurs then the input geometry is returned unchanged.
+func Snap(input, target geopb.EWKB, tolerance float64) (geopb.EWKB, error) {
+	g, err := ensureInitInternal()
+	if err != nil {
+		return nil, err
+	}
+	var cEWKB C.CR_GEOS_String
+	if err := statusToError(
+		C.CR_GEOS_Snap(g, goToCSlice(input), goToCSlice(target), C.double(tolerance), &cEWKB),
+	); err != nil {
+		return nil, err
+	}
+	return cStringToSafeGoBytes(cEWKB), nil
+}

pkg/geo/geos/geos.h

@@ -192,6 +192,8 @@ CR_GEOS_Status CR_GEOS_RelatePattern(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_Slic
 CR_GEOS_Status CR_GEOS_VoronoiDiagram(CR_GEOS* lib, CR_GEOS_Slice g, CR_GEOS_Slice env,
                                       double tolerance, int onlyEdges, CR_GEOS_String* ret);
 
+CR_GEOS_Status CR_GEOS_Snap(CR_GEOS* lib, CR_GEOS_Slice input, CR_GEOS_Slice target, double tolerance, CR_GEOS_String* ret);
+
 #ifdef __cplusplus
 }  // extern "C"
 #endif

pkg/sql/logictest/testdata/logic_test/geospatial

@@ -5768,3 +5768,17 @@ SRID=4326;POINT (359 90)
 SRID=4326;POINT (-179 -23)
 SRID=4326;LINESTRING (2 2, 5 -60, -160 0)
 SRID=4326;POLYGON ((354 -10, 6 -10, 0 20, 354 -10), (3 2, 359 -1, 1 -5, 3 2))
+
+query T
+SELECT ST_AsText(
+    ST_Snap(poly,line, ST_Distance(poly, line)*1.25)
+) AS polysnapped
+FROM (SELECT
+    ST_GeomFromText('MULTIPOLYGON(
+      (( 26 125, 26 200, 126 200, 126 125, 26 125 ),
+      ( 51 150, 101 150, 76 175, 51 150 )),
+      (( 151 100, 151 200, 176 175, 151 100 )))') As poly,
+    ST_GeomFromText('LINESTRING (5 107, 54 84, 101 100)') As line
+	) tbl(poly, line);
+----
+MULTIPOLYGON (((5 107, 26 200, 126 200, 126 125, 101 100, 54 84, 5 107), (51 150, 101 150, 76 175, 51 150)), ((151 100, 151 200, 176 175, 151 100)))

pkg/sql/sem/builtins/geo_builtins.go

@@ -5011,6 +5011,33 @@ The calculations are done on a sphere.`,
 			Volatility: tree.VolatilityImmutable,
 		},
 	),
+	"st_snap": makeBuiltin(
+		defProps(),
+		tree.Overload{
+			Types: tree.ArgTypes{
+				{"input", types.Geometry},
+				{"target", types.Geometry},
+				{"tolerance", types.Float},
+			},
+			ReturnType: tree.FixedReturnType(types.Geometry),
+			Fn: func(_ *tree.EvalContext, args tree.Datums) (tree.Datum, error) {
+				g1 := tree.MustBeDGeometry(args[0])
+				g2 := tree.MustBeDGeometry(args[1])
+				tolerance := tree.MustBeDFloat(args[2])
+				ret, err := geomfn.Snap(g1.Geometry, g2.Geometry, float64(tolerance))
+				if err != nil {
+					return nil, err
+				}
+				return tree.NewDGeometry(ret), nil
+			},
+			Info: infoBuilder{
+				info: `Snaps the vertices and segments of input geometry the target geometry's vertices.
+Tolerance is used to control where snapping is performed. The result geometry is the input geometry with the vertices snapped. 
+If no snapping occurs then the input geometry is returned unchanged.`,
+			}.String(),
+			Volatility: tree.VolatilityImmutable,
+		},
+	),
 	"st_buffer": makeBuiltin(
 		defProps(),
 		tree.Overload{
@@ -6404,7 +6431,6 @@ May return a Point or LineString in the case of degenerate inputs.`,
 	"st_quantizecoordinates":   makeBuiltin(tree.FunctionProperties{UnsupportedWithIssue: 49012}),
 	"st_seteffectivearea":      makeBuiltin(tree.FunctionProperties{UnsupportedWithIssue: 49030}),
 	"st_simplifyvw":            makeBuiltin(tree.FunctionProperties{UnsupportedWithIssue: 49039}),
-	"st_snap":                  makeBuiltin(tree.FunctionProperties{UnsupportedWithIssue: 49040}),
 	"st_split":                 makeBuiltin(tree.FunctionProperties{UnsupportedWithIssue: 49045}),
 	"st_tileenvelope":          makeBuiltin(tree.FunctionProperties{UnsupportedWithIssue: 49053}),
 	"st_wrapx":                 makeBuiltin(tree.FunctionProperties{UnsupportedWithIssue: 49068}),