DRAFT: JTS Port

.ci/compose-unit.yaml

@@ -2,7 +2,12 @@ services:
   tests:
     image: golang:1.18
     working_dir: /go/src/github.com/peterstace/simplefeatures
-    entrypoint: go test -covermode=count -coverprofile=coverage.out -test.count=1 -test.run=. ./geom ./rtree
+    entrypoint: >-
+      go test -covermode=count -coverprofile=coverage.out -test.count=1 -test.run=.
+      ./carto
+      ./geom
+      ./internal/jtsport/...
+      ./rtree
     volumes:
       - ..:/go/src/github.com/peterstace/simplefeatures
     environment:

.golangci.yaml

@@ -1,5 +1,7 @@
 run:
   timeout: 5m
+  skip-dirs:
+    - internal/jtsport/jts
 
 issues:
   exclude-rules:

CLAUDE.md

@@ -0,0 +1 @@
+# CLAUDE.md

THIRD_PARTY_LICENSES

@@ -0,0 +1,4 @@
+This software includes code derived from JTS (Java Topology Suite)
+Copyright (c) 2016-2024 LocationTech and Contributors
+Licensed under the Eclipse Distribution License 1.0
+See: https://github.com/locationtech/jts.

geom/alg_relate.go

@@ -1,5 +1,7 @@
 package geom
 
+import "github.com/peterstace/simplefeatures/internal/jtsport/jts"
+
 // Relate calculates the DE-9IM matrix between two [Geometry] values,
 // describing how they relate to each other.
 //
@@ -21,42 +23,64 @@ package geom
 // The matrix is represented by a 9 character string, with entries in row-major
 // order (i.e. entries are ordered II, IB, IE, BI, BB, BE, EI, EB, EE).
 func Relate(a, b Geometry) (string, error) {
-	// TODO: don't need to return an error from this function
 	if a.IsEmpty() || b.IsEmpty() {
-		im := newMatrix()
-		im.set(imExterior, imExterior, '2')
-		if a.IsEmpty() && b.IsEmpty() {
-			return im.code(), nil
-		}
+		return relateWithEmptyInput(a, b), nil
+	}
+	// TODO: Optimize for when the envelopes don't intersect.
+	return jtsRelateNG(a, b)
+}
 
