Add ClipByRect2D with Sutherland-Hodgman and Liang-Barsky algorithms

CHANGELOG.md

@@ -5,6 +5,13 @@
 - Add `NewEnvelopeXY` constructor for building an `Envelope` from variadic x/y
   coordinate pairs, following the existing `New*XY` constructor pattern.
 
+- Add `ClipByRect2D` function that clips a geometry to a 2D axis-aligned
+  rectangle (defined by an `Envelope`). It uses the Sutherland-Hodgman
+  algorithm for polygons and the Liang-Barsky algorithm for line strings. The
+  result is always 2D: any Z or M values on the input are discarded, avoiding
+  the ambiguity of synthesising Z/M values at rectangle corners introduced by
+  the clip.
+
 ## v0.59.0
 
 2026-03-27

geom/alg_clip_by_rect.go

@@ -0,0 +1,114 @@
+package geom
+
+// ClipByRect2D clips a geometry to the 2D axis-aligned rectangle defined by
+// the given [Envelope], returning the portion of the geometry that lies within
+// the rectangle. It uses the Sutherland-Hodgman algorithm for polygon clipping
+// and the Liang-Barsky algorithm for line clipping.
+//
+// The result is always 2D ([DimXY]): any Z or M values on the input are
+// discarded. Linear interpolation of Z/M at clip-edge intersections is
+// reasonably well-defined, but values that would have to be synthesised at
+// rectangle corners (where the clipped boundary must traverse the rect) are
+// not, so this function avoids the problem by reducing to two dimensions
+// throughout.
+func ClipByRect2D(g Geometry, rect Envelope) Geometry {
+	return clipByRect2D(g.Force2D(), rect)
+}
+
+// clipByRect2D dispatches by geometry type. The input must already have been
+// reduced to [DimXY] by the caller.
+func clipByRect2D(g Geometry, rect Envelope) Geometry {
+	switch g.Type() {
+	case TypePoint:
+		return clipPointByRect(g.MustAsPoint(), rect).AsGeometry()
+	case TypeMultiPoint:
+		return clipMultiPointByRect(g.MustAsMultiPoint(), rect).AsGeometry()
+	case TypeLineString:
+		return clipLineStringByRect(g.MustAsLineString(), rect)
+	case TypeMultiLineString:
+		return clipMultiLineStringByRect(g.MustAsMultiLineString(), rect).AsGeometry()
+	case TypePolygon:
+		return clipPolygonByRect(g.MustAsPolygon(), rect)
+	case TypeMultiPolygon:
+		return clipMultiPolygonByRect(g.MustAsMultiPolygon(), rect).AsGeometry()
+	case TypeGeometryCollection:
+		return clipGeometryCollectionByRect(g.MustAsGeometryCollection(), rect).AsGeometry()
+	default:
+		panic("unknown geometry: " + g.Type().String())
+	}
+}
+
+func clipPointByRect(p Point, rect Envelope) Point {
+	xy, ok := p.XY()
+	if !ok {
+		return p
+	}
+	if rect.Contains(xy) {
+		return p
+	}
+	return Point{}
+}
+
+func clipMultiPointByRect(mp MultiPoint, rect Envelope) MultiPoint {
+	n := mp.NumPoints()
+	var pts []Point
+	for i := 0; i < n; i++ {
+		clipped := clipPointByRect(mp.PointN(i), rect)
+		if !clipped.IsEmpty() {
+			pts = append(pts, clipped)
+		}
+	}
+	return NewMultiPoint(pts)
+}
+
+func clipMultiLineStringByRect(mls MultiLineString, rect Envelope) MultiLineString {
+	n := mls.NumLineStrings()
+	var lines []LineString
+	for i := 0; i < n; i++ {
+		clipped := clipLineStringByRect(mls.LineStringN(i), rect)
+		switch clipped.Type() {
+		case TypeLineString:
+			ls := clipped.MustAsLineString()
+			if !ls.IsEmpty() {
+				lines = append(lines, ls)
+			}
+		case TypeMultiLineString:
+			lines = append(lines, clipped.MustAsMultiLineString().Dump()...)
+		default:
+			panic("unexpected type from clipLineStringByRect: " + clipped.Type().String())
+		}
+	}
+	return NewMultiLineString(lines)
+}
+
+func clipMultiPolygonByRect(mp MultiPolygon, rect Envelope) MultiPolygon {
+	n := mp.NumPolygons()
+	var polys []Polygon
+	for i := 0; i < n; i++ {
+		clipped := clipPolygonByRect(mp.PolygonN(i), rect)
+		switch clipped.Type() {
+		case TypePolygon:
+			p := clipped.MustAsPolygon()
+			if !p.IsEmpty() {
+				polys = append(polys, p)
+			}
+		case TypeMultiPolygon:
+			polys = append(polys, clipped.MustAsMultiPolygon().Dump()...)
+		default:
+			panic("unexpected type from clipPolygonByRect: " + clipped.Type().String())
+		}
+	}
+	return NewMultiPolygon(polys)
+}
+
+func clipGeometryCollectionByRect(gc GeometryCollection, rect Envelope) GeometryCollection {
+	n := gc.NumGeometries()
+	var geoms []Geometry
+	for i := 0; i < n; i++ {
+		clipped := clipByRect2D(gc.GeometryN(i), rect)
+		if !clipped.IsEmpty() {
+			geoms = append(geoms, clipped)
+		}
+	}
+	return NewGeometryCollection(geoms)
+}

