builtins: implement ST_Boundary, ST_Difference, and ST_Relate with BNR

docs/generated/sql/functions.md

@@ -1377,6 +1377,9 @@ has no relationship with the commit order of concurrent transactions.</p>
 <tr><td><a name="st_azimuth"></a><code>st_azimuth(geometry_a: geometry, geometry_b: geometry) &rarr; <a href="float.html">float</a></code></td><td><span class="funcdesc"><p>Returns the azimuth in radians of the segment defined by the given point geometries, or NULL if the two points are coincident.</p>
 <p>The azimuth is angle is referenced from north, and is positive clockwise: North = 0; East = π/2; South = π; West = 3π/2.</p>
 </span></td></tr>
+<tr><td><a name="st_boundary"></a><code>st_boundary(geometry: geometry) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns the closure of the combinatorial boundary of this Geometry.</p>
+<p>This function utilizes the GEOS module.</p>
+</span></td></tr>
 <tr><td><a name="st_box2dfromgeohash"></a><code>st_box2dfromgeohash(geohash: <a href="string.html">string</a>) &rarr; box2d</code></td><td><span class="funcdesc"><p>Return a Box2D from a GeoHash string with max precision.</p>
 </span></td></tr>
 <tr><td><a name="st_box2dfromgeohash"></a><code>st_box2dfromgeohash(geohash: <a href="string.html">string</a>, precision: <a href="int.html">int</a>) &rarr; box2d</code></td><td><span class="funcdesc"><p>Return a Box2D from a GeoHash string with supplied precision.</p>
@@ -1549,6 +1552,9 @@ from the given Geometry.</p>
 <tr><td><a name="st_dfullywithinexclusive"></a><code>st_dfullywithinexclusive(geometry_a: geometry, geometry_b: geometry, distance: <a href="float.html">float</a>) &rarr; <a href="bool.html">bool</a></code></td><td><span class="funcdesc"><p>Returns true if every pair of points comprising geometry_a and geometry_b are within distance units, exclusive. In other words, the ST_MaxDistance between geometry_a and geometry_b is less than distance units.</p>
 <p>This function variant will attempt to utilize any available spatial index.</p>
 </span></td></tr>
+<tr><td><a name="st_difference"></a><code>st_difference(geometry_a: geometry, geometry_b: geometry) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns the difference of two Geometries.</p>
+<p>This function utilizes the GEOS module.</p>
+</span></td></tr>
 <tr><td><a name="st_dimension"></a><code>st_dimension(geometry: geometry) &rarr; <a href="int.html">int</a></code></td><td><span class="funcdesc"><p>Returns the number of topological dimensions of a given Geometry.</p>
 </span></td></tr>
 <tr><td><a name="st_disjoint"></a><code>st_disjoint(geometry_a: geometry, geometry_b: geometry) &rarr; <a href="bool.html">bool</a></code></td><td><span class="funcdesc"><p>Returns true if geometry_a does not overlap, touch or is within geometry_b.</p>
@@ -1987,6 +1993,9 @@ Negative azimuth values and values greater than 2π (360 degrees) are supported.
 <tr><td><a name="st_relate"></a><code>st_relate(geometry_a: geometry, geometry_b: geometry) &rarr; <a href="string.html">string</a></code></td><td><span class="funcdesc"><p>Returns the DE-9IM spatial relation between geometry_a and geometry_b.</p>
 <p>This function utilizes the GEOS module.</p>
 </span></td></tr>
+<tr><td><a name="st_relate"></a><code>st_relate(geometry_a: geometry, geometry_b: geometry, bnr: <a href="int.html">int</a>) &rarr; <a href="string.html">string</a></code></td><td><span class="funcdesc"><p>Returns the DE-9IM spatial relation between geometry_a and geometry_b using the given boundary node rule (1:OGC/MOD2, 2:Endpoint, 3:MultivalentEndpoint, 4:MonovalentEndpoint).</p>
+<p>This function utilizes the GEOS module.</p>
+</span></td></tr>
 <tr><td><a name="st_relate"></a><code>st_relate(geometry_a: geometry, geometry_b: geometry, pattern: <a href="string.html">string</a>) &rarr; <a href="bool.html">bool</a></code></td><td><span class="funcdesc"><p>Returns whether the DE-9IM spatial relation between geometry_a and geometry_b matches the DE-9IM pattern.</p>
 <p>This function utilizes the GEOS module.</p>
 </span></td></tr>

