Add PreparedGeometry for efficient repeated spatial predicate evaluation

CHANGELOG.md

@@ -2,6 +2,13 @@
 
 ## Unreleased
 
+- Add `PreparedGeometry` type and `Prepare` function for efficient repeated
+  spatial predicate evaluation. When a geometry is prepared, it caches spatial
+  indices so that subsequent calls to `Intersects`, `Contains`,
+  `ContainsProperly`, `CoveredBy`, `Covers`, `Crosses`, `Disjoint`, `Overlaps`,
+  `Touches`, and `Within` against different test geometries are fast. The
+  implementation is based on a port of JTS.
+
 - Add `Buffer` function that computes the buffer of a geometry at a given
   distance. Options are available for controlling quad segments, end cap style
   (round, flat, square), join style (round, mitre, bevel), single-sided mode,

geom/alg_prepared.go

@@ -0,0 +1,96 @@
+package geom
+
+import "github.com/peterstace/simplefeatures/internal/jtsport/jts"
+
+// PreparedGeometry is a geometry that has been preprocessed to efficiently
+// enable evaluation of spatial predicates against many other geometries.
+//
+// It is created by calling [Prepare] with the geometry to be prepared. The
+// prepared geometry caches spatial indices and other structures so that
+// repeated predicate evaluations (e.g. [PreparedGeometry.Intersects],
+// [PreparedGeometry.Contains]) against different test geometries are fast.
+type PreparedGeometry struct {
+	prep jts.GeomPrep_PreparedGeometry
+}
+
+// Prepare preprocesses a geometry for efficient repeated predicate evaluation.
+func Prepare(g Geometry) (PreparedGeometry, error) {
+	var pg PreparedGeometry
+	err := catch(func() error {
+		jtsG, err := toJTS(g)
+		if err != nil {
+			return err
+		}
+		pg.prep = jts.GeomPrep_PreparedGeometryFactory_Prepare(jtsG)
+		return nil
+	})
+	return pg, err
+}
+
+func toJTS(g Geometry) (*jts.Geom_Geometry, error) {
+	jtsG, err := jts.Io_NewWKBReader().ReadBytes(g.AsBinary())
+	if err != nil {
+		return nil, wrap(err, "converting geometry to JTS")
+	}
+	return jtsG, nil
+}
+
+func (p PreparedGeometry) eval(g Geometry, pred func(*jts.Geom_Geometry) bool) (bool, error) {
+	var result bool
+	err := catch(func() error {
+		jtsG, err := toJTS(g)
+		if err != nil {
+			return err
+		}
+		result = pred(jtsG)
+		return nil
+	})
+	return result, err
+}
+
+// Intersects reports whether the prepared geometry intersects g.
+func (p PreparedGeometry) Intersects(g Geometry) (bool, error) {
+	return p.eval(g, p.prep.Intersects)
+}
+
+// Contains reports whether the prepared geometry contains g.
+func (p PreparedGeometry) Contains(g Geometry) (bool, error) {
+	return p.eval(g, p.prep.Contains)
+}
+
+// ContainsProperly reports whether the prepared geometry properly contains g.
+// A geometry properly contains another if it contains it and the other
+// geometry has no points on the boundary of the prepared geometry.
+func (p PreparedGeometry) ContainsProperly(g Geometry) (bool, error) {
+	return p.eval(g, p.prep.ContainsProperly)
+}
+
+// CoveredBy reports whether the prepared geometry is covered by g.
+func (p PreparedGeometry) CoveredBy(g Geometry) (bool, error) {
+	return p.eval(g, p.prep.CoveredBy)
+}
+
+// Covers reports whether the prepared geometry covers g.
+func (p PreparedGeometry) Covers(g Geometry) (bool, error) {
+	return p.eval(g, p.prep.Covers)
+}
+
+// Disjoint reports whether the prepared geometry is disjoint from g.
+func (p PreparedGeometry) Disjoint(g Geometry) (bool, error) {
+	return p.eval(g, p.prep.Disjoint)
+}
+
+// Overlaps reports whether the prepared geometry overlaps g.
+func (p PreparedGeometry) Overlaps(g Geometry) (bool, error) {
+	return p.eval(g, p.prep.Overlaps)
+}
+
+// Touches reports whether the prepared geometry touches g.
+func (p PreparedGeometry) Touches(g Geometry) (bool, error) {
+	return p.eval(g, p.prep.Touches)
+}
+
+// Within reports whether the prepared geometry is within g.
+func (p PreparedGeometry) Within(g Geometry) (bool, error) {
+	return p.eval(g, p.prep.Within)
+}