geom/alg_clip_by_rect_liang_barsky.go

@@ -0,0 +1,98 @@
+package geom
+
+func clipLineStringByRect(ls LineString, rect Envelope) Geometry {
+	emptyLine := LineString{}.AsGeometry()
+
+	seq := ls.Coordinates()
+	n := seq.Length()
+
+	lo, hi, ok := rect.MinMaxXYs()
+	if !ok {
+		return emptyLine
+	}
+
+	var chains [][]XY
+	var cur []XY
+
+	for i := 0; i < n-1; i++ {
+		a := seq.GetXY(i)
+		b := seq.GetXY(i + 1)
+
+		tMin, tMax, ok := clipSegmentParams(a, b, lo, hi)
+		if !ok {
+			if len(cur) > 0 {
+				chains = append(chains, cur)
+				cur = nil
+			}
+			continue
+		}
+
+		ca := lerpXY(a, b, tMin)
+		cb := lerpXY(a, b, tMax)
+
+		if len(cur) > 0 && cur[len(cur)-1] == ca {
+			cur = append(cur, cb)
+		} else {
+			if len(cur) > 0 {
+				chains = append(chains, cur)
+			}
+			cur = []XY{ca, cb}
+		}
+	}
+	if len(cur) > 0 {
+		chains = append(chains, cur)
+	}
+
+	if len(chains) == 0 {
+		return emptyLine
+	}
+
+	lines := make([]LineString, len(chains))
+	for i, c := range chains {
+		lines[i] = NewLineString(xysToSeq(c))
+	}
+
+	if len(lines) == 1 {
+		return lines[0].AsGeometry()
+	}
+	return NewMultiLineString(lines).AsGeometry()
+}
+
+// clipSegmentParams uses the Liang-Barsky algorithm to compute the parametric
+// range [tMin, tMax] of segment a->b that lies inside the axis-aligned
+// rectangle from lo to hi. It returns false if no segment of positive length
+// survives.
+func clipSegmentParams(a, b, lo, hi XY) (float64, float64, bool) {
+	tMin := 0.0
+	tMax := 1.0
+	dx := b.X - a.X
+	dy := b.Y - a.Y
+	for _, pq := range [4][2]float64{
+		{-dx, a.X - lo.X}, // left
+		{dx, hi.X - a.X},  // right
+		{-dy, a.Y - lo.Y}, // bottom
+		{dy, hi.Y - a.Y},  // top
+	} {
+		p, q := pq[0], pq[1]
+		if p == 0 {
+			if q < 0 {
+				return 0, 0, false
+			}
+			continue
+		}
+		t := q / p
+		if p < 0 {
+			if t > tMin {
+				tMin = t
+			}
+		} else {
+			if t < tMax {
+				tMax = t
+			}
+		}
+		if tMin >= tMax {
+			return 0, 0, false
+		}
+	}
+	return tMin, tMax, true
+}