diff --git a/CHANGELOG.md b/CHANGELOG.md
index e8c49802..7ea502bf 100644
--- a/CHANGELOG.md
+++ b/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
diff --git a/geom/alg_clip_by_rect.go b/geom/alg_clip_by_rect.go
new file mode 100644
index 00000000..9f9715f6
--- /dev/null
+++ b/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)
+}
diff --git a/geom/alg_clip_by_rect_liang_barsky.go b/geom/alg_clip_by_rect_liang_barsky.go
new file mode 100644
index 00000000..8a54ede7
--- /dev/null
+++ b/geom/alg_clip_by_rect_liang_barsky.go
@@ -0,0 +1,138 @@
+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)
+
+		ca, cb, ok := clipSegment(a, b, lo, hi)
+		if !ok {
+			if len(cur) > 0 {
+				chains = append(chains, cur)
+				cur = nil
+			}
+			continue
+		}
+
+		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()
+}
+
+// clipSegment uses the Liang-Barsky algorithm to compute the portion of segment
+// a->b that lies inside the axis-aligned rectangle from lo to hi. It returns the
+// clipped endpoints ca and cb, and false if no segment of positive length
+// survives.
+//
+// An endpoint that a rect edge introduced (i.e. where the segment was cut by
+// that edge) is snapped exactly onto the edge's coordinate, rather than being
+// left a few ULPs off by interpolation. This matches the exact-boundary
+// guarantee that the polygon clipper provides, and keeps results stable under
+// this package's zero-tolerance equality.
+func clipSegment(a, b, lo, hi XY) (ca, cb XY, ok bool) {
+	tMin := 0.0
+	tMax := 1.0
+	dx := b.X - a.X
+	dy := b.Y - a.Y
+
+	// Each rect edge contributes a Liang-Barsky constraint.
+	edges := [4]clipEdge{
+		{-dx, a.X - lo.X, true, lo.X},  // left
+		{dx, hi.X - a.X, true, hi.X},   // right
+		{-dy, a.Y - lo.Y, false, lo.Y}, // bottom
+		{dy, hi.Y - a.Y, false, hi.Y},  // top
+	}
+
+	// minEdge/maxEdge record which edge last bound tMin/tMax, or -1 if the
+	// endpoint is still the original a/b (inside the rect) and must not be
+	// snapped.
+	minEdge, maxEdge := -1, -1
+	for i, e := range edges {
+		if e.p == 0 {
+			if e.q < 0 {
+				return XY{}, XY{}, false
+			}
+			continue
+		}
+		t := e.q / e.p
+		if e.p < 0 {
+			if t > tMin {
+				tMin = t
+				minEdge = i
+			}
+		} else {
+			if t < tMax {
+				tMax = t
+				maxEdge = i
+			}
+		}
+		if tMin >= tMax {
+			return XY{}, XY{}, false
+		}
+	}
+
+	ca = snapToEdge(lerpXY(a, b, tMin), edges, minEdge)
+	cb = snapToEdge(lerpXY(a, b, tMax), edges, maxEdge)
+	return ca, cb, true
+}
+
+// clipEdge describes one Liang-Barsky constraint, corresponding to a single rect
+// edge. p and q are the Liang-Barsky numerator/denominator for the edge. axisX
+// reports whether the edge fixes the X coordinate (left/right) or the Y
+// coordinate (bottom/top); val is the boundary coordinate to snap that axis to.
+type clipEdge struct {
+	p, q  float64
+	axisX bool
+	val   float64
+}
+
+// snapToEdge fixes the axis controlled by edge i exactly onto its boundary
+// coordinate, leaving the interpolated point xy unchanged when i is -1 (the
+// endpoint lies inside the rect and was not introduced by an edge).
+func snapToEdge(xy XY, edges [4]clipEdge, i int) XY {
+	if i < 0 {
+		return xy
+	}
+	if edges[i].axisX {
+		xy.X = edges[i].val
+	} else {
+		xy.Y = edges[i].val
+	}
+	return xy
+}
diff --git a/geom/alg_clip_by_rect_sutherland_hodgman.go b/geom/alg_clip_by_rect_sutherland_hodgman.go
new file mode 100644
index 00000000..816fdee1
--- /dev/null
+++ b/geom/alg_clip_by_rect_sutherland_hodgman.go
@@ -0,0 +1,587 @@
+package geom
+
+import (
+	"fmt"
+	"sort"
+)
+
+// polygonClipper holds precomputed values for clipping a polygon against an
+// axis-aligned rectangle.
+type polygonClipper struct {
+	lo, hi  XY
+	w, h    float64
+	perim   float64
+	corners [4]float64 // CCW corner parameters: BL, BR, TR, TL
+}
+
+func newPolygonClipper(lo, hi XY) polygonClipper {
+	w := hi.X - lo.X
+	h := hi.Y - lo.Y
+	return polygonClipper{
+		lo:      lo,
+		hi:      hi,
+		w:       w,
+		h:       h,
+		perim:   2*w + 2*h,
+		corners: [4]float64{0, w, w + h, 2*w + h},
+	}
+}
+
+func clipPolygonByRect(p Polygon, rect Envelope) Geometry {
+	emptyPoly := Polygon{}.AsGeometry()
+
+	if p.IsEmpty() {
+		return emptyPoly
+	}
+
+	// Degenerate rect: point or line envelope → empty polygon.
+	if !rect.IsRectangle() {
+		return emptyPoly
+	}
+
+	lo, hi, ok := rect.MinMaxXYs()
+	if !ok {
+		return emptyPoly
+	}
+
+	pEnv := p.Envelope()
+
+	// Fast path: envelopes disjoint.
+	if !rect.Intersects(pEnv) {
+		return emptyPoly
+	}
+
+	// Fast path: polygon entirely inside rect.
+	if rect.Covers(pEnv) {
+		return p.AsGeometry()
+	}
+
+	// Normalise to CCW exterior / CW holes. This ensures all clipped rings
+	// have known winding, which the topology resolution depends on.
+	p = p.ForceCCW()
+
+	c := newPolygonClipper(lo, hi)
+
+	// Clip exterior ring.
+	clippedExt := c.clipRingSH(p.ExteriorRing())
+	if len(clippedExt) == 0 {
+		return emptyPoly
+	}
+
+	// Classify holes.
+	var freeHoles []LineString
+	var clippedHoles [][]XY
+	for i := 0; i < p.NumInteriorRings(); i++ {
+		hole := p.InteriorRingN(i)
+		holeEnv := hole.Envelope()
+
+		if !rect.Intersects(holeEnv) {
+			// Hole entirely outside rect → discard.
+			continue
+		}
+
+		// A hole is "free" (entirely inside the rect with no boundary
+		// contact) only if its envelope is strictly inside the rect. If any
+		// vertex lies on the rect boundary, it must be clipped to avoid
+		// producing invalid polygons with shared edges.
+		if rect.Covers(holeEnv) && !c.ringTouchesRectBoundary(hole) {
+			freeHoles = append(freeHoles, hole)
+			continue
+		}
+
+		// Hole touches or crosses rect boundary — clip it.
+		clippedHole := c.clipRingSH(hole)
+		if len(clippedHole) > 0 {
+			clippedHoles = append(clippedHoles, clippedHole)
+		}
+	}
+
+	return c.resolveClippedPolygon(clippedExt, freeHoles, clippedHoles)
+}
+
+// ringTouchesRectBoundary reports whether any vertex of the ring lies on the rect
+// boundary.
+func (c *polygonClipper) ringTouchesRectBoundary(ring LineString) bool {
+	seq := ring.Coordinates()
+	for i := 0; i < seq.Length(); i++ {
+		xy := seq.GetXY(i)
+		if xy.X == c.lo.X || xy.X == c.hi.X || xy.Y == c.lo.Y || xy.Y == c.hi.Y {
+			return true
+		}
+	}
+	return false
+}
+
+// resolveClippedPolygon takes the clipped exterior ring, classified holes, and
+// produces the output Geometry (Polygon or MultiPolygon). The input polygon
+// must have been normalised to CCW (exterior CCW, holes CW) before clipping,
+// so that all clipped rings have known winding.
+func (c *polygonClipper) resolveClippedPolygon(
+	clippedExterior []XY,
+	freeHoles []LineString,
+	clippedHoles [][]XY,
+) Geometry {
+	extArcs := c.extractInteriorArcs(clippedExterior)
+
+	// Collect all arcs.
+	var allArcs []interiorArc
+	allArcs = append(allArcs, extArcs...)
+
+	emptyPoly := Polygon{}.AsGeometry()
+	for _, hole := range clippedHoles {
+		arcs := c.extractInteriorArcs(hole)
+		if len(arcs) == 0 {
+			// Clipped hole is entirely on the rect boundary — it covers
+			// the entire rect. No polygon area survives.
+			return emptyPoly
+		}
+		allArcs = append(allArcs, arcs...)
+	}
+
+	var outputRings [][]XY
+	if len(allArcs) == 0 {
+		// No interior arcs: the polygon contains the entire rect.
+		// Use the clipped exterior directly — it is the rect itself.
+		outputRings = append(outputRings, clippedExterior)
+	} else {
+		outputRings = c.walkArcs(allArcs)
+	}
+
+	// Build output polygons. Each output ring is an exterior (already closed).
+	var polys []Polygon
+	for _, ring := range outputRings {
+		exterior := NewLineString(xysToSeq(ring))
+		polys = append(polys, NewPolygon([]LineString{exterior}))
+	}
+
+	// Assign free holes to the correct exterior ring.
+	for _, hole := range freeHoles {
+		holeXY, ok := hole.StartPoint().XY()
+		if !ok {
+			continue
+		}
+		for i, ring := range outputRings {
+			if c.pointInRing(holeXY, ring) {
+				existingRings := polys[i].DumpRings()
+				existingRings = append(existingRings, hole)
+				polys[i] = NewPolygon(existingRings)
+				break
+			}
+		}
+	}
+
+	// Return Polygon or MultiPolygon.
+	if len(polys) == 0 {
+		return emptyPoly
+	}
+	if len(polys) == 1 {
+		return polys[0].AsGeometry()
+	}
+	return NewMultiPolygon(polys).AsGeometry()
+}
+
+// walkArcs performs the topology resolution walk, producing output rings from
+// the collected interior arcs. It pairs arc endpoints along the rect boundary
+// (CCW) and traces complete rings.
+func (c *polygonClipper) walkArcs(arcs []interiorArc) [][]XY {
+	// Build a sorted list of arc "end" events and a map from "end param" to
+	// the index of the arc that ends there. Also build a map from
+	// "start param" to arc index.
+	type endpoint struct {
+		param float64
+		idx   int
+	}
+	var endPoints []endpoint // sorted by param
+	startByParam := make(map[float64]int)
+	for i, a := range arcs {
+		endPoints = append(endPoints, endpoint{a.endParam, i})
+		startByParam[a.startParam] = i
+	}
+
+	// Sort endpoints by param.
+	sort.Slice(endPoints, func(i, j int) bool {
+		return endPoints[i].param < endPoints[j].param
+	})
+
+	// Collect all start params sorted.
+	var startParams []float64
+	for p := range startByParam {
+		startParams = append(startParams, p)
+	}
+	sort.Float64s(startParams)
+
+	// For a given end param, find the next start param going CCW. The ok
+	// result is false only when no successor can be found, which requires
+	// non-finite boundary params (e.g. NaN coordinates in the input): the
+	// preference loop rejects every param and the coincident fallback can't
+	// match either. The caller closes the current ring as-is in that case,
+	// producing a nonsensical-but-bounded result rather than looping forever.
+	findNextStart := func(endParam float64) (float64, int, bool) {
+		// Prefer the nearest start param strictly ahead of endParam going CCW.
+		best := -1.0
+		bestDist := c.perim + 1
+		for _, sp := range startParams {
+			d := c.ccwDist(endParam, sp)
+			if d > 0 && d < bestDist {
+				bestDist = d
+				best = sp
+			}
+		}
+		if best >= 0 {
+			return best, startByParam[best], true
+		}
+		// No start param is strictly ahead: fall back to one coincident with
+		// endParam, where the next arc begins exactly where this one ended.
+		for _, sp := range startParams {
+			if sp == endParam {
+				return sp, startByParam[sp], true
+			}
+		}
+		return 0, 0, false
+	}
+
+	used := make([]bool, len(arcs))
+	var rings [][]XY
+
+	for {
+		// Find first unused arc.
+		firstIdx := -1
+		for i := range arcs {
+			if !used[i] {
+				firstIdx = i
+				break
+			}
+		}
+		if firstIdx < 0 {
+			break
+		}
+
+		var ring []XY
+		curIdx := firstIdx
+
+		for {
+			used[curIdx] = true
+			a := arcs[curIdx]
+
+			// Append arc coordinates.
+			ring = append(ring, a.coords...)
+
+			// Find next arc via boundary.
+			nextStartParam, nextIdx, ok := findNextStart(a.endParam)
+			if !ok {
+				break // No successor arc; close the ring as-is.
+			}
+
+			// Build boundary path from this arc's end to the next arc's start.
+			ring = c.appendBoundaryPath(ring, a.endParam, nextStartParam)
+
+			if nextIdx == firstIdx {
+				break // Ring complete.
+			}
+			curIdx = nextIdx
+		}
+
+		// Close the ring.
+		ring = append(ring, ring[0])
+		ring = removeDupConsecutive(ring)
+		if len(ring) >= 4 {
+			rings = append(rings, ring)
+		}
+	}
+
+	return rings
+}
+
+// appendBoundaryPath appends coordinates along the rect boundary going CCW
+// from startParam to endParam. It includes rect corners between the two
+// parameters but does NOT include the start or end points themselves (those are
+// the last/first coordinates of the adjacent arcs).
+func (c *polygonClipper) appendBoundaryPath(dst []XY, startParam, endParam float64) []XY {
+	// Find the first corner after startParam going CCW. The corners are in
+	// ascending parameter order, so this is the first with cp > startParam.
+	// If startParam is past the last corner (on the left edge), no corner
+	// qualifies and firstIdx stays at 0 — the bottom-left corner is the
+	// next one going CCW.
+	firstIdx := 0
+	for i, cp := range c.corners {
+		if cp > startParam {
+			firstIdx = i
+			break
+		}
+	}
+
+	// Iterate corners in CCW order from firstIdx, collecting those strictly
+	// between startParam and endParam.
+	totalDist := c.ccwDist(startParam, endParam)
+	for k := 0; k < 4; k++ {
+		cp := c.corners[(firstIdx+k)%4]
+		d := c.ccwDist(startParam, cp)
+		if d >= totalDist {
+			break
+		}
+		dst = append(dst, c.paramToXY(cp))
+	}
+	return dst
+}
+
+// extractInteriorArcs decomposes a clipped ring into interior arcs. The ring
+// must be explicitly closed (first == last). Each arc starts and ends at a
+// point on the rect boundary. Returns an empty slice if the ring has no
+// interior arcs (entirely on the boundary or entirely interior with no
+// boundary contact).
+func (c *polygonClipper) extractInteriorArcs(ring []XY) []interiorArc {
+	n := len(ring)
+	if n < 4 {
+		return nil
+	}
+
+	// There are n-1 edges in a closed ring (the last vertex == first vertex,
+	// so edge n-1 from ring[n-1] to ring[0] is degenerate and skipped).
+	numEdges := n - 1
+
+	// Classify each edge as boundary (both endpoints on same rect edge) or interior.
+	isBdry := make([]bool, numEdges)
+	for i := 0; i < numEdges; i++ {
+		isBdry[i] = c.isSameRectEdge(ring[i], ring[i+1])
+	}
+
+	// Find a starting boundary edge so we can walk from there. If there are
+	// no boundary edges, the entire ring is interior (no boundary contact).
+	start := -1
+	for i := 0; i < numEdges; i++ {
+		if isBdry[i] {
+			start = i
+			break
+		}
+	}
+	if start < 0 {
+		return nil
+	}
+
+	// Walk around the ring starting from a boundary edge.
+	var arcs []interiorArc
+	i := start
+	for steps := 0; steps < numEdges; {
+		// Skip boundary edges.
+		for steps < numEdges && isBdry[i%numEdges] {
+			i++
+			steps++
+		}
+		if steps >= numEdges {
+			break
+		}
+		// Start of an interior arc at ring[i%numEdges].
+		var arcCoords []XY
+		arcCoords = append(arcCoords, ring[i%numEdges])
+		i++
+		steps++
+		for steps < numEdges && !isBdry[i%numEdges] {
+			arcCoords = append(arcCoords, ring[i%numEdges])
+			i++
+			steps++
+		}
+		if steps < numEdges {
+			// The arc ends at ring[i%numEdges] (the start of the next boundary edge).
+			arcCoords = append(arcCoords, ring[i%numEdges])
+		} else {
+			// Wrapped around; the arc ends at ring[start] (where we began).
+			arcCoords = append(arcCoords, ring[start])
+		}
+
+		sp := c.rectBoundaryParam(arcCoords[0])
+		ep := c.rectBoundaryParam(arcCoords[len(arcCoords)-1])
+		arcs = append(arcs, interiorArc{
+			coords:     arcCoords,
+			startParam: sp,
+			endParam:   ep,
+		})
+	}
+	return arcs
+}
+
+// interiorArc represents a portion of a clipped ring that passes through the
+// interior of the clipping rectangle (not along its boundary). The first and
+// last coordinates are on the rectangle boundary; everything in between is in
+// the interior.
+type interiorArc struct {
+	coords     []XY
+	startParam float64 // boundary parameter of first point
+	endParam   float64 // boundary parameter of last point
+}
+
+// clipRingSH clips a closed ring against an axis-aligned rectangle using the
+// Sutherland-Hodgman algorithm. It returns the clipped ring as a closed []XY
+// (first point == last point), or nil if the ring is entirely outside the
+// rectangle. The input [LineString] must be a closed ring.
+func (c *polygonClipper) clipRingSH(ring LineString) []XY {
+	coords := ring.Coordinates().asXYs()
+	if len(coords) < 4 {
+		return nil
+	}
+
+	// Clip against each of the 4 edges.
+	coords = c.clipToEdge(coords,
+		func(xy XY) bool { return xy.X >= c.lo.X },
+		func(dst, src []XY, ai, bi int) []XY { return appendInterpX(dst, src, ai, bi, c.lo.X) })
+	coords = c.clipToEdge(coords,
+		func(xy XY) bool { return xy.X <= c.hi.X },
+		func(dst, src []XY, ai, bi int) []XY { return appendInterpX(dst, src, ai, bi, c.hi.X) })
+	coords = c.clipToEdge(coords,
+		func(xy XY) bool { return xy.Y >= c.lo.Y },
+		func(dst, src []XY, ai, bi int) []XY { return appendInterpY(dst, src, ai, bi, c.lo.Y) })
+	coords = c.clipToEdge(coords,
+		func(xy XY) bool { return xy.Y <= c.hi.Y },
+		func(dst, src []XY, ai, bi int) []XY { return appendInterpY(dst, src, ai, bi, c.hi.Y) })
+
+	coords = removeDupConsecutive(coords)
+	if len(coords) < 4 {
+		return nil
+	}
+	return coords
+}
+
+// clipToEdge performs one pass of the Sutherland-Hodgman algorithm, clipping a
+// closed ring against a single half-plane defined by isInside. The
+// appendIntersection function appends the intersection point between points at
+// indices ai and bi to dst. The output is also an explicitly closed ring.
+func (c *polygonClipper) clipToEdge(
+	coords []XY,
+	isInside func(xy XY) bool,
+	appendIntersection func(dst, coords []XY, ai, bi int) []XY,
+) []XY {
+	n := len(coords)
+	if n == 0 {
+		return nil
+	}
+	var out []XY
+	ai := n - 1
+	for bi := 0; bi < n; bi++ {
+		aIn := isInside(coords[ai])
+		bIn := isInside(coords[bi])
+		switch {
+		case aIn && bIn:
+			out = append(out, coords[bi])
+		case aIn && !bIn:
+			out = appendIntersection(out, coords, ai, bi)
+		case !aIn && bIn:
+			out = appendIntersection(out, coords, ai, bi)
+			out = append(out, coords[bi])
+		}
+		ai = bi
+	}
+	// Ensure the output is explicitly closed. When the input's closing vertex
+	// is outside the clip region, the degenerate closing edge emits nothing,
+	// leaving the output open.
+	if len(out) > 0 && out[0] != out[len(out)-1] {
+		out = append(out, out[0])
+	}
+	return out
+}
+
+// appendInterpX appends the intersection of the segment between points ai and
+// bi with the vertical line x=k. The X coordinate is set exactly; Y is
+// linearly interpolated.
+func appendInterpX(dst, coords []XY, ai, bi int, x float64) []XY {
+	a := coords[ai]
+	b := coords[bi]
+	t := (x - a.X) / (b.X - a.X)
+	return append(dst, XY{X: x, Y: lerp(a.Y, b.Y, t)})
+}
+
+// appendInterpY appends the intersection of the segment between points ai and
+// bi with the horizontal line y=k. The Y coordinate is set exactly; X is
+// linearly interpolated.
+func appendInterpY(dst, coords []XY, ai, bi int, y float64) []XY {
+	a := coords[ai]
+	b := coords[bi]
+	t := (y - a.Y) / (b.Y - a.Y)
+	return append(dst, XY{X: lerp(a.X, b.X, t), Y: y})
+}
+
+// rectBoundaryParam returns the CCW boundary parameter for a point on the rect
+// boundary. The parameterisation starts at the bottom-left corner and goes
+// counter-clockwise: bottom→right→top→left.
+func (c *polygonClipper) rectBoundaryParam(xy XY) float64 {
+	switch {
+	case xy.Y == c.lo.Y: // bottom edge
+		return xy.X - c.lo.X
+	case xy.X == c.hi.X: // right edge
+		return c.w + (xy.Y - c.lo.Y)
+	case xy.Y == c.hi.Y: // top edge
+		return c.w + c.h + (c.hi.X - xy.X)
+	case xy.X == c.lo.X: // left edge
+		return 2*c.w + c.h + (c.hi.Y - xy.Y)
+	default:
+		panic(fmt.Sprintf("point %v not on rect boundary [%v, %v]", xy, c.lo, c.hi))
+	}
+}
+
+// paramToXY converts a CCW boundary parameter back to an XY coordinate.
+func (c *polygonClipper) paramToXY(param float64) XY {
+	switch {
+	case param < c.w: // bottom edge
+		return XY{X: c.lo.X + param, Y: c.lo.Y}
+	case param < c.w+c.h: // right edge
+		return XY{X: c.hi.X, Y: c.lo.Y + param - c.w}
+	case param < 2*c.w+c.h: // top edge
+		return XY{X: c.hi.X - (param - c.w - c.h), Y: c.hi.Y}
+	case param < c.perim: // left edge
+		return XY{X: c.lo.X, Y: c.hi.Y - (param - 2*c.w - c.h)}
+	default:
+		panic(fmt.Sprintf("boundary parameter %v out of range [0, %v)", param, c.perim))
+	}
+}
+
+// ccwDist returns the CCW distance from param a to param b on the rect
+// boundary.
+func (c *polygonClipper) ccwDist(a, b float64) float64 {
+	d := b - a
+	if d < 0 {
+		d += c.perim
+	}
+	return d
+}
+
+// isSameRectEdge returns true if both points lie on the same edge of the
+// rectangle.
+func (c *polygonClipper) isSameRectEdge(a, b XY) bool {
+	return (a.X == c.lo.X && b.X == c.lo.X) ||
+		(a.X == c.hi.X && b.X == c.hi.X) ||
+		(a.Y == c.lo.Y && b.Y == c.lo.Y) ||
+		(a.Y == c.hi.Y && b.Y == c.hi.Y)
+}
+
+// pointInRing returns true if xy is inside the ring defined by the given
+// closed []XY (first point == last point), using the ray-casting algorithm.
+func (c *polygonClipper) pointInRing(xy XY, ring []XY) bool {
+	n := len(ring)
+	inside := false
+	for i := 0; i < n; i++ {
+		j := (i + 1) % n
+		yi, yj := ring[i].Y, ring[j].Y
+		xi, xj := ring[i].X, ring[j].X
+		if (yi > xy.Y) != (yj > xy.Y) {
+			slope := (xy.Y - yi) / (yj - yi)
+			xIntersect := xi + slope*(xj-xi)
+			if xy.X < xIntersect {
+				inside = !inside
+			}
+		}
+	}
+	return inside
+}
+
+// removeDupConsecutive removes consecutive points with identical X,Y.
+// For closed rings (first == last), the closing point is preserved.
+func removeDupConsecutive(coords []XY) []XY {
+	if len(coords) == 0 {
+		return nil
+	}
+	out := coords[:1]
+	for i := 1; i < len(coords); i++ {
+		if coords[i] != out[len(out)-1] {
+			out = append(out, coords[i])
+		}
+	}
+	return out
+}
diff --git a/geom/alg_clip_by_rect_sutherland_hodgman_test.go b/geom/alg_clip_by_rect_sutherland_hodgman_test.go
new file mode 100644
index 00000000..d7930c6c
--- /dev/null
+++ b/geom/alg_clip_by_rect_sutherland_hodgman_test.go
@@ -0,0 +1,570 @@
+package geom_test
+
+import (
+	"testing"
+
+	"github.com/peterstace/simplefeatures/geom"
+	"github.com/peterstace/simplefeatures/internal/test"
+)
+
+var d4Transforms = []struct {
+	name string
+	fn   func(geom.XY) geom.XY
+}{
+	{"identity", func(xy geom.XY) geom.XY { return xy }},
+	{"rot90", func(xy geom.XY) geom.XY { return geom.XY{X: -xy.Y, Y: xy.X} }},
+	{"rot180", func(xy geom.XY) geom.XY { return geom.XY{X: -xy.X, Y: -xy.Y} }},
+	{"rot270", func(xy geom.XY) geom.XY { return geom.XY{X: xy.Y, Y: -xy.X} }},
+	{"reflect_x", func(xy geom.XY) geom.XY { return geom.XY{X: -xy.X, Y: xy.Y} }},
+	{"reflect_y", func(xy geom.XY) geom.XY { return geom.XY{X: xy.X, Y: -xy.Y} }},
+	{"reflect_diag", func(xy geom.XY) geom.XY { return geom.XY{X: xy.Y, Y: xy.X} }},
+	{"reflect_anti", func(xy geom.XY) geom.XY { return geom.XY{X: -xy.Y, Y: -xy.X} }},
+}
+
+func TestClipByRect2D(t *testing.T) {
+	// R is a non-square rectangle so that the D4 transforms produce distinct
+	// configurations.
+	rect := geom.NewEnvelope(geom.XY{X: 1, Y: 2}, geom.XY{X: 5, Y: 4})
+
+	for _, tt := range []struct {
+		name  string
+		input string
+		want  string
+		opts  []geom.ExactEqualsOption
+	}{
+		// PT1: Empty Point
+		{"PT1", "POINT EMPTY", "POINT EMPTY", nil},
+		// PT2: Point strictly inside R
+		{"PT2", "POINT(3 3)", "POINT(3 3)", nil},
+		// PT3: Point strictly outside R
+		{"PT3", "POINT(0 0)", "POINT EMPTY", nil},
+		// PT4: Point on edge of R (on left edge, x=1)
+		{"PT4", "POINT(1 3)", "POINT(1 3)", nil},
+		// PT5: Point on corner of R (bottom-left corner)
+		{"PT5", "POINT(1 2)", "POINT(1 2)", nil},
+
+		// MP1: Empty MultiPoint
+		{"MP1", "MULTIPOINT EMPTY", "MULTIPOINT EMPTY", nil},
+		// MP2: All points inside R
+		{"MP2", "MULTIPOINT(2 3,3 3)", "MULTIPOINT(2 3,3 3)", nil},
+		// MP3: All points outside R
+		{"MP3", "MULTIPOINT(0 0,6 6)", "MULTIPOINT EMPTY", nil},
+		// MP4: Some points inside, some outside
+		{"MP4", "MULTIPOINT(3 3,0 0)", "MULTIPOINT(3 3)", nil},
+		// MP5: Single point inside R
+		{"MP5", "MULTIPOINT(3 3)", "MULTIPOINT(3 3)", nil},
+		// MP6: Points on edges and corners of R
+		{"MP6", "MULTIPOINT(1 3,1 2)", "MULTIPOINT(1 3,1 2)", nil},
+		// MP7: Mix of inside, on-boundary, and outside
+		{"MP7", "MULTIPOINT(3 3,1 3,0 0)", "MULTIPOINT(3 3,1 3)", nil},
+		// MP8: MultiPoint containing empty points
+		{"MP8", "MULTIPOINT(3 3,EMPTY)", "MULTIPOINT(3 3)", nil},
+		// MP9: XYZ MultiPoint, all points outside R — output is XY EMPTY
+		{"MP9", "MULTIPOINT Z(0 0 7,6 6 8)", "MULTIPOINT EMPTY", nil},
+
+		// LS1: Empty LineString
+		{"LS1", "LINESTRING EMPTY", "LINESTRING EMPTY", nil},
+		// LS2: Entirely inside R
+		{"LS2", "LINESTRING(2 3,4 3)", "LINESTRING(2 3,4 3)", nil},
+		// LS3: Entirely outside R (no overlap)
+		{"LS3", "LINESTRING(6 5,7 6)", "LINESTRING EMPTY", nil},
+		// LS4: Entirely outside R on one side (all left of R)
+		{"LS4", "LINESTRING(0 3,0 3.5)", "LINESTRING EMPTY", nil},
+		// LS5: Crosses R, entering and exiting once
+		{"LS5", "LINESTRING(0 3,2 3,6 5)", "LINESTRING(1 3,2 3,4 4)", nil},
+		// LS6: Crosses R, entering and exiting multiple times
+		{"LS6", "LINESTRING(0 3,3 3,3 5,4 5,4 3,6 3)", "MULTILINESTRING((1 3,3 3,3 4),(4 4,4 3,5 3))", nil},
+		// LS7: One endpoint inside, one outside
+		{"LS7", "LINESTRING(3 3,6 3)", "LINESTRING(3 3,5 3)", nil},
+		// LS8: One endpoint outside, one inside
+		{"LS8", "LINESTRING(0 3,3 3)", "LINESTRING(1 3,3 3)", nil},
+		// LS9: Both endpoints outside, segment passes through R
+		{"LS9", "LINESTRING(0 3,6 3)", "LINESTRING(1 3,5 3)", nil},
+		// LS10: Endpoint exactly on edge of R, other inside
+		{"LS10", "LINESTRING(1 3,3 3)", "LINESTRING(1 3,3 3)", nil},
+		// LS11: Endpoint exactly on corner of R, other inside
+		{"LS11", "LINESTRING(1 2,3 3)", "LINESTRING(1 2,3 3)", nil},
+		// LS12: Both endpoints on boundary of R
+		{"LS12", "LINESTRING(1 3,5 3)", "LINESTRING(1 3,5 3)", nil},
+		// LS13: LineString lies entirely along one edge of R
+		{"LS13", "LINESTRING(2 2,4 2)", "LINESTRING(2 2,4 2)", nil},
+		// LS14: LineString lies entirely along two adjacent edges (L-shaped)
+		{"LS14", "LINESTRING(1 3,1 2,3 2)", "LINESTRING(1 3,1 2,3 2)", nil},
+		// LS15: Segment touches corner of R but does not enter (V-shape)
+		{"LS15", "LINESTRING(0 1,1 2,0 3)", "LINESTRING EMPTY", nil},
+		// LS16: Segment touches edge of R tangentially
+		{"LS16", "LINESTRING(2 5,3 4,4 5)", "LINESTRING EMPTY", nil},
+		// LS17: Diagonal line crossing two edges of R
+		{"LS17", "LINESTRING(0 0,6 6)", "LINESTRING(2 2,4 4)", nil},
+		// LS18: Axis-aligned line crossing two opposite edges of R
+		{"LS18", "LINESTRING(3 0,3 6)", "LINESTRING(3 2,3 4)", nil},
+		// LS19: Closed LineString (ring) inside R
+		{"LS19", "LINESTRING(2 2.5,4 2.5,3 3.5,2 2.5)", "LINESTRING(2 2.5,4 2.5,3 3.5,2 2.5)", nil},
+		// LS20: Closed LineString (ring) partially overlapping R
+		{"LS20", "LINESTRING(2 1,4 1,4 5,2 5,2 1)", "MULTILINESTRING((4 2,4 4),(2 4,2 2))", nil},
+		// LS21: Multi-segment LineString with some segments inside, some outside
+		{"LS21", "LINESTRING(3 3,4 3,6 3)", "LINESTRING(3 3,4 3,5 3)", nil},
+		// LS22: Zigzag LineString entering and exiting R many times
+		{"LS22", "LINESTRING(0 3,2 3,2 5,3 5,3 3,4 3,4 5,5 5,5 3,7 3)", "MULTILINESTRING((1 3,2 3,2 4),(3 4,3 3,4 3,4 4),(5 4,5 3))", nil},
+		// LS23: Vertex exactly on R boundary, adjacent edges inside
+		{"LS23", "LINESTRING(2 3,3 2,4 3)", "LINESTRING(2 3,3 2,4 3)", nil},
+		// LS24: Vertex exactly on R boundary, adjacent edges outside
+		{"LS24", "LINESTRING(0 1,1 3,0 5)", "LINESTRING EMPTY", nil},
+		// LS25: LineString passes through two opposite corners of R
+		{"LS25", "LINESTRING(0 1.5,6 4.5)", "LINESTRING(1 2,5 4)", nil},
+
+		// MLS1: Empty MultiLineString
+		{"MLS1", "MULTILINESTRING EMPTY", "MULTILINESTRING EMPTY", nil},
+		// MLS2: All component LineStrings inside R
+		{"MLS2", "MULTILINESTRING((2 3,4 3),(3 2.5,3 3.5))", "MULTILINESTRING((2 3,4 3),(3 2.5,3 3.5))", nil},
+		// MLS3: All component LineStrings outside R
+		{"MLS3", "MULTILINESTRING((6 5,7 6),(0 0,0 1))", "MULTILINESTRING EMPTY", nil},
+		// MLS4: Some components inside, some outside
+		{"MLS4", "MULTILINESTRING((2 3,4 3),(6 5,7 6))", "MULTILINESTRING((2 3,4 3))", nil},
+		// MLS5: Single component, clipped to a single segment
+		{"MLS5", "MULTILINESTRING((0 3,6 3))", "MULTILINESTRING((1 3,5 3))", nil},
+		// MLS6: One component crosses R, another is inside R
+		{"MLS6", "MULTILINESTRING((0 3,6 3),(3 2.5,3 3.5))", "MULTILINESTRING((1 3,5 3),(3 2.5,3 3.5))", nil},
+		// MLS7: Component that produces multiple fragments when clipped
+		{"MLS7", "MULTILINESTRING((0 3,3 3,3 5,4 5,4 3,6 3))", "MULTILINESTRING((1 3,3 3,3 4),(4 4,4 3,5 3))", nil},
+		// MLS8: MultiLineString containing empty LineStrings
+		{"MLS8", "MULTILINESTRING((2 3,4 3),EMPTY)", "MULTILINESTRING((2 3,4 3))", nil},
+		// MLS9: XYZ MultiLineString, all components outside R — output is XY EMPTY
+		{"MLS9", "MULTILINESTRING Z((6 5 1,7 6 2),(0 0 3,0 1 4))", "MULTILINESTRING EMPTY", nil},
+
+		// PG1: Empty Polygon
+		{"PG1", "POLYGON EMPTY", "POLYGON EMPTY", nil},
+		// PG2: Entirely inside R
+		{"PG2", "POLYGON((2 2.5,4 2.5,4 3.5,2 3.5,2 2.5))", "POLYGON((2 2.5,4 2.5,4 3.5,2 3.5,2 2.5))", []geom.ExactEqualsOption{geom.IgnoreOrder}},
+		// PG3: Entirely outside R (no overlap)
+		{"PG3", "POLYGON((6 5,8 5,8 7,6 7,6 5))", "POLYGON EMPTY", nil},
+		// PG4: R entirely inside Polygon
+		{"PG4", "POLYGON((0 0,6 0,6 6,0 6,0 0))", "POLYGON((1 2,5 2,5 4,1 4,1 2))", []geom.ExactEqualsOption{geom.IgnoreOrder}},
+		// PG5: Partially overlapping one edge (left) of R
+		{"PG5", "POLYGON((0 2.5,3 2.5,3 3.5,0 3.5,0 2.5))", "POLYGON((1 2.5,3 2.5,3 3.5,1 3.5,1 2.5))", []geom.ExactEqualsOption{geom.IgnoreOrder}},
+		// PG6: Partially overlapping two adjacent edges (bottom-left corner clip)
+		{"PG6", "POLYGON((0 1,3 1,3 3,0 3,0 1))", "POLYGON((1 2,3 2,3 3,1 3,1 2))", []geom.ExactEqualsOption{geom.IgnoreOrder}},
+		// PG7: Partially overlapping two opposite edges (strip left-right)
+		{"PG7", "POLYGON((0 2.5,6 2.5,6 3.5,0 3.5,0 2.5))", "POLYGON((1 2.5,5 2.5,5 3.5,1 3.5,1 2.5))", []geom.ExactEqualsOption{geom.IgnoreOrder}},
+		// PG8: Partially overlapping three edges (left, bottom, right)
+		{"PG8", "POLYGON((0 1,6 1,6 3,0 3,0 1))", "POLYGON((1 2,5 2,5 3,1 3,1 2))", []geom.ExactEqualsOption{geom.IgnoreOrder}},
+		// PG9: Polygon shares entire bottom edge with R
+		{"PG9", "POLYGON((1 2,5 2,5 3,1 3,1 2))", "POLYGON((1 2,5 2,5 3,1 3,1 2))", []geom.ExactEqualsOption{geom.IgnoreOrder}},
+		// PG10: Polygon vertex exactly on R left edge (x=1)
+		{"PG10", "POLYGON((1 3,2 2.5,2 3.5,1 3))", "POLYGON((1 3,2 2.5,2 3.5,1 3))", []geom.ExactEqualsOption{geom.IgnoreOrder}},
+		// PG11: Polygon vertex exactly on R corner (1,2)
+		{"PG11", "POLYGON((1 2,3 2.5,3 3.5,1 2))", "POLYGON((1 2,3 2.5,3 3.5,1 2))", []geom.ExactEqualsOption{geom.IgnoreOrder}},
+		// PG12: Polygon edge collinear with R bottom edge, polygon inside R
+		{"PG12", "POLYGON((2 2,4 2,3 3,2 2))", "POLYGON((2 2,4 2,3 3,2 2))", []geom.ExactEqualsOption{geom.IgnoreOrder}},
+		// PG13: Polygon edge collinear with R bottom edge, polygon outside R
+		{"PG13", "POLYGON((2 2,4 2,3 1,2 2))", "POLYGON EMPTY", nil},
+		// PG14: Polygon touches R at single point on left edge (vertex-to-edge)
+		{"PG14", "POLYGON((0 2.5,1 3,0 3.5,0 2.5))", "POLYGON EMPTY", nil},
+		// PG15: Polygon touches R at single corner point (1,2)
+		{"PG15", "POLYGON((0 1,1 2,0 3,0 1))", "POLYGON EMPTY", nil},
+		// PG16: Convex polygon (large square) clipped at all 4 edges
+		{"PG16", "POLYGON((-1 0,7 0,7 6,-1 6,-1 0))", "POLYGON((1 2,5 2,5 4,1 4,1 2))", []geom.ExactEqualsOption{geom.IgnoreOrder}},
+		// PG17: Concave polygon, concavity inside R
+		{"PG17", "POLYGON((2 2.5,4 2.5,4 3.5,3 3,2 3.5,2 2.5))", "POLYGON((2 2.5,4 2.5,4 3.5,3 3,2 3.5,2 2.5))", []geom.ExactEqualsOption{geom.IgnoreOrder}},
+		// PG18: Concave polygon, concavity facing left R boundary
+		{"PG18", "POLYGON((0 2.5,3 2.5,3 3,2 3,2 3.5,0 3.5,0 2.5))", "POLYGON((1 2.5,3 2.5,3 3,2 3,2 3.5,1 3.5,1 2.5))", []geom.ExactEqualsOption{geom.IgnoreOrder}},
+		// PG19: C-shaped polygon, concavity crosses left edge → MultiPolygon
+		{
+			"PG19",
+			"POLYGON((0 2.5,3 2.5,3 2.8,0.5 2.8,0.5 3.2,3 3.2,3 3.5,0 3.5,0 2.5))",
+			"MULTIPOLYGON(((1 2.5,3 2.5,3 2.8,1 2.8,1 2.5)),((1 3.2,3 3.2,3 3.5,1 3.5,1 3.2)))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// PG20: Very thin sliver polygon partially inside R
+		{"PG20", "POLYGON((0 2.999,6 2.999,6 3.001,0 3.001,0 2.999))", "POLYGON((1 2.999,5 2.999,5 3.001,1 3.001,1 2.999))", []geom.ExactEqualsOption{geom.IgnoreOrder}},
+		// PG21: Triangle clipped at all 4 rect edges producing a polygon
+		{"PG21", "POLYGON((3 0,7 3,3 6,3 0))", "POLYGON((3 4,3 2,5 2,5 4,3 4))", []geom.ExactEqualsOption{geom.IgnoreOrder}},
+		// PG22: Counter-clockwise exterior ring (explicitly CCW, same as PG2 but winding verified by D4 reflections)
+		{"PG22", "POLYGON((2 2.5,2 3.5,4 3.5,4 2.5,2 2.5))", "POLYGON((2 2.5,2 3.5,4 3.5,4 2.5,2 2.5))", []geom.ExactEqualsOption{geom.IgnoreOrder}},
+
+		// PH1: Exterior and hole both inside R
+		{
+			"PH1",
+			"POLYGON((2 2.5,4 2.5,4 3.5,2 3.5,2 2.5),(2.5 2.8,2.5 3.2,3.5 3.2,3.5 2.8,2.5 2.8))",
+			"POLYGON((2 2.5,4 2.5,4 3.5,2 3.5,2 2.5),(2.5 2.8,2.5 3.2,3.5 3.2,3.5 2.8,2.5 2.8))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// PH2: Exterior clipped, hole entirely inside clipped region
+		{
+			"PH2",
+			"POLYGON((0 2.5,4 2.5,4 3.5,0 3.5,0 2.5),(2 2.8,2 3.2,3 3.2,3 2.8,2 2.8))",
+			"POLYGON((1 2.5,4 2.5,4 3.5,1 3.5,1 2.5),(2 2.8,2 3.2,3 3.2,3 2.8,2 2.8))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// PH3: Multiple holes, all inside R
+		{
+			"PH3",
+			"POLYGON((0 0,6 0,6 6,0 6,0 0),(2 2.5,2 3,3 3,3 2.5,2 2.5),(3.5 2.5,3.5 3,4.5 3,4.5 2.5,3.5 2.5))",
+			"POLYGON((1 2,5 2,5 4,1 4,1 2),(2 2.5,2 3,3 3,3 2.5,2 2.5),(3.5 2.5,3.5 3,4.5 3,4.5 2.5,3.5 2.5))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// PH4: Hole entirely outside R (inside exterior outside R)
+		{
+			"PH4",
+			"POLYGON((0 0,6 0,6 6,0 6,0 0),(0.2 0.2,0.2 0.8,0.8 0.8,0.8 0.2,0.2 0.2))",
+			"POLYGON((1 2,5 2,5 4,1 4,1 2))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// PH5: Hole in part of exterior that is clipped away
+		{
+			"PH5",
+			"POLYGON((2 0,4 0,4 3,2 3,2 0),(2.5 0.5,2.5 1.5,3.5 1.5,3.5 0.5,2.5 0.5))",
+			"POLYGON((2 2,4 2,4 3,2 3,2 2))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// PH6: Hole crosses one edge of R (left edge)
+		{
+			"PH6",
+			"POLYGON((-2 -2,8 -2,8 8,-2 8,-2 -2),(0 2.5,0 3.5,2 3.5,2 2.5,0 2.5))",
+			"POLYGON((1 4,1 3.5,2 3.5,2 2.5,1 2.5,1 2,5 2,5 4,1 4))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// PH7: Hole crosses two adjacent edges of R (bottom-left corner)
+		{
+			"PH7",
+			"POLYGON((-2 -2,8 -2,8 8,-2 8,-2 -2),(0 1,0 3,2 3,2 1,0 1))",
+			"POLYGON((1 3,2 3,2 2,5 2,5 4,1 4,1 3))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// PH8: Hole crosses two opposite edges of R (left and right) — splits polygon
+		{
+			"PH8",
+			"POLYGON((-2 -2,8 -2,8 8,-2 8,-2 -2),(0 2.8,0 3.2,6 3.2,6 2.8,0 2.8))",
+			"MULTIPOLYGON(((1 2,5 2,5 2.8,1 2.8,1 2)),((1 3.2,5 3.2,5 4,1 4,1 3.2)))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// PH9: Hole crosses all four edges of R, leaving top and bottom strips
+		{
+			"PH9",
+			"POLYGON((-2 -2,8 -2,8 8,-2 8,-2 -2),(0.5 2.5,0.5 3.5,5.5 3.5,5.5 2.5,0.5 2.5))",
+			"MULTIPOLYGON(((1 2,5 2,5 2.5,1 2.5,1 2)),((1 3.5,5 3.5,5 4,1 4,1 3.5)))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// PH10: Hole on R boundary, exterior extends beyond R — becomes concavity
+		{
+			"PH10",
+			"POLYGON((-2 -2,8 -2,8 8,-2 8,-2 -2),(1 2.5,1 3.5,2 3.5,2 2.5,1 2.5))",
+			"POLYGON((1 4,1 3.5,2 3.5,2 2.5,1 2.5,1 2,5 2,5 4,1 4))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// PH11: One hole inside R, another outside R
+		{
+			"PH11",
+			"POLYGON((-2 -2,8 -2,8 8,-2 8,-2 -2),(2.5 2.8,2.5 3.2,3.5 3.2,3.5 2.8,2.5 2.8),(-1.5 -1.5,-1.5 -0.5,-0.5 -0.5,-0.5 -1.5,-1.5 -1.5))",
+			"POLYGON((1 2,5 2,5 4,1 4,1 2),(2.5 2.8,2.5 3.2,3.5 3.2,3.5 2.8,2.5 2.8))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// PH12: One hole inside R, another splits polygon
+		{
+			"PH12",
+			"POLYGON((-2 -2,8 -2,8 8,-2 8,-2 -2),(2.5 2.2,2.5 2.4,3.5 2.4,3.5 2.2,2.5 2.2),(0 2.5,0 3.5,6 3.5,6 2.5,0 2.5))",
+			"MULTIPOLYGON(((1 2,5 2,5 2.5,1 2.5,1 2),(2.5 2.2,2.5 2.4,3.5 2.4,3.5 2.2,2.5 2.2)),((1 3.5,5 3.5,5 4,1 4,1 3.5)))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// PH13: Multiple holes each crossing one edge of R (left edge)
+		{
+			"PH13",
+			"POLYGON((-2 -2,8 -2,8 8,-2 8,-2 -2),(0 2.2,0 2.6,2 2.6,2 2.2,0 2.2),(0 3.4,0 3.8,2 3.8,2 3.4,0 3.4))",
+			"POLYGON((1 4,1 3.8,2 3.8,2 3.4,1 3.4,1 2.6,2 2.6,2 2.2,1 2.2,1 2,5 2,5 4,1 4))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// PH14: Two holes each crossing two opposite edges, splitting into three pieces
+		{
+			"PH14",
+			"POLYGON((-2 -2,8 -2,8 8,-2 8,-2 -2),(0 2.5,0 2.8,6 2.8,6 2.5,0 2.5),(0 3.2,0 3.5,6 3.5,6 3.2,0 3.2))",
+			"MULTIPOLYGON(((1 2,5 2,5 2.5,1 2.5,1 2)),((1 2.8,5 2.8,5 3.2,1 3.2,1 2.8)),((1 3.5,5 3.5,5 4,1 4,1 3.5)))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// PH15: Hole covers entire clipped area
+		{
+			"PH15",
+			"POLYGON((-2 -2,8 -2,8 8,-2 8,-2 -2),(0.5 1.5,0.5 4.5,5.5 4.5,5.5 1.5,0.5 1.5))",
+			"POLYGON EMPTY",
+			nil,
+		},
+
+		// PG_Z1: XYZ polygon containing R — Z is dropped, output is XY only
+		{
+			"PG_Z1",
+			"POLYGON Z((0 0 10,6 0 10,6 6 10,0 6 10,0 0 10))",
+			"POLYGON((1 2,5 2,5 4,1 4,1 2))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+
+		// MPG1: Empty MultiPolygon
+		{"MPG1", "MULTIPOLYGON EMPTY", "MULTIPOLYGON EMPTY", nil},
+		// MPG2: All component Polygons inside R
+		{
+			"MPG2",
+			"MULTIPOLYGON(((2 2.5,3 2.5,3 3,2 3,2 2.5)),((3.5 2.5,4.5 2.5,4.5 3,3.5 3,3.5 2.5)))",
+			"MULTIPOLYGON(((2 2.5,3 2.5,3 3,2 3,2 2.5)),((3.5 2.5,4.5 2.5,4.5 3,3.5 3,3.5 2.5)))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// MPG3: All component Polygons outside R
+		{
+			"MPG3",
+			"MULTIPOLYGON(((6 5,7 5,7 6,6 6,6 5)),((8 8,9 8,9 9,8 9,8 8)))",
+			"MULTIPOLYGON EMPTY",
+			nil,
+		},
+		// MPG4: Some components inside, some outside
+		{
+			"MPG4",
+			"MULTIPOLYGON(((2 2.5,3 2.5,3 3,2 3,2 2.5)),((6 5,7 5,7 6,6 6,6 5)))",
+			"MULTIPOLYGON(((2 2.5,3 2.5,3 3,2 3,2 2.5)))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// MPG5: One component partially clipped, another fully inside
+		{
+			"MPG5",
+			"MULTIPOLYGON(((0 2.5,3 2.5,3 3.5,0 3.5,0 2.5)),((3.5 2.5,4.5 2.5,4.5 3.5,3.5 3.5,3.5 2.5)))",
+			"MULTIPOLYGON(((1 2.5,3 2.5,3 3.5,1 3.5,1 2.5)),((3.5 2.5,4.5 2.5,4.5 3.5,3.5 3.5,3.5 2.5)))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// MPG6: One component becomes multiple polygons when clipped (C-shape)
+		{
+			"MPG6",
+			"MULTIPOLYGON(((0 2.5,3 2.5,3 2.8,0.5 2.8,0.5 3.2,3 3.2,3 3.5,0 3.5,0 2.5)),((4 2.5,4.5 2.5,4.5 3,4 3,4 2.5)))",
+			"MULTIPOLYGON(((1 2.5,3 2.5,3 2.8,1 2.8,1 2.5)),((1 3.2,3 3.2,3 3.5,1 3.5,1 3.2)),((4 2.5,4.5 2.5,4.5 3,4 3,4 2.5)))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// MPG7: Multiple components, each partially clipped
+		{
+			"MPG7",
+			"MULTIPOLYGON(((0 2.3,3 2.3,3 2.7,0 2.7,0 2.3)),((0 3.3,3 3.3,3 3.7,0 3.7,0 3.3)))",
+			"MULTIPOLYGON(((1 2.3,3 2.3,3 2.7,1 2.7,1 2.3)),((1 3.3,3 3.3,3 3.7,1 3.7,1 3.3)))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// MPG8: Components with holes, some holes clipped
+		{
+			"MPG8",
+			"MULTIPOLYGON(((0 2.5,3 2.5,3 3.5,0 3.5,0 2.5),(1.5 2.8,1.5 3.2,2.5 3.2,2.5 2.8,1.5 2.8)),((3.5 2.5,4.5 2.5,4.5 3.5,3.5 3.5,3.5 2.5)))",
+			"MULTIPOLYGON(((1 2.5,3 2.5,3 3.5,1 3.5,1 2.5),(1.5 2.8,1.5 3.2,2.5 3.2,2.5 2.8,1.5 2.8)),((3.5 2.5,4.5 2.5,4.5 3.5,3.5 3.5,3.5 2.5)))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// MPG9: Single component fully inside R
+		{
+			"MPG9",
+			"MULTIPOLYGON(((2 2.5,4 2.5,4 3.5,2 3.5,2 2.5)))",
+			"MULTIPOLYGON(((2 2.5,4 2.5,4 3.5,2 3.5,2 2.5)))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// MPG10: MultiPolygon containing empty Polygons
+		{
+			"MPG10",
+			"MULTIPOLYGON(((2 2.5,4 2.5,4 3.5,2 3.5,2 2.5)),EMPTY)",
+			"MULTIPOLYGON(((2 2.5,4 2.5,4 3.5,2 3.5,2 2.5)))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// MPG11: XYZ MultiPolygon, all components outside R — output is XY EMPTY
+		{
+			"MPG11",
+			"MULTIPOLYGON Z(((6 5 1,7 5 2,7 6 3,6 6 4,6 5 1)),((8 8 5,9 8 6,9 9 7,8 9 8,8 8 5)))",
+			"MULTIPOLYGON EMPTY",
+			nil,
+		},
+
+		// GC1: Empty GeometryCollection
+		{"GC1", "GEOMETRYCOLLECTION EMPTY", "GEOMETRYCOLLECTION EMPTY", nil},
+		// GC2: Contains only Points, all inside R
+		{"GC2", "GEOMETRYCOLLECTION(POINT(3 3),POINT(4 3))", "GEOMETRYCOLLECTION(POINT(3 3),POINT(4 3))", nil},
+		// GC3: Contains only Points, all outside R
+		{"GC3", "GEOMETRYCOLLECTION(POINT(0 0),POINT(6 6))", "GEOMETRYCOLLECTION EMPTY", nil},
+		// GC4: Contains mixed types, all inside R
+		{
+			"GC4",
+			"GEOMETRYCOLLECTION(POINT(3 3),LINESTRING(2 3,4 3))",
+			"GEOMETRYCOLLECTION(POINT(3 3),LINESTRING(2 3,4 3))",
+			nil,
+		},
+		// GC5: Contains mixed types, all outside R
+		{
+			"GC5",
+			"GEOMETRYCOLLECTION(POINT(0 0),LINESTRING(6 5,7 6))",
+			"GEOMETRYCOLLECTION EMPTY",
+			nil,
+		},
+		// GC6: Contains mixed types, some inside, some outside
+		{
+			"GC6",
+			"GEOMETRYCOLLECTION(POINT(3 3),POINT(0 0),LINESTRING(2 3,4 3))",
+			"GEOMETRYCOLLECTION(POINT(3 3),LINESTRING(2 3,4 3))",
+			nil,
+		},
+		// GC7: Contains Point, LineString, and Polygon (each clipped independently)
+		{
+			"GC7",
+			"GEOMETRYCOLLECTION(POINT(3 3),LINESTRING(0 3,6 3),POLYGON((2 2.5,4 2.5,4 3.5,2 3.5,2 2.5)))",
+			"GEOMETRYCOLLECTION(POINT(3 3),LINESTRING(1 3,5 3),POLYGON((2 2.5,4 2.5,4 3.5,2 3.5,2 2.5)))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// GC8: Contains nested GeometryCollection
+		{
+			"GC8",
+			"GEOMETRYCOLLECTION(POINT(3 3),GEOMETRYCOLLECTION(POINT(4 3)))",
+			"GEOMETRYCOLLECTION(POINT(3 3),GEOMETRYCOLLECTION(POINT(4 3)))",
+			nil,
+		},
+		// GC9: Contains nested GeometryCollection with mixed types
+		{
+			"GC9",
+			"GEOMETRYCOLLECTION(POINT(3 3),GEOMETRYCOLLECTION(POINT(0 0),LINESTRING(2 3,4 3)))",
+			"GEOMETRYCOLLECTION(POINT(3 3),GEOMETRYCOLLECTION(LINESTRING(2 3,4 3)))",
+			nil,
+		},
+		// GC10: All child geometries are empty
+		{
+			"GC10",
+			"GEOMETRYCOLLECTION(POINT EMPTY,LINESTRING EMPTY)",
+			"GEOMETRYCOLLECTION EMPTY",
+			nil,
+		},
+		// GC11: Contains MultiPoint, MultiLineString, MultiPolygon
+		{
+			"GC11",
+			"GEOMETRYCOLLECTION(MULTIPOINT(3 3,0 0),MULTILINESTRING((2 3,4 3)),MULTIPOLYGON(((2 2.5,4 2.5,4 3.5,2 3.5,2 2.5))))",
+			"GEOMETRYCOLLECTION(MULTIPOINT(3 3),MULTILINESTRING((2 3,4 3)),MULTIPOLYGON(((2 2.5,4 2.5,4 3.5,2 3.5,2 2.5))))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// GC12: Deeply nested GeometryCollections (3+ levels)
+		{
+			"GC12",
+			"GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 3))))",
+			"GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 3))))",
+			nil,
+		},
+		// GC13: Nested GC whose children all clip to empty
+		{
+			"GC13",
+			"GEOMETRYCOLLECTION(POINT(3 3),GEOMETRYCOLLECTION(POINT(0 0),POINT(6 6)))",
+			"GEOMETRYCOLLECTION(POINT(3 3))",
+			nil,
+		},
+		// GC14: XYZ GeometryCollection, all children outside R — output is XY EMPTY
+		{
+			"GC14",
+			"GEOMETRYCOLLECTION Z(POINT Z(0 0 1),LINESTRING Z(6 5 2,7 6 3))",
+			"GEOMETRYCOLLECTION EMPTY",
+			nil,
+		},
+
+		// NE1: Very large coordinates (outside R)
+		{"NE1", "LINESTRING(1000000000000000 3,1000000000000006 3)", "LINESTRING EMPTY", nil},
+		// NE2: Very small coordinates near float64 epsilon (inside R)
+		{"NE2", "POINT(3 3.0000000000000004)", "POINT(3 3.0000000000000004)", nil},
+		// NE3: Negative coordinates (outside R)
+		{"NE3", "LINESTRING(-3 -3,-1 -1)", "LINESTRING EMPTY", nil},
+		// NE4: Intersection parameter t very close to 0
+		{"NE4", "LINESTRING(0.999999 3,4 3)", "LINESTRING(1 3,4 3)", nil},
+		// NE5: Intersection parameter t very close to 1
+		{"NE5", "LINESTRING(2 3,5.000001 3)", "LINESTRING(2 3,5 3)", nil},
+		// NE6: Polygon vertex exactly at intersection point with R edge
+		{"NE6", "LINESTRING(1 3,5 3)", "LINESTRING(1 3,5 3)", nil},
+		// NE7: Segment nearly parallel to R edge (small angle)
+		{"NE7", "LINESTRING(0 2.0001,6 2.0001)", "LINESTRING(1 2.0001,5 2.0001)", nil},
+		// NE8: Zero-length segments in input (duplicate consecutive vertices)
+		{"NE8", "LINESTRING(2 3,2 3,4 3)", "LINESTRING(2 3,2 3,4 3)", nil},
+
+		// CD1: XY geometry clipped
+		{"CD1", "LINESTRING(0 3,6 3)", "LINESTRING(1 3,5 3)", nil},
+		// CD2: XYZ input → XY output (Z dropped)
+		{"CD2", "LINESTRING Z(0 3 0,6 3 6)", "LINESTRING(1 3,5 3)", nil},
+		// CD3: XYM input → XY output (M dropped)
+		{"CD3", "LINESTRING M(0 3 0,6 3 6)", "LINESTRING(1 3,5 3)", nil},
+		// CD4: XYZM input → XY output (Z and M dropped)
+		{"CD4", "LINESTRING ZM(0 3 0 12,6 3 6 24)", "LINESTRING(1 3,5 3)", nil},
+		// CD5: XYZ Point input → XY output
+		{"CD5", "POINT Z(3 3 7)", "POINT(3 3)", nil},
+		// CD6: XYZM Polygon input → XY output, polygon clipped
+		{
+			"CD6",
+			"POLYGON ZM((0 2.5 1 11,6 2.5 2 12,6 3.5 3 13,0 3.5 4 14,0 2.5 1 11))",
+			"POLYGON((1 2.5,5 2.5,5 3.5,1 3.5,1 2.5))",
+			[]geom.ExactEqualsOption{geom.IgnoreOrder},
+		},
+		// CD7: XYM MultiPoint input → XY output, partial survival
+		{"CD7", "MULTIPOINT M(3 3 1,0 0 2)", "MULTIPOINT(3 3)", nil},
+		// CD8: XYZ GeometryCollection input → XY output, mixed survival
+		{
+			"CD8",
+			"GEOMETRYCOLLECTION Z(POINT Z(3 3 1),LINESTRING Z(0 3 0,6 3 6))",
+			"GEOMETRYCOLLECTION(POINT(3 3),LINESTRING(1 3,5 3))",
+			nil,
+		},
+	} {
+		t.Run(tt.name, func(t *testing.T) {
+			for _, tr := range d4Transforms {
+				t.Run(tr.name, func(t *testing.T) {
+					input := test.FromWKT(t, tt.input).TransformXY(tr.fn)
+					want := test.FromWKT(t, tt.want).TransformXY(tr.fn)
+					r := rect.TransformXY(tr.fn)
+					got := geom.ClipByRect2D(input, r)
+					test.ExactEquals(t, got, want, tt.opts...)
+				})
+			}
+		})
+	}
+}
+
+func TestClipByRect2DDegenerateRect(t *testing.T) {
+	emptyRect := geom.Envelope{}
+	pointRect := geom.NewEnvelope(geom.XY{X: 3, Y: 3}, geom.XY{X: 3, Y: 3})
+	lineRect := geom.NewEnvelope(geom.XY{X: 1, Y: 3}, geom.XY{X: 5, Y: 3})
+
+	for _, tt := range []struct {
+		name  string
+		rect  geom.Envelope
+		input string
+		want  string
+	}{
+		// DR1: Empty envelope, Point
+		{"DR1", emptyRect, "POINT(3 3)", "POINT EMPTY"},
+		// DR2: Empty envelope, LineString
+		{"DR2", emptyRect, "LINESTRING(2 3,4 3)", "LINESTRING EMPTY"},
+		// DR3: Empty envelope, Polygon
+		{"DR3", emptyRect, "POLYGON((2 2,4 2,4 4,2 4,2 2))", "POLYGON EMPTY"},
+		// DR4: Point envelope, Point at same location
+		{"DR4", pointRect, "POINT(3 3)", "POINT(3 3)"},
+		// DR5: Point envelope, Point at different location
+		{"DR5", pointRect, "POINT(2 2)", "POINT EMPTY"},
+		// DR6: Point envelope, MultiPoint with some points at location
+		{"DR6", pointRect, "MULTIPOINT(3 3,2 2)", "MULTIPOINT(3 3)"},
+		// DR7: Point envelope, MultiPoint with no points at location
+		{"DR7", pointRect, "MULTIPOINT(2 2,4 4)", "MULTIPOINT EMPTY"},
+		// DR8: Point envelope, LineString through that point
+		{"DR8", pointRect, "LINESTRING(0 0,6 6)", "LINESTRING EMPTY"},
+		// DR9: Point envelope, Polygon containing that point
+		{"DR9", pointRect, "POLYGON((2 2,4 2,4 4,2 4,2 2))", "POLYGON EMPTY"},
+		// DR10: Line envelope, Point on the line
+		{"DR10", lineRect, "POINT(3 3)", "POINT(3 3)"},
+		// DR11: Line envelope, Point off the line
+		{"DR11", lineRect, "POINT(3 2)", "POINT EMPTY"},
+		// DR12: Line envelope, MultiPoint with some points on the line
+		{"DR12", lineRect, "MULTIPOINT(3 3,3 2)", "MULTIPOINT(3 3)"},
+		// DR13: Line envelope, MultiPoint with no points on the line
+		{"DR13", lineRect, "MULTIPOINT(0 0,6 6)", "MULTIPOINT EMPTY"},
+		// DR14: Line envelope, LineString crossing the line
+		{"DR14", lineRect, "LINESTRING(3 0,3 6)", "LINESTRING EMPTY"},
+		// DR15: Line envelope, LineString collinear with line
+		{"DR15", lineRect, "LINESTRING(1 3,5 3)", "LINESTRING(1 3,5 3)"},
+		// DR16: Line envelope, Polygon crossing the line
+		{"DR16", lineRect, "POLYGON((0 0,6 0,6 6,0 6,0 0))", "POLYGON EMPTY"},
+	} {
+		t.Run(tt.name, func(t *testing.T) {
+			for _, tr := range d4Transforms {
+				t.Run(tr.name, func(t *testing.T) {
+					input := test.FromWKT(t, tt.input).TransformXY(tr.fn)
+					want := test.FromWKT(t, tt.want).TransformXY(tr.fn)
+					r := tt.rect.TransformXY(tr.fn)
+					got := geom.ClipByRect2D(input, r)
+					test.ExactEquals(t, got, want)
+				})
+			}
+		})
+	}
+}
diff --git a/geom/alg_linear_interpolation.go b/geom/alg_linear_interpolation.go
index 5bd5eb46..adce8703 100644
--- a/geom/alg_linear_interpolation.go
+++ b/geom/alg_linear_interpolation.go
@@ -74,6 +74,14 @@ func lerp(a, b, t float64) float64 {
 	return math.Min(b, x)
 }
 
+// lerpXY linearly interpolates between [XY] points a and b at parameter t.
+func lerpXY(a, b XY, t float64) XY {
+	return XY{
+		X: lerp(a.X, b.X, t),
+		Y: lerp(a.Y, b.Y, t),
+	}
+}
+
 func interpolateCoords(c0, c1 Coordinates, frac float64) Coordinates {
 	return Coordinates{
 		XY: XY{
diff --git a/geom/type_sequence.go b/geom/type_sequence.go
index c2b320f8..e3521f4c 100644
--- a/geom/type_sequence.go
+++ b/geom/type_sequence.go
@@ -96,6 +96,25 @@ func (s Sequence) GetXY(i int) XY {
 	}
 }
 
+// asXYs returns the [XY] of every point location in the [Sequence], in order.
+func (s Sequence) asXYs() []XY {
+	stride := s.ctype.Dimension()
+	xys := make([]XY, 0, len(s.floats)/stride)
+	for off := 0; off < len(s.floats); off += stride {
+		xys = append(xys, XY{X: s.floats[off], Y: s.floats[off+1]})
+	}
+	return xys
+}
+
+// xysToSeq builds a [DimXY] [Sequence] from a slice of [XY].
+func xysToSeq(xys []XY) Sequence {
+	floats := make([]float64, 0, 2*len(xys))
+	for _, xy := range xys {
+		floats = append(floats, xy.X, xy.Y)
+	}
+	return NewSequence(floats, DimXY)
+}
+
 // Reverse returns a new [Sequence] containing the same point locations, but in
 // reversed order.
 func (s Sequence) Reverse() Sequence {