geo/geomfn: implement ST_SymDifference

Added ST_SymDifference builtin.

Release justification: low-risk update to new functionality
Release note (sql change): Implemented the geometry based builtin
`ST_SymDifference`.

Resolves: #49051

docs/generated/sql/functions.md

@@ -1638,6 +1638,9 @@ The paths themselves are given in the direction of the first geometry.</p>
 <li>S: has spatial reference system</li>
 </ul>
 </span></td></tr>
+<tr><td><a name="st_symdifference"></a><code>st_symdifference(geometry_a: geometry, geometry_b: geometry) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns the symmetric difference of both geometries.</p>
+<p>This function utilizes the GEOS module.</p>
+</span></td></tr>
 <tr><td><a name="st_touches"></a><code>st_touches(geometry_a: geometry, geometry_b: geometry) &rarr; <a href="bool.html">bool</a></code></td><td><span class="funcdesc"><p>Returns true if the only points in common between geometry_a and geometry_b are on the boundary. Note points do not touch other points.</p>
 <p>This function utilizes the GEOS module.</p>
 <p>This function variant will attempt to utilize any available geospatial index.</p>

pkg/geo/geomfn/topology_operations.go

@@ -66,7 +66,7 @@ func Intersection(a geo.Geometry, b geo.Geometry) (geo.Geometry, error) {
 	return geo.ParseGeometryFromEWKB(retEWKB)
 }
 