-		var flip bool
-		nonEmpty := b
-		if b.IsEmpty() {
-			nonEmpty = a
-			flip = true
-		}
-		switch nonEmpty.Dimension() {
-		case 0:
-			im.set(imExterior, imInterior, '0')
-			im.set(imExterior, imBoundary, 'F')
-		case 1:
-			im.set(imExterior, imInterior, '1')
-			if !nonEmpty.Boundary().IsEmpty() {
-				im.set(imExterior, imBoundary, '0')
-			}
-		case 2:
-			im.set(imExterior, imInterior, '2')
-			im.set(imExterior, imBoundary, '1')
-		}
-		if flip {
-			im.transpose()
+// relateWithEmptyInput computes the DE-9IM matrix when at least one input is empty.
+func relateWithEmptyInput(a, b Geometry) string {
+	im := newMatrix()
+	im.set(imExterior, imExterior, '2')
+	if a.IsEmpty() && b.IsEmpty() {
+		return im.code()
+	}
+
+	var flip bool
+	nonEmpty := b
+	if b.IsEmpty() {
+		nonEmpty = a
+		flip = true
+	}
+	switch nonEmpty.Dimension() {
+	case 0:
+		im.set(imExterior, imInterior, '0')
+		im.set(imExterior, imBoundary, 'F')
+	case 1:
+		im.set(imExterior, imInterior, '1')
+		if !nonEmpty.Boundary().IsEmpty() {
+			im.set(imExterior, imBoundary, '0')
 		}
-		return im.code(), nil
+	case 2:
+		im.set(imExterior, imInterior, '2')
+		im.set(imExterior, imBoundary, '1')
 	}
+	if flip {
+		im.transpose()
+	}
+	return im.code()
+}
 
-	overlay := newDCELFromGeometries(a, b)
-	im := overlay.extractIntersectionMatrix()
-	return im.code(), nil
+// jtsRelateNG invokes the JTS port's RelateNG operation.
+func jtsRelateNG(a, b Geometry) (string, error) {
+	var result string
+	err := catch(func() error {
+		wkbReader := jts.Io_NewWKBReader()
+		jtsA, err := wkbReader.ReadBytes(a.AsBinary())
+		if err != nil {
+			return wrap(err, "converting geometry A to JTS")
+		}
+		jtsB, err := wkbReader.ReadBytes(b.AsBinary())
+		if err != nil {
+			return wrap(err, "converting geometry B to JTS")
+		}
+		im := jts.OperationRelateng_RelateNG_RelateMatrix(jtsA, jtsB)
+		result = im.String()
+		return validateIntersectionMatrix(result)
+	})
+	return result, err
 }
 
 func relateMatchesAnyPattern(a, b Geometry, patterns ...string) (bool, error) {

geom/alg_set_op.go

@@ -1,20 +1,37 @@
 package geom
 
+import "github.com/peterstace/simplefeatures/internal/jtsport/jts"
+
+// TODO: Optimise operations by comparing bounding boxes first.
+
+// hasGC determines if either argument is a GeometryCollection. The JTS port
+// doesn't have full support for GeometryCollections, so we need to handle them
+// specially.
+//
+// Specifically, the JTS port has the following restrictions for binary overlay
+// operations:
+//
+// 1. All input elements must have the same dimension.
+//
+// 2. Polygons in a GeometryCollection must not overlap.
+func hasGC(a, b Geometry) bool {
+	return a.IsGeometryCollection() || b.IsGeometryCollection()
+}
+
 // Union returns a geometry that represents the parts from either geometry A or
 // geometry B (or both). An error may be returned in pathological cases of
 // numerical degeneracy.
 func Union(a, b Geometry) (Geometry, error) {
-	if a.IsEmpty() && b.IsEmpty() {
-		return Geometry{}, nil
+	if hasGC(a, b) {
+		return gcAwareUnion(a, b)
 	}
-	if a.IsEmpty() {
-		return UnaryUnion(b)
-	}
-	if b.IsEmpty() {
-		return UnaryUnion(a)
-	}
-	g, err := setOp(a, or, b)
-	return g, wrap(err, "executing union")
+	return jtsOverlayOp(a, b, jts.OperationOverlayng_OverlayNG_UNION)
+}
+
+func gcAwareUnion(a, b Geometry) (Geometry, error) {
+	// UnaryUnion supports arbitrary GeometryCollections.
+	gc := NewGeometryCollection([]Geometry{a, b})
+	return UnaryUnion(gc.AsGeometry())
 }
 
 // Intersection returns a geometry that represents the parts that are common to
@@ -24,8 +41,64 @@ func Intersection(a, b Geometry) (Geometry, error) {
 	if a.IsEmpty() || b.IsEmpty() {
 		return Geometry{}, nil
 	}
-	g, err := setOp(a, and, b)
-	return g, wrap(err, "executing intersection")
+	if hasGC(a, b) {
+		return gcAwareIntersection(a, b)
+	}
+	return jtsOverlayOp(a, b, jts.OperationOverlayng_OverlayNG_INTERSECTION)
+}
+
+func gcAwareIntersection(a, b Geometry) (Geometry, error) {
+	partsA, partsB, err := prepareOverlayInputParts(a, b)
+	if err != nil {
+		return Geometry{}, err
+	}
+
+	// The total result is the union of the intersections across the Cartesian
+	// product of parts.
+	var results []Geometry
+	for _, partA := range partsA {
+		for _, partB := range partsB {
+			result, err := jtsOverlayOp(partA, partB, jts.OperationOverlayng_OverlayNG_INTERSECTION)
+			if err != nil {
+				return Geometry{}, err
+			}
+			results = append(results, result)
+		}
+	}
+	return UnaryUnion(NewGeometryCollection(results).AsGeometry())
+}
+
+func explodeGeometryCollections(dst []Geometry, g Geometry) []Geometry {
+	if gc, ok := g.AsGeometryCollection(); ok {
+		for i := 0; i < gc.NumGeometries(); i++ {
+			dst = explodeGeometryCollections(dst, gc.GeometryN(i))
+		}
+		return dst
+	}
+	return append(dst, g)
+}
+
+func prepareOverlayInputParts(a, b Geometry) ([]Geometry, []Geometry, error) {
+	// Normalize GC inputs by unioning their parts.
+	if a.IsGeometryCollection() {
+		var err error
+		a, err = UnaryUnion(a)
+		if err != nil {
+			return nil, nil, err
+		}
+	}
+	if b.IsGeometryCollection() {
+		var err error
+		b, err = UnaryUnion(b)
+		if err != nil {
+			return nil, nil, err
+		}
+	}
+
+	// Extract non-GC parts from each input.
+	partsA := explodeGeometryCollections(nil, a)
+	partsB := explodeGeometryCollections(nil, b)
+	return partsA, partsB, nil
 }
 
 // Difference returns a geometry that represents the parts of input geometry A
@@ -35,11 +108,57 @@ func Difference(a, b Geometry) (Geometry, error) {
 	if a.IsEmpty() {
 		return Geometry{}, nil
 	}
-	if b.IsEmpty() {
-		return UnaryUnion(a)
+	if hasGC(a, b) {
+		return gcAwareDifference(a, b)
 	}
-	g, err := setOp(a, andNot, b)
-	return g, wrap(err, "executing difference")
+	return jtsOverlayOp(a, b, jts.OperationOverlayng_OverlayNG_DIFFERENCE)
+}
+
+func gcAwareDifference(a, b Geometry) (Geometry, error) {
+	partsA, partsB, err := prepareOverlayInputParts(a, b)
+	if err != nil {
+		return Geometry{}, err
+	}
+
+	// The total result is the union of each part of A after each part of B has
+	// been removed (sequentially).
+	var results []Geometry
+	for _, partA := range partsA {
+		result := partA
+		for _, partB := range partsB {
+			var err error
+			result, err = jtsOverlayOp(result, partB, jts.OperationOverlayng_OverlayNG_DIFFERENCE)
+			if err != nil {
+				return Geometry{}, err
+			}
+			if result.IsEmpty() {
+				break
+			}
+		}
+		results = append(results, result)
+	}
+	return UnaryUnion(NewGeometryCollection(results).AsGeometry())
+}
+
+// jtsOverlayOp invokes the JTS port's overlay operation with the given opCode.
+func jtsOverlayOp(a, b Geometry, opCode int) (Geometry, error) {
+	var result Geometry
+	err := catch(func() error {
+		wkbReader := jts.Io_NewWKBReader()
+		jtsA, err := wkbReader.ReadBytes(a.AsBinary())
+		if err != nil {
+			return wrap(err, "converting geometry A to JTS")
+		}
+		jtsB, err := wkbReader.ReadBytes(b.AsBinary())
+		if err != nil {
+			return wrap(err, "converting geometry B to JTS")
+		}
+		jtsResult := jts.OperationOverlayng_OverlayNGRobust_Overlay(jtsA, jtsB, opCode)
+		wkbWriter := jts.Io_NewWKBWriter()
+		result, err = UnmarshalWKB(wkbWriter.Write(jtsResult), NoValidate{})
+		return wrap(err, "converting JTS overlay result to simplefeatures")
+	})
+	return result, err
 }
 
 // SymmetricDifference returns a geometry that represents the parts of geometry
@@ -55,14 +174,25 @@ func SymmetricDifference(a, b Geometry) (Geometry, error) {
 	if b.IsEmpty() {
 		return UnaryUnion(a)
 	}
-	g, err := setOp(a, xor, b)
-	return g, wrap(err, "executing symmetric difference")
+
+	diffAB, err := Difference(a, b)
+	if err != nil {
+		return Geometry{}, err
+	}
+	diffBA, err := Difference(b, a)
+	if err != nil {
+		return Geometry{}, err
+	}
+	return Union(diffAB, diffBA)
 }
 
 // UnaryUnion is a single input variant of the Union function, unioning
 // together the components of the input geometry.
 func UnaryUnion(g Geometry) (Geometry, error) {
-	return setOp(g, or, Geometry{})
+	if g.IsEmpty() {
+		return Geometry{}, nil
+	}
+	return jtsUnaryUnion(g)
 }
 
 // UnionMany unions together the input geometries.
@@ -71,19 +201,19 @@ func UnionMany(gs []Geometry) (Geometry, error) {
 	return UnaryUnion(gc.AsGeometry())
 }
 
-func setOp(a Geometry, include func([2]bool) bool, b Geometry) (Geometry, error) {
-	overlay := newDCELFromGeometries(a, b)
-	g, err := overlay.extractGeometry(include)
-	if err != nil {
-		return Geometry{}, wrap(err, "internal error extracting geometry")
-	}
-	if err := g.Validate(); err != nil {
-		return Geometry{}, wrap(err, "invalid geometry produced by overlay")
-	}
-	return g, nil
+// jtsUnaryUnion invokes the JTS port's unary union operation.
+func jtsUnaryUnion(g Geometry) (Geometry, error) {
+	var result Geometry
+	err := catch(func() error {
+		wkbReader := jts.Io_NewWKBReader()
+		jtsG, err := wkbReader.ReadBytes(g.AsBinary())
+		if err != nil {
+			return wrap(err, "converting geometry to JTS")
+		}
+		jtsResult := jts.OperationOverlayng_OverlayNGRobust_Union(jtsG)
+		wkbWriter := jts.Io_NewWKBWriter()
+		result, err = UnmarshalWKB(wkbWriter.Write(jtsResult), NoValidate{})
+		return wrap(err, "converting JTS union result to simplefeatures")
+	})
+	return result, err
 }
-
-func or(b [2]bool) bool     { return b[0] || b[1] }
-func and(b [2]bool) bool    { return b[0] && b[1] }
-func xor(b [2]bool) bool    { return b[0] != b[1] }
-func andNot(b [2]bool) bool { return b[0] && !b[1] }