pkg/geo/geomfn/de9im.go

@@ -25,6 +25,15 @@ func Relate(a geo.Geometry, b geo.Geometry) (string, error) {
  return geos.Relate(a.EWKB(), b.EWKB())
 }
 
+// RelateBoundaryNodeRule returns the DE-9IM relation between A and B using
+// the given boundary node rule (as specified by GEOS).
+func RelateBoundaryNodeRule(a, b geo.Geometry, bnr int) (string, error) {
+ if a.SRID() != b.SRID() {
+ return "", geo.NewMismatchingSRIDsError(a.SpatialObject(), b.SpatialObject())
+ }
+ return geos.RelateBoundaryNodeRule(a.EWKB(), b.EWKB(), bnr)
+}
+
 // RelatePattern returns whether the DE-9IM relation between A and B matches.
 func RelatePattern(a geo.Geometry, b geo.Geometry, pattern string) (bool, error) {
  if a.SRID() != b.SRID() {

pkg/geo/geomfn/de9im_test.go

@@ -36,6 +36,40 @@ func TestRelate(t *testing.T) {
  }
 }
 
+func TestRelateBoundaryNodeRule(t *testing.T) {
+ testCases := []struct {
+ a geo.Geometry
+ b geo.Geometry
+ bnr int
+ expected string
+ }{
+ {leftRect, rightRect, 1, "FF2F11212"},
+ {leftRect, rightRect, 2, "FF2F11212"},
+ {leftRect, rightRect, 3, "1F2F002F2"},
+ {leftRect, rightRect, 4, "FF2F11212"},
+ }
+
+ for i, tc := range testCases {
+ t.Run(fmt.Sprintf("tc:%d", i), func(t *testing.T) {
+ ret, err := RelateBoundaryNodeRule(tc.a, tc.b, tc.bnr)
+ require.NoError(t, err)
+ require.Equal(t, tc.expected, ret)
+ })
+ }
+
+ t.Run("errors on invalid BNR", func(t *testing.T) {
+ _, err := RelateBoundaryNodeRule(leftRect, rightRect, 0)
+ require.Error(t, err)
+ _, err = RelateBoundaryNodeRule(leftRect, rightRect, 5)
+ require.Error(t, err)
+ })
+
+ t.Run("errors if SRIDs mismatch", func(t *testing.T) {
+ _, err := RelateBoundaryNodeRule(mismatchingSRIDGeometryA, mismatchingSRIDGeometryB, 1)
+ requireMismatchingSRIDError(t, err)
+ })
+}
+
 func TestMatchesDE9IM(t *testing.T) {
  testCases := []struct {
  str string

pkg/geo/geomfn/topology_operations.go

@@ -15,6 +15,19 @@ import (
  "github.com/cockroachdb/cockroach/pkg/geo/geos"
 )
 
+// Boundary returns the boundary of a given Geometry.
+func Boundary(g geo.Geometry) (geo.Geometry, error) {
+ // follow PostGIS behavior
+ if g.Empty() {
+ return g, nil
+ }
+ boundaryEWKB, err := geos.Boundary(g.EWKB())
+ if err != nil {
+ return geo.Geometry{}, err
+ }
+ return geo.ParseGeometryFromEWKB(boundaryEWKB)
+}
+
 // Centroid returns the Centroid of a given Geometry.
 func Centroid(g geo.Geometry) (geo.Geometry, error) {
  centroidEWKB, err := geos.Centroid(g.EWKB())
@@ -45,6 +58,22 @@ func ConvexHull(g geo.Geometry) (geo.Geometry, error) {
  return geo.ParseGeometryFromEWKB(convexHullEWKB)
 }
 
+// Difference returns the difference between two given Geometries.
+func Difference(a, b geo.Geometry) (geo.Geometry, error) {
+ // follow PostGIS behavior
+ if a.Empty() || b.Empty() {
+ return a, nil
+ }
+ if a.SRID() != b.SRID() {
+ return geo.Geometry{}, geo.NewMismatchingSRIDsError(a.SpatialObject(), b.SpatialObject())
+ }
+ diffEWKB, err := geos.Difference(a.EWKB(), b.EWKB())
+ if err != nil {
+ return geo.Geometry{}, err
+ }
+ return geo.ParseGeometryFromEWKB(diffEWKB)
+}
+
 // Simplify returns a simplified Geometry.
 func Simplify(g geo.Geometry, tolerance float64) (geo.Geometry, error) {
  simplifiedEWKB, err := geos.Simplify(g.EWKB(), tolerance)

pkg/geo/geomfn/topology_operations_test.go

@@ -20,6 +20,36 @@ import (
  "github.com/twpayne/go-geom"
 )
 
+func TestBoundary(t *testing.T) {
+ testCases := []struct {
+ wkt string
+ expected string
+ }{
+ {"POINT EMPTY", "POINT EMPTY"},
+ {"POINT (1 1)", "GEOMETRYCOLLECTION EMPTY"},
+ {"LINESTRING EMPTY", "LINESTRING EMPTY"},
+ {"LINESTRING (100 150, 50 60, 70 80, 160 170)", "MULTIPOINT (100 150, 160 170)"},
+ {"SRID=4000;LINESTRING (100 150, 50 60, 70 80, 160 170)", "SRID=4000;MULTIPOINT (100 150, 160 170)"},
+ {
+ "POLYGON ((10 130, 50 190, 110 190, 140 150, 150 80, 100 10, 20 40, 10 130), (70 40, 100 50, 120 80, 80 110, 50 90, 70 40))",
+ "MULTILINESTRING ((10 130, 50 190, 110 190, 140 150, 150 80, 100 10, 20 40, 10 130), (70 40, 100 50, 120 80, 80 110, 50 90, 70 40))",
+ },
+ }
+
+ for _, tc := range testCases {
+ t.Run(tc.wkt, func(t *testing.T) {
+ g, err := geo.ParseGeometry(tc.wkt)
+ require.NoError(t, err)
+ ret, err := Boundary(g)
+ require.NoError(t, err)
+
+ wkt, err := geo.SpatialObjectToEWKT(ret.SpatialObject(), 0)
+ require.NoError(t, err)
+ require.EqualValues(t, tc.expected, wkt)
+ })
+ }
+}
+
 func TestCentroid(t *testing.T) {
  testCases := []struct {
  wkt string
@@ -93,6 +123,39 @@ func TestConvexHull(t *testing.T) {
  }
 }
 
+func TestDifference(t *testing.T) {
+ testCases := []struct {
+ wkt1 string
+ wkt2 string
+ expected string
+ }{
+ {"POINT EMPTY", "LINESTRING EMPTY", "POINT EMPTY"},
+ {"LINESTRING EMPTY", "POINT EMPTY", "LINESTRING EMPTY"},
+ {"LINESTRING (50 100, 50 200)", "LINESTRING(50 50, 50 150)", "LINESTRING (50 150, 50 200)"},
+ {"SRID=4000;LINESTRING (50 100, 50 200)", "SRID=4000;LINESTRING(50 50, 50 150)", "SRID=4000;LINESTRING (50 150, 50 200)"},
+ }
+
+ for _, tc := range testCases {
+ t.Run(fmt.Sprintf("%v - %v", tc.wkt2, tc.wkt1), func(t *testing.T) {
+ g1, err := geo.ParseGeometry(tc.wkt1)
+ require.NoError(t, err)
+ g2, err := geo.ParseGeometry(tc.wkt2)
+ require.NoError(t, err)
+ ret, err := Difference(g1, g2)
+ require.NoError(t, err)
+
+ wkt, err := geo.SpatialObjectToEWKT(ret.SpatialObject(), 0)
+ require.NoError(t, err)
+ require.EqualValues(t, tc.expected, wkt)
+ })
+ }
+
+ t.Run("errors if SRIDs mismatch", func(t *testing.T) {
+ _, err := Difference(mismatchingSRIDGeometryA, mismatchingSRIDGeometryB)
+ requireMismatchingSRIDError(t, err)
+ })
+}
+
 func TestSimplify(t *testing.T) {
  testCases := []struct {
  wkt string

pkg/geo/geos/geos.cc

@@ -97,8 +97,11 @@ typedef int (*CR_GEOS_Length_r)(CR_GEOS_Handle, CR_GEOS_Geometry, double*);
 typedef int (*CR_GEOS_Normalize_r)(CR_GEOS_Handle, CR_GEOS_Geometry);
 typedef CR_GEOS_Geometry (*CR_GEOS_LineMerge_r)(CR_GEOS_Handle, CR_GEOS_Geometry);
 
+typedef CR_GEOS_Geometry (*CR_GEOS_Boundary_r)(CR_GEOS_Handle, CR_GEOS_Geometry);
 typedef CR_GEOS_Geometry (*CR_GEOS_Centroid_r)(CR_GEOS_Handle, CR_GEOS_Geometry);
 typedef CR_GEOS_Geometry (*CR_GEOS_ConvexHull_r)(CR_GEOS_Handle, CR_GEOS_Geometry);
+typedef CR_GEOS_Geometry (*CR_GEOS_Difference_r)(CR_GEOS_Handle, CR_GEOS_Geometry,
+ CR_GEOS_Geometry);
 typedef CR_GEOS_Geometry (*CR_GEOS_Simplify_r)(CR_GEOS_Handle, CR_GEOS_Geometry, double);
 typedef CR_GEOS_Geometry (*CR_GEOS_TopologyPreserveSimplify_r)(CR_GEOS_Handle, CR_GEOS_Geometry,
  double);
@@ -125,6 +128,8 @@ typedef char (*CR_GEOS_Touches_r)(CR_GEOS_Handle, CR_GEOS_Geometry, CR_GEOS_Geom
 typedef char (*CR_GEOS_Within_r)(CR_GEOS_Handle, CR_GEOS_Geometry, CR_GEOS_Geometry);
 
 typedef char* (*CR_GEOS_Relate_r)(CR_GEOS_Handle, CR_GEOS_Geometry, CR_GEOS_Geometry);
+typedef char* (*CR_GEOS_RelateBoundaryNodeRule_r)(CR_GEOS_Handle, CR_GEOS_Geometry,
+ CR_GEOS_Geometry, int);
 typedef char (*CR_GEOS_RelatePattern_r)(CR_GEOS_Handle, CR_GEOS_Geometry, CR_GEOS_Geometry,
  const char*);
 
@@ -190,8 +195,10 @@ struct CR_GEOS {
  CR_GEOS_Normalize_r GEOSNormalize_r;
  CR_GEOS_LineMerge_r GEOSLineMerge_r;
 
+ CR_GEOS_Boundary_r GEOSBoundary_r;
  CR_GEOS_Centroid_r GEOSGetCentroid_r;
  CR_GEOS_ConvexHull_r GEOSConvexHull_r;
+ CR_GEOS_Difference_r GEOSDifference_r;
  CR_GEOS_Simplify_r GEOSSimplify_r;
  CR_GEOS_TopologyPreserveSimplify_r GEOSTopologyPreserveSimplify_r;
  CR_GEOS_Union_r GEOSUnion_r;
@@ -215,6 +222,7 @@ struct CR_GEOS {
  CR_GEOS_Within_r GEOSWithin_r;
 
  CR_GEOS_Relate_r GEOSRelate_r;
+ CR_GEOS_RelateBoundaryNodeRule_r GEOSRelateBoundaryNodeRule_r;
  CR_GEOS_RelatePattern_r GEOSRelatePattern_r;
 
  CR_GEOS_WKBWriter_create_r GEOSWKBWriter_create_r;
@@ -276,6 +284,8 @@ struct CR_GEOS {
  INIT(GEOSNormalize_r);
  INIT(GEOSLineMerge_r);
  INIT(GEOSisSimple_r);
+ INIT(GEOSBoundary_r);
+ INIT(GEOSDifference_r);
  INIT(GEOSGetCentroid_r);
  INIT(GEOSConvexHull_r);
  INIT(GEOSSimplify_r);
@@ -297,6 +307,7 @@ struct CR_GEOS {
  INIT(GEOSTouches_r);
  INIT(GEOSWithin_r);
  INIT(GEOSRelate_r);
+ INIT(GEOSRelateBoundaryNodeRule_r);
  INIT(GEOSRelatePattern_r);
  INIT(GEOSSharedPaths_r);
  INIT(GEOSWKTReader_create_r);
@@ -705,6 +716,24 @@ CR_GEOS_Status CR_GEOS_MakeValid(CR_GEOS* lib, CR_GEOS_Slice g, CR_GEOS_String*
 // Topology operators.
 //
 
+CR_GEOS_Status CR_GEOS_Boundary(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_String* boundaryEWKB) {
+ std::string error;
+ auto handle = initHandleWithErrorBuffer(lib, &error);
+ auto geom = CR_GEOS_GeometryFromSlice(lib, handle, a);
+ *boundaryEWKB = {.data = NULL, .len = 0};
+ if (geom != nullptr) {
+ auto boundaryGeom = lib->GEOSBoundary_r(handle, geom);
+ if (boundaryGeom != nullptr) {
+ auto srid = lib->GEOSGetSRID_r(handle, geom);
+ CR_GEOS_writeGeomToEWKB(lib, handle, boundaryGeom, boundaryEWKB, srid);
+ lib->GEOSGeom_destroy_r(handle, boundaryGeom);
+ }
+ lib->GEOSGeom_destroy_r(handle, geom);
+ }
+ lib->GEOS_finish_r(handle);
+ return toGEOSString(error.data(), error.length());
+}
+
 CR_GEOS_Status CR_GEOS_Centroid(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_String* centroidEWKB) {
  std::string error;
  auto handle = initHandleWithErrorBuffer(lib, &error);
@@ -741,6 +770,31 @@ CR_GEOS_Status CR_GEOS_ConvexHull(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_String*
  return toGEOSString(error.data(), error.length());
 }
 
+CR_GEOS_Status CR_GEOS_Difference(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_Slice b,
+ CR_GEOS_String* diffEWKB) {
+ std::string error;
+ auto handle = initHandleWithErrorBuffer(lib, &error);
+ auto geomA = CR_GEOS_GeometryFromSlice(lib, handle, a);
+ auto geomB = CR_GEOS_GeometryFromSlice(lib, handle, b);
+ *diffEWKB = {.data = NULL, .len = 0};
+ if (geomA != nullptr && geomB != nullptr) {
+ auto diffGeom = lib->GEOSDifference_r(handle, geomA, geomB);
+ if (diffGeom != nullptr) {
+ auto srid = lib->GEOSGetSRID_r(handle, geomA);
+ CR_GEOS_writeGeomToEWKB(lib, handle, diffGeom, diffEWKB, srid);
+ lib->GEOSGeom_destroy_r(handle, diffGeom);
+ }
+ }
+ if (geomA != nullptr) {
+ lib->GEOSGeom_destroy_r(handle, geomA);
+ }
+ if (geomB != nullptr) {
+ lib->GEOSGeom_destroy_r(handle, geomB);
+ }
+ lib->GEOS_finish_r(handle);
+ return toGEOSString(error.data(), error.length());
+}
+
 CR_GEOS_Status CR_GEOS_Simplify(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_String* simplifyEWKB,
  double tolerance) {
  std::string error;
@@ -1017,6 +1071,33 @@ CR_GEOS_Status CR_GEOS_Relate(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_Slice b, CR
  return toGEOSString(error.data(), error.length());
 }
 
+CR_GEOS_Status CR_GEOS_RelateBoundaryNodeRule(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_Slice b,
+ int bnr, CR_GEOS_String* ret) {
+ std::string error;
+ auto handle = initHandleWithErrorBuffer(lib, &error);
+
+ auto wkbReader = lib->GEOSWKBReader_create_r(handle);
+ auto geomA = lib->GEOSWKBReader_read_r(handle, wkbReader, a.data, a.len);
+ auto geomB = lib->GEOSWKBReader_read_r(handle, wkbReader, b.data, b.len);
+ lib->GEOSWKBReader_destroy_r(handle, wkbReader);
+
+ if (geomA != nullptr && geomB != nullptr) {
+ auto r = lib->GEOSRelateBoundaryNodeRule_r(handle, geomA, geomB, bnr);
+ if (r != NULL) {
+ *ret = toGEOSString(r, strlen(r));
+ lib->GEOSFree_r(handle, r);
+ }
+ }
+ if (geomA != nullptr) {
+ lib->GEOSGeom_destroy_r(handle, geomA);
+ }
+ if (geomB != nullptr) {
+ lib->GEOSGeom_destroy_r(handle, geomB);
+ }
+ lib->GEOS_finish_r(handle);
+ return toGEOSString(error.data(), error.length());
+}
+
 CR_GEOS_Status CR_GEOS_RelatePattern(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_Slice b,
  CR_GEOS_Slice pattern, char* ret) {
  std::string error;