Skip to content

Commit abb3854

Browse files
authored
Merge pull request #705 from peterstace/rect_clip
Add ClipByRect operation to geos package
2 parents c592301 + 5a722f1 commit abb3854

4 files changed

Lines changed: 123 additions & 0 deletions

File tree

CHANGELOG.md

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,11 @@
22

33
## Unreleased
44

5+
- Add `geos.ClipByRect` function that clips a geometry to an axis-aligned
6+
rectangle (defined by a `geom.Envelope`). This wraps the GEOS
7+
`GEOSClipByRect` operation, which is faster than computing a full
8+
`Intersection` with a rectangular polygon.
9+
510
## v0.58.0
611

712
2026-02-15

geos/entrypoints.go

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -331,3 +331,16 @@ func UnaryUnion(g geom.Geometry) (geom.Geometry, error) {
331331
func ConcaveHull(g geom.Geometry, concavenessRatio float64, allowHoles bool) (geom.Geometry, error) {
332332
return rawgeos.ConcaveHull(g, concavenessRatio, allowHoles)
333333
}
334+
335+
// ClipByRect clips a geometry to an axis-aligned rectangle defined by the
336+
// given [geom.Envelope]. If the envelope is empty, then an empty
337+
// [geom.GeometryCollection] is returned.
338+
//
339+
// The validity of the result is not checked.
340+
func ClipByRect(g geom.Geometry, rect geom.Envelope) (geom.Geometry, error) {
341+
lo, hi, ok := rect.MinMaxXYs()
342+
if !ok {
343+
return geom.Geometry{}, nil
344+
}
345+
return rawgeos.ClipByRect(g, lo.X, lo.Y, hi.X, hi.Y)
346+
}

geos/entrypoints_test.go

Lines changed: 83 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1016,3 +1016,86 @@ func TestCoverageIsValid(t *testing.T) {
10161016
})
10171017
}
10181018
}
1019+
1020+
func TestClipByRect(t *testing.T) {
1021+
for _, tc := range []struct {
1022+
name string
1023+
input string
1024+
rect geom.Envelope
1025+
want string
1026+
}{
1027+
{
1028+
name: "polygon fully inside rect",
1029+
input: "POLYGON((1 1,1 2,2 2,2 1,1 1))",
1030+
rect: geom.NewEnvelope(geom.XY{X: 0, Y: 0}, geom.XY{X: 3, Y: 3}),
1031+
want: "POLYGON((1 1,1 2,2 2,2 1,1 1))",
1032+
},
1033+
{
1034+
name: "polygon partially overlapping rect",
1035+
input: "POLYGON((0 0,0 4,4 4,4 0,0 0))",
1036+
rect: geom.NewEnvelope(geom.XY{X: 1, Y: 1}, geom.XY{X: 3, Y: 3}),
1037+
want: "POLYGON((1 1,1 3,3 3,3 1,1 1))",
1038+
},
1039+
{
1040+
name: "polygon fully outside rect",
1041+
input: "POLYGON((0 0,0 1,1 1,1 0,0 0))",
1042+
rect: geom.NewEnvelope(geom.XY{X: 5, Y: 5}, geom.XY{X: 6, Y: 6}),
1043+
want: "GEOMETRYCOLLECTION EMPTY",
1044+
},
1045+
{
1046+
name: "linestring clipped by rect",
1047+
input: "LINESTRING(0 0,4 4)",
1048+
rect: geom.NewEnvelope(geom.XY{X: 1, Y: 1}, geom.XY{X: 3, Y: 3}),
1049+
want: "LINESTRING(1 1,3 3)",
1050+
},
1051+
{
1052+
name: "point inside rect",
1053+
input: "POINT(2 2)",
1054+
rect: geom.NewEnvelope(geom.XY{X: 1, Y: 1}, geom.XY{X: 3, Y: 3}),
1055+
want: "POINT(2 2)",
1056+
},
1057+
{
1058+
name: "point outside rect",
1059+
input: "POINT(0 0)",
1060+
rect: geom.NewEnvelope(geom.XY{X: 1, Y: 1}, geom.XY{X: 3, Y: 3}),
1061+
want: "GEOMETRYCOLLECTION EMPTY",
1062+
},
1063+
{
1064+
name: "empty input geometry",
1065+
input: "GEOMETRYCOLLECTION EMPTY",
1066+
rect: geom.NewEnvelope(geom.XY{X: 0, Y: 0}, geom.XY{X: 1, Y: 1}),
1067+
want: "GEOMETRYCOLLECTION EMPTY",
1068+
},
1069+
{
1070+
name: "u-shaped polygon clipped through both arms produces multipolygon",
1071+
input: "POLYGON((0 0,4 0,4 3,3 3,3 1,1 1,1 3,0 3,0 0))",
1072+
rect: geom.NewEnvelope(geom.XY{X: 0, Y: 2}, geom.XY{X: 4, Y: 4}),
1073+
want: "MULTIPOLYGON(((0 2,0 3,1 3,1 2,0 2)),((3 2,3 3,4 3,4 2,3 2)))",
1074+
},
1075+
{
1076+
name: "polygon with hole inside rect",
1077+
input: "POLYGON((0 0,0 6,6 6,6 0,0 0),(2 2,4 2,4 4,2 4,2 2))",
1078+
rect: geom.NewEnvelope(geom.XY{X: 1, Y: 1}, geom.XY{X: 5, Y: 5}),
1079+
want: "POLYGON((1 1,1 5,5 5,5 1,1 1),(2 2,4 2,4 4,2 4,2 2))",
1080+
},
1081+
{
1082+
name: "polygon with hole partially outside rect removes hole",
1083+
input: "POLYGON((0 0,0 6,6 6,6 0,0 0),(1 1,3 1,3 3,1 3,1 1))",
1084+
rect: geom.NewEnvelope(geom.XY{X: 2, Y: 2}, geom.XY{X: 5, Y: 5}),
1085+
want: "POLYGON((2 3,2 5,5 5,5 2,3 2,3 3,2 3))",
1086+
},
1087+
{
1088+
name: "empty envelope",
1089+
input: "POLYGON((0 0,0 1,1 1,1 0,0 0))",
1090+
rect: geom.Envelope{},
1091+
want: "GEOMETRYCOLLECTION EMPTY",
1092+
},
1093+
} {
1094+
t.Run(tc.name, func(t *testing.T) {
1095+
got, err := geos.ClipByRect(geomFromWKT(t, tc.input), tc.rect)
1096+
skipIfUnsupported(t, err)
1097+
expectNoErr(t, err)
1098+
expectGeomEq(t, got, geomFromWKT(t, tc.want), geom.IgnoreOrder)
1099+
})
1100+
}
1101+
}