geom/alg_prepared_test.go

@@ -0,0 +1,249 @@
+package geom_test
+
+import (
+	"strconv"
+	"testing"
+
+	"github.com/peterstace/simplefeatures/geom"
+	"github.com/peterstace/simplefeatures/internal/test"
+)
+
+func TestPreparedGeometry(t *testing.T) {
+	for i, tt := range []struct {
+		wkt1, wkt2 string
+	}{
+		// Point vs point
+		{"POINT(1 2)", "POINT(1 2)"},
+		{"POINT(1 2)", "POINT(3 4)"},
+
+		// Point vs linestring
+		{"POINT(0.5 0)", "LINESTRING(0 0,1 0)"},
+		{"POINT(0 0)", "LINESTRING(0 0,1 0)"},
+		{"POINT(0 1)", "LINESTRING(0 0,1 0)"},
+
+		// Point vs polygon
+		{"POINT(0.5 0.5)", "POLYGON((0 0,1 0,1 1,0 1,0 0))"},
+		{"POINT(0 0)", "POLYGON((0 0,1 0,1 1,0 1,0 0))"},
+		{"POINT(5 5)", "POLYGON((0 0,1 0,1 1,0 1,0 0))"},
+
+		// Linestring vs linestring
+		{"LINESTRING(0 0,1 1)", "LINESTRING(0 1,1 0)"},
+		{"LINESTRING(0 0,1 0)", "LINESTRING(1 0,2 0)"},
+		{"LINESTRING(0 0,2 0)", "LINESTRING(1 0,3 0)"},
+		{"LINESTRING(0 0,1 0)", "LINESTRING(2 0,3 0)"},
+
+		// Linestring vs polygon
+		{"LINESTRING(0 0,2 2)", "POLYGON((0 0,1 0,1 1,0 1,0 0))"},
+		{"LINESTRING(0.25 0.25,0.75 0.75)", "POLYGON((0 0,1 0,1 1,0 1,0 0))"},
+		{"LINESTRING(5 5,6 6)", "POLYGON((0 0,1 0,1 1,0 1,0 0))"},
+		{"LINESTRING(0 0,1 0)", "POLYGON((0 0,1 0,1 1,0 1,0 0))"},
+
+		// Polygon vs polygon
+		{"POLYGON((0 0,3 0,3 3,0 3,0 0))", "POLYGON((1 1,2 1,2 2,1 2,1 1))"},
+		{"POLYGON((1 1,2 1,2 2,1 2,1 1))", "POLYGON((0 0,3 0,3 3,0 3,0 0))"},
+		{"POLYGON((0 0,2 0,2 2,0 2,0 0))", "POLYGON((1 0,3 0,3 2,1 2,1 0))"},
+		{"POLYGON((0 0,1 0,1 1,0 1,0 0))", "POLYGON((1 0,2 0,2 1,1 1,1 0))"},
+		{"POLYGON((0 0,1 0,1 1,0 1,0 0))", "POLYGON((2 2,3 2,3 3,2 3,2 2))"},
+		{"POLYGON((0 0,1 0,1 1,0 1,0 0))", "POLYGON((0 0,1 0,1 1,0 1,0 0))"},
+
+		// Empty geometries
+		{"POINT EMPTY", "POINT EMPTY"},
+		{"POINT EMPTY", "POINT(1 2)"},
+		{"POLYGON((0 0,1 0,1 1,0 1,0 0))", "POINT EMPTY"},
+
+		// Polygon with hole
+		{"POLYGON((0 0,10 0,10 10,0 10,0 0),(3 3,7 3,7 7,3 7,3 3))", "POINT(5 5)"},
+		{"POLYGON((0 0,10 0,10 10,0 10,0 0),(3 3,7 3,7 7,3 7,3 3))", "POINT(1 1)"},
+
+		// Multi-geometries
+		{"MULTIPOINT((0 0),(1 1))", "POINT(0 0)"},
+		{"MULTILINESTRING((0 0,1 0),(0 1,1 1))", "POINT(0.5 0)"},
+		{"MULTIPOLYGON(((0 0,1 0,1 1,0 1,0 0)),((2 2,3 2,3 3,2 3,2 2)))", "POINT(0.5 0.5)"},
+
+		// Geometry collection as prepared geometry
+		{"GEOMETRYCOLLECTION(POINT(0 0),LINESTRING(1 1,2 2),POLYGON((3 3,6 3,6 6,3 6,3 3)))", "POINT(0 0)"},
+		{"GEOMETRYCOLLECTION(POINT(0 0),LINESTRING(1 1,2 2),POLYGON((3 3,6 3,6 6,3 6,3 3)))", "POINT(4 4)"},
+		{"GEOMETRYCOLLECTION(POINT(0 0),LINESTRING(1 1,2 2),POLYGON((3 3,6 3,6 6,3 6,3 3)))", "POINT(9 9)"},
+		{"GEOMETRYCOLLECTION(POINT(0 0),LINESTRING(1 1,2 2),POLYGON((3 3,6 3,6 6,3 6,3 3)))", "LINESTRING(4 4,5 5)"},
+
+		// Geometry collection as test geometry
+		{"POLYGON((0 0,10 0,10 10,0 10,0 0))", "GEOMETRYCOLLECTION(POINT(1 1),LINESTRING(2 2,3 3))"},
+		{"POLYGON((0 0,10 0,10 10,0 10,0 0))", "GEOMETRYCOLLECTION(POINT(1 1),POINT(20 20))"},
+
+		// Both geometry collections
+		{"GEOMETRYCOLLECTION(POINT(0 0),POLYGON((1 1,4 1,4 4,1 4,1 1)))", "GEOMETRYCOLLECTION(POINT(2 2),LINESTRING(2 2,3 3))"},
+		{"GEOMETRYCOLLECTION(POINT(0 0),POLYGON((1 1,4 1,4 4,1 4,1 1)))", "GEOMETRYCOLLECTION(POINT(20 20),LINESTRING(20 20,30 30))"},
+
+		// Empty geometry collection
+		{"GEOMETRYCOLLECTION EMPTY", "POINT(1 1)"},
+		{"POINT(1 1)", "GEOMETRYCOLLECTION EMPTY"},
+		{"GEOMETRYCOLLECTION EMPTY", "GEOMETRYCOLLECTION EMPTY"},
+	} {
+		t.Run(strconv.Itoa(i), func(t *testing.T) {
+			g1 := test.FromWKT(t, tt.wkt1)
+			g2 := test.FromWKT(t, tt.wkt2)
+
+			pg, err := geom.Prepare(g1)
+			test.NoErr(t, err)
+
+			type predicate struct {
+				name string
+				got  func() (bool, error)
+				want func() (bool, error)
+			}
+
+			predicates := []predicate{
+				{
+					name: "Intersects",
+					got:  func() (bool, error) { return pg.Intersects(g2) },
+					want: func() (bool, error) { return geom.Intersects(g1, g2), nil },
+				},
+				{
+					name: "Contains",
+					got:  func() (bool, error) { return pg.Contains(g2) },
+					want: func() (bool, error) { return geom.Contains(g1, g2) },
+				},
+				{
+					name: "CoveredBy",
+					got:  func() (bool, error) { return pg.CoveredBy(g2) },
+					want: func() (bool, error) { return geom.CoveredBy(g1, g2) },
+				},
+				{
+					name: "Covers",
+					got:  func() (bool, error) { return pg.Covers(g2) },
+					want: func() (bool, error) { return geom.Covers(g1, g2) },
+				},
+				{
+					name: "Disjoint",
+					got:  func() (bool, error) { return pg.Disjoint(g2) },
+					want: func() (bool, error) { return geom.Disjoint(g1, g2) },
+				},
+				{
+					name: "Overlaps",
+					got:  func() (bool, error) { return pg.Overlaps(g2) },
+					want: func() (bool, error) { return geom.Overlaps(g1, g2) },
+				},
+				{
+					name: "Touches",
+					got:  func() (bool, error) { return pg.Touches(g2) },
+					want: func() (bool, error) { return geom.Touches(g1, g2) },
+				},
+				{
+					name: "Within",
+					got:  func() (bool, error) { return pg.Within(g2) },
+					want: func() (bool, error) { return geom.Within(g1, g2) },
+				},
+			}
+
+			for _, pred := range predicates {
+				t.Run(pred.name, func(t *testing.T) {
+					got, gotErr := pred.got()
+					test.NoErr(t, gotErr)
+					want, wantErr := pred.want()
+					test.NoErr(t, wantErr)
+					test.Eq(t, got, want)
+				})
+			}
+		})
+	}
+}
+
+func TestPreparedGeometryContainsProperly(t *testing.T) {
+	for i, tt := range []struct {
+		wkt1, wkt2 string
+		want       bool
+	}{
+		// Point in polygon interior: true
+		{
+			wkt1: "POLYGON((0 0,1 0,1 1,0 1,0 0))",
+			wkt2: "POINT(0.5 0.5)",
+			want: true,
+		},
+		// Point on polygon boundary: false (key distinction from Contains)
+		{
+			wkt1: "POLYGON((0 0,1 0,1 1,0 1,0 0))",
+			wkt2: "POINT(0 0)",
+			want: false,
+		},
+		// Point on polygon edge: false
+		{
+			wkt1: "POLYGON((0 0,1 0,1 1,0 1,0 0))",
+			wkt2: "POINT(0.5 0)",
+			want: false,
+		},
+		// Polygon properly containing another polygon
+		{
+			wkt1: "POLYGON((0 0,10 0,10 10,0 10,0 0))",
+			wkt2: "POLYGON((1 1,2 1,2 2,1 2,1 1))",
+			want: true,
+		},
+		// Equal polygons: false (boundary points shared)
+		{
+			wkt1: "POLYGON((0 0,1 0,1 1,0 1,0 0))",
+			wkt2: "POLYGON((0 0,1 0,1 1,0 1,0 0))",
+			want: false,
+		},
+		// Both empty
+		{
+			wkt1: "POINT EMPTY",
+			wkt2: "POINT EMPTY",
+			want: false,
+		},
+		// One empty
+		{
+			wkt1: "POLYGON((0 0,1 0,1 1,0 1,0 0))",
+			wkt2: "POINT EMPTY",
+			want: false,
+		},
+		// Disjoint
+		{
+			wkt1: "POLYGON((0 0,1 0,1 1,0 1,0 0))",
+			wkt2: "POINT(5 5)",
+			want: false,
+		},
+		// Polygon with shared boundary edge: false
+		{
+			wkt1: "POLYGON((0 0,2 0,2 2,0 2,0 0))",
+			wkt2: "POLYGON((1 0,2 0,2 1,1 1,1 0))",
+			want: false,
+		},
+	} {
+		t.Run(strconv.Itoa(i), func(t *testing.T) {
+			g1 := test.FromWKT(t, tt.wkt1)
+			g2 := test.FromWKT(t, tt.wkt2)
+
+			pg, err := geom.Prepare(g1)
+			test.NoErr(t, err)
+
+			got, err := pg.ContainsProperly(g2)
+			test.NoErr(t, err)
+			test.Eq(t, got, tt.want)
+		})
+	}
+}
+
+func TestPreparedGeometryMultipleEvaluations(t *testing.T) {
+	pg, err := geom.Prepare(test.FromWKT(t, "POLYGON((0 0,10 0,10 10,0 10,0 0))"))
+	test.NoErr(t, err)
+
+	tests := []struct {
+		wkt  string
+		want bool
+	}{
+		{"POINT(5 5)", true},
+		{"POINT(15 15)", false},
+		{"LINESTRING(1 1,2 2)", true},
+		{"LINESTRING(11 11,12 12)", false},
+		{"POLYGON((1 1,2 1,2 2,1 2,1 1))", true},
+		{"POLYGON((20 20,21 20,21 21,20 21,20 20))", false},
+	}
+
+	for i, tt := range tests {
+		t.Run(strconv.Itoa(i), func(t *testing.T) {
+			g := test.FromWKT(t, tt.wkt)
+			got, err := pg.Intersects(g)
+			test.NoErr(t, err)
+			test.Eq(t, got, tt.want)
+		})
+	}
+}