-// Union returns the geometries of intersection between A and B.
+// Union returns the geometries of union between A and B.
 func Union(a geo.Geometry, b geo.Geometry) (geo.Geometry, error) {
 	if a.SRID() != b.SRID() {
 		return geo.Geometry{}, geo.NewMismatchingSRIDsError(a.SpatialObject(), b.SpatialObject())
@@ -78,6 +78,18 @@ func Union(a geo.Geometry, b geo.Geometry) (geo.Geometry, error) {
 	return geo.ParseGeometryFromEWKB(retEWKB)
 }
 
+// SymDifference returns the geometries of symmetric difference between A and B.
+func SymDifference(a geo.Geometry, b geo.Geometry) (geo.Geometry, error) {
+	if a.SRID() != b.SRID() {
+		return geo.Geometry{}, geo.NewMismatchingSRIDsError(a.SpatialObject(), b.SpatialObject())
+	}
+	retEWKB, err := geos.SymDifference(a.EWKB(), b.EWKB())
+	if err != nil {
+		return geo.Geometry{}, err
+	}
+	return geo.ParseGeometryFromEWKB(retEWKB)
+}
+
 // SharedPaths Returns a geometry collection containing paths shared by the two input geometries.
 func SharedPaths(a geo.Geometry, b geo.Geometry) (geo.Geometry, error) {
 	if a.SRID() != b.SRID() {

pkg/geo/geomfn/topology_operations_test.go

@@ -183,6 +183,32 @@ func TestUnion(t *testing.T) {
 	})
 }
 
+func TestSymDifference(t *testing.T) {
+	testCases := []struct {
+		a        geo.Geometry
+		b        geo.Geometry
+		expected geo.Geometry
+	}{
+		{rightRect, rightRect, emptyRect},
+		{leftRect, rightRect, geo.MustParseGeometry("POLYGON((0 0, -1 0, -1 1, 0 1, 1 1, 1 0, 0 0))")},
+		{leftRect, overlappingRightRect, geo.MustParseGeometry("MULTIPOLYGON(((-0.1 0, -1 0, -1 1, -0.1 1, -0.1 0)), ((0 0, 0 1, 1 1, 1 0, 0 0)))")},
+		{rightRect, rightRectPoint, geo.MustParseGeometry("POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))")},
+	}
+
+	for i, tc := range testCases {
+		t.Run(fmt.Sprintf("tc:%d", i), func(t *testing.T) {
+			g, err := SymDifference(tc.a, tc.b)
+			require.NoError(t, err)
+			require.Equal(t, tc.expected, g)
+		})
+	}
+
+	t.Run("errors if SRIDs mismatch", func(t *testing.T) {
+		_, err := SymDifference(mismatchingSRIDGeometryA, mismatchingSRIDGeometryB)
+		requireMismatchingSRIDError(t, err)
+	})
+}
+
 func TestSharedPaths(t *testing.T) {
 	type args struct {
 		a geo.Geometry

pkg/geo/geos/geos.cc

@@ -100,6 +100,8 @@ typedef CR_GEOS_Geometry (*CR_GEOS_ConvexHull_r)(CR_GEOS_Handle, CR_GEOS_Geometr
 typedef CR_GEOS_Geometry (*CR_GEOS_Union_r)(CR_GEOS_Handle, CR_GEOS_Geometry, CR_GEOS_Geometry);
 typedef CR_GEOS_Geometry (*CR_GEOS_Intersection_r)(CR_GEOS_Handle, CR_GEOS_Geometry,
                                                    CR_GEOS_Geometry);
+typedef CR_GEOS_Geometry (*CR_GEOS_SymDifference_r)(CR_GEOS_Handle, CR_GEOS_Geometry,
+                                                    CR_GEOS_Geometry);
 typedef CR_GEOS_Geometry (*CR_GEOS_PointOnSurface_r)(CR_GEOS_Handle, CR_GEOS_Geometry);
 
 typedef CR_GEOS_Geometry (*CR_GEOS_Interpolate_r)(CR_GEOS_Handle, CR_GEOS_Geometry, double);
@@ -185,6 +187,7 @@ struct CR_GEOS {
   CR_GEOS_Union_r GEOSUnion_r;
   CR_GEOS_PointOnSurface_r GEOSPointOnSurface_r;
   CR_GEOS_Intersection_r GEOSIntersection_r;
+  CR_GEOS_SymDifference_r GEOSSymDifference_r;
 
   CR_GEOS_Interpolate_r GEOSInterpolate_r;
 
@@ -266,6 +269,7 @@ struct CR_GEOS {
     INIT(GEOSUnion_r);
     INIT(GEOSPointOnSurface_r);
     INIT(GEOSIntersection_r);
+    INIT(GEOSSymDifference_r);
     INIT(GEOSInterpolate_r);
     INIT(GEOSDistance_r);
     INIT(GEOSCovers_r);
@@ -798,6 +802,33 @@ CR_GEOS_Status CR_GEOS_Intersection(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_Slice
   return toGEOSString(error.data(), error.length());
 }
 
+CR_GEOS_Status CR_GEOS_SymDifference(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_Slice b,
+                                    CR_GEOS_String* symdifferenceEWKB) {
+  std::string error;
+  auto handle = initHandleWithErrorBuffer(lib, &error);
+  *symdifferenceEWKB = {.data = NULL, .len = 0};
+
+  auto geomA = CR_GEOS_GeometryFromSlice(lib, handle, a);
+  auto geomB = CR_GEOS_GeometryFromSlice(lib, handle, b);
+  if (geomA != nullptr && geomB != nullptr) {
+    auto symdifferenceGeom = lib->GEOSSymDifference_r(handle, geomA, geomB);
+    if (symdifferenceGeom != nullptr) {
+      auto srid = lib->GEOSGetSRID_r(handle, geomA);
+      CR_GEOS_writeGeomToEWKB(lib, handle, symdifferenceGeom, symdifferenceEWKB, srid);
+      lib->GEOSGeom_destroy_r(handle, symdifferenceGeom);
+    }
+  }
+  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());
+}
+
 //
 // Linear Reference
 //

pkg/geo/geos/geos.go

@@ -417,6 +417,19 @@ func Union(a geopb.EWKB, b geopb.EWKB) (geopb.EWKB, error) {
 	return cStringToSafeGoBytes(cEWKB), nil
 }
 
+// SymDifference returns an EWKB which is the symmetric difference of shapes A and B.
+func SymDifference(a geopb.EWKB, b geopb.EWKB) (geopb.EWKB, error) {
+	g, err := ensureInitInternal()
+	if err != nil {
+		return nil, err
+	}
+	var cEWKB C.CR_GEOS_String
+	if err := statusToError(C.CR_GEOS_SymDifference(g, goToCSlice(a), goToCSlice(b), &cEWKB)); err != nil {
+		return nil, err
+	}
+	return cStringToSafeGoBytes(cEWKB), nil
+}
+
 // InterpolateLine returns the point along the given LineString which is at
 // a given distance from starting point.
 // Note: For distance less than 0 it returns start point similarly for distance

pkg/geo/geos/geos.h

@@ -103,6 +103,8 @@ CR_GEOS_Status CR_GEOS_PointOnSurface(CR_GEOS* lib, CR_GEOS_Slice a,
                                       CR_GEOS_String* pointOnSurfaceEWKB);
 CR_GEOS_Status CR_GEOS_Intersection(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_Slice b,
                                     CR_GEOS_String* intersectionEWKB);
+CR_GEOS_Status CR_GEOS_SymDifference(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_Slice b,
+                                     CR_GEOS_String* symdifferenceEWKB);
 CR_GEOS_Status CR_GEOS_SharedPaths(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_Slice b, CR_GEOS_String* ret);
 
 //