internal/rawgeos/entrypoints.go

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -44,6 +44,16 @@ GEOSGeometry *GEOSCoverageSimplifyVW_r(GEOSContextHandle_t handle, const GEOSGeo
4444
int GEOSCoverageIsValid_r(GEOSContextHandle_t handle, const GEOSGeometry* g, double gapWidth, GEOSGeometry** invalidEdges) { return 2; }
4545
#endif
4646
47+
#define CLIP_BY_RECT_MIN_VERSION "3.5.0"
48+
#define CLIP_BY_RECT_MISSING ( \
49+
GEOS_VERSION_MAJOR < 3 || \
50+
(GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR < 5) \
51+
)
52+
#if CLIP_BY_RECT_MISSING
53+
// This stub implementation always fails:
54+
GEOSGeometry *GEOSClipByRect_r(GEOSContextHandle_t handle, const GEOSGeometry* g, double xmin, double ymin, double xmax, double ymax) { return NULL; }
55+
#endif
56+
4757
#define CONCAVE_HULL_MIN_VERSION "3.11.0"
4858
#define CONCAVE_HULL_MISSING ( \
4959
GEOS_VERSION_MAJOR < 3 || \
@@ -441,6 +451,18 @@ func Envelope(g geom.Geometry) (geom.Geometry, error) {
441451
return result, wrap(err, "executing GEOSEnvelope_r")
442452
}
443453

454+
func ClipByRect(g geom.Geometry, xmin, ymin, xmax, ymax float64) (geom.Geometry, error) {
455+
if C.CLIP_BY_RECT_MISSING != 0 {
456+
return geom.Geometry{}, UnsupportedGEOSVersionError{
457+
C.CLIP_BY_RECT_MIN_VERSION, "ClipByRect",
458+
}
459+
}
460+
result, err := unaryOpG(g, func(ctx C.GEOSContextHandle_t, g *C.GEOSGeometry) *C.GEOSGeometry {
461+
return C.GEOSClipByRect_r(ctx, g, C.double(xmin), C.double(ymin), C.double(xmax), C.double(ymax))
462+
})
463+
return result, wrap(err, "executing GEOSClipByRect_r")
464+
}
465+
444466
func Area(g geom.Geometry) (float64, error) {
445467
result, err := unaryOpF(g, func(h C.GEOSContextHandle_t, g *C.GEOSGeometry, d *C.double) C.int {
446468
return C.GEOSArea_r(h, g, d)

0 commit comments

Comments
 (0)