Add SnapToGrid methods

These methods transform the receiver geometry by snapping its control
points to a base 10 grid.

CHANGELOG.md

@@ -1,5 +1,13 @@
 # Changelog
 
+## Unreleased
+
+YYYY-MM-DD
+
+- Adds a new method `SnapToGrid` to each concrete geometry type. The method
+  returns a copy of the geometry with its control points snapped to a base 10
+  grid.
+
 ## v0.47.0
 
 2024-01-19

geom/snap_to_grid.go

@@ -0,0 +1,52 @@
+package geom
+
+import "math"
+
+func snapToGridXY(dp int) func(XY) XY {
+	return func(pt XY) XY {
+		return XY{
+			snapToGridFloat64(pt.X, dp),
+			snapToGridFloat64(pt.Y, dp),
+		}
+	}
+}
+
+// NOTE: This function would naively be implemented as:
+//
+//	return math.Round(f * math.Pow10(dp)) / math.Pow10(dp)
+//
+// However, this approach would causes two problems (which are fixed by having
+// a slightly more complex implementation):
+//
+// 1. In floating point math, numbers of the form 10^dp can be represented
+// exactly for values of dp such that 0 <= dp <= 15 (i.e. those less than
+// 2^53). Numbers of the form 10^dp can never be represented exactly for
+// negative values of dp (since the fractional part is recurring in base 2). To
+// remedy this, the function is split into "positive", "negative", and "zero"
+// dp cases.
+//
+// 2. For large values of dp, the input could overflow or underflow after being
+// multiplied by the scale factor. This causes the wrong result when the scale
+// factor is multiplied or divided out after rounding. This can be remedied by
+// detecting this case and returning the input unaltered (for an overflow) or
+// as zero (for an underflow).
+func snapToGridFloat64(f float64, dp int) float64 {
+	switch {
+	case dp > 0:
+		scale := math.Pow10(dp)
+		scaled := f * scale
+		if scaled > math.MaxFloat64 {
+			return f
+		}
+		return math.Round(scaled) / scale
+	case dp < 0:
+		scale := math.Pow10(-dp)
+		scaled := f / scale
+		if scaled == 0 {
+			return 0
+		}
+		return math.Round(scaled) * scale
+	default:
+		return math.Round(f)
+	}
+}

geom/snap_to_grid_test.go

@@ -0,0 +1,100 @@
+package geom_test
+
+import (
+	"math"
+	"strconv"
+	"testing"
+
+	"github.com/peterstace/simplefeatures/geom"
+)
+
+func TestSnapPiFloat64ToGrid(t *testing.T) {
+	for _, tc := range []struct {
+		dp   int
+		want float64
+	}{
+		{-999, 0},
+		{-99, 0},
+		{-9, 0},
+		{-2, 0},
+		{-1, 0},
+		{0, 3},
+		{1, 3.1},
+		{2, 3.14},
+		{3, 3.142},
+		{4, 3.1416},
+		{5, 3.14159},
+		{6, 3.141593},
+		{7, 3.1415927},
+		{8, 3.14159265},
+		{9, 3.141592654},
+		{10, 3.1415926536},
+		{11, 3.14159265359},
+		{12, 3.141592653590},
+		{13, 3.1415926535898},
+		{14, 3.14159265358979},
+		{15, 3.141592653589793},
+		{16, 3.141592653589793},
+		{17, 3.141592653589793},
+		{99, 3.141592653589793},
+		{999, 3.141592653589793},
+	} {
+		t.Run(strconv.Itoa(tc.dp), func(t *testing.T) {
+			pt := geom.XY{0, math.Pi}.AsPoint()
+			pt = pt.SnapToGrid(tc.dp)
+			xy, ok := pt.XY()
+			expectTrue(t, ok)
+			expectFloat64Eq(t, xy.Y, tc.want)
+		})
+	}
+}
+
+func TestSnapToGrid(t *testing.T) {
+	for i, tc := range []struct {
+		input  string
+		output string
+	}{
+		{"GEOMETRYCOLLECTION EMPTY", "GEOMETRYCOLLECTION EMPTY"},
+		{"POINT EMPTY", "POINT EMPTY"},
+		{"LINESTRING EMPTY", "LINESTRING EMPTY"},
+		{"POLYGON EMPTY", "POLYGON EMPTY"},
+		{"MULTIPOINT EMPTY", "MULTIPOINT EMPTY"},
+		{"MULTILINESTRING EMPTY", "MULTILINESTRING EMPTY"},
+		{"MULTIPOLYGON EMPTY", "MULTIPOLYGON EMPTY"},
+
+		{
+			"GEOMETRYCOLLECTION(POINT(1.11 2.22))",
+			"GEOMETRYCOLLECTION(POINT(1.1 2.2))",
+		},
+		{
+			"POINT(1.11 2.22)",
+			"POINT(1.1 2.2)",
+		},
+		{
+			"LINESTRING(1.11 2.22,3.33 4.44)",
+			"LINESTRING(1.1 2.2,3.3 4.4)",
+		},
+		{
+			"POLYGON((0.00 0.00,0.00 1.11,1.11 0.00,0.00 0.00))",
+			"POLYGON((0.0 0.0,0.0 1.1,1.1 0.0,0.0 0.0))",
+		},
+		{
+			"MULTIPOINT(1.11 2.22,3.33 4.44)",
+			"MULTIPOINT(1.1 2.2,3.3 4.4)",
+		},
+		{
+			"MULTILINESTRING((1.11 2.22,3.33 4.44),(5.55 6.66,7.77 8.88))",
+			"MULTILINESTRING((1.1 2.2,3.3 4.4),(5.6 6.7,7.8 8.9))",
+		},
+		{
+			"MULTIPOLYGON(((0.00 0.00,0.00 1.11,1.11 0.00,0.00 0.00)),((2.22 3.33,2.22 4.44,3.33 3.33,2.22 3.33)))",
+			"MULTIPOLYGON(((0.0 0.0,0.0 1.1,1.1 0.0,0.0 0.0)),((2.2 3.3,2.2 4.4,3.3 3.3,2.2 3.3)))",
+		},
+	} {
+		t.Run(strconv.Itoa(i), func(t *testing.T) {
+			in := geomFromWKT(t, tc.input)
+			got := in.SnapToGrid(1)
+			expectGeomEqWKT(t, got, tc.output)
+		})
+	}
+}

geom/type_geometry.go

@@ -973,3 +973,35 @@ func (g Geometry) Validate() error {
 		panic("unknown type: " + g.Type().String())
 	}
 }
+
+// SnapToGrid returns a copy of the geometry with all coordinates snapped to a
+// base 10 grid.
+//
+// The grid spacing is specified by the number of decimal places to round to
+// (with negative decimal places being allowed). E.g., a decimalPlaces value of
+// 2 would cause all coordinates to be rounded to the nearest 0.01, and a
+// decimalPlaces of -1 would cause all coordinates to be rounded to the nearest
+// 10.
+//
+// Returned geometries may be invalid due to snapping, even if the input
+// geometry was valid.
+func (g Geometry) SnapToGrid(decimalPlaces int) Geometry {
+	switch g.gtype {
+	case TypeGeometryCollection:
+		return g.MustAsGeometryCollection().SnapToGrid(decimalPlaces).AsGeometry()
+	case TypePoint:
+		return g.MustAsPoint().SnapToGrid(decimalPlaces).AsGeometry()
+	case TypeLineString:
+		return g.MustAsLineString().SnapToGrid(decimalPlaces).AsGeometry()
+	case TypePolygon:
+		return g.MustAsPolygon().SnapToGrid(decimalPlaces).AsGeometry()
+	case TypeMultiPoint:
+		return g.MustAsMultiPoint().SnapToGrid(decimalPlaces).AsGeometry()
+	case TypeMultiLineString:
+		return g.MustAsMultiLineString().SnapToGrid(decimalPlaces).AsGeometry()
+	case TypeMultiPolygon:
+		return g.MustAsMultiPolygon().SnapToGrid(decimalPlaces).AsGeometry()
+	default:
+		panic("unknown type: " + g.Type().String())
+	}
+}

geom/type_geometry_collection.go

@@ -552,3 +552,18 @@ func (c GeometryCollection) Simplify(threshold float64, nv ...NoValidate) (Geome
 	}
 	return NewGeometryCollection(geoms), nil
 }
+
+// SnapToGrid returns a copy of the GeometryCollection with all coordinates
+// snapped to a base 10 grid.
+//
+// The grid spacing is specified by the number of decimal places to round to
+// (with negative decimal places being allowed). E.g., a decimalPlaces value of
+// 2 would cause all coordinates to be rounded to the nearest 0.01, and a
+// decimalPlaces of -1 would cause all coordinates to be rounded to the nearest
+// 10.
+//
+// Returned GeometryCollections may be invalid due to snapping, even if the
+// input geometry was valid.
+func (c GeometryCollection) SnapToGrid(decimalPlaces int) GeometryCollection {
+	return c.TransformXY(snapToGridXY(decimalPlaces))
+}

geom/type_line_string.go

@@ -473,3 +473,18 @@ func (s LineString) InterpolateEvenlySpacedPoints(n int) MultiPoint {
 	}
 	return NewMultiPoint(pts)
 }
+
+// SnapToGrid returns a copy of the LineString with all coordinates snapped to
+// a base 10 grid.
+//
+// The grid spacing is specified by the number of decimal places to round to
+// (with negative decimal places being allowed). E.g., a decimalPlaces value of
+// 2 would cause all coordinates to be rounded to the nearest 0.01, and a
+// decimalPlaces of -1 would cause all coordinates to be rounded to the nearest
+// 10.
+//
+// Returned LineStrings may be invalid due to snapping, even if the input
+// geometry was valid.
+func (s LineString) SnapToGrid(decimalPlaces int) LineString {
+	return s.TransformXY(snapToGridXY(decimalPlaces))
+}

geom/type_multi_line_string.go

@@ -510,3 +510,18 @@ func (m MultiLineString) Simplify(threshold float64) MultiLineString {
 	}
 	return NewMultiLineString(lss)
 }
+
+// SnapToGrid returns a copy of the MultiLineString with all coordinates
+// snapped to a base 10 grid.
+//
+// The grid spacing is specified by the number of decimal places to round to
+// (with negative decimal places being allowed). E.g., a decimalPlaces value of
+// 2 would cause all coordinates to be rounded to the nearest 0.01, and a
+// decimalPlaces of -1 would cause all coordinates to be rounded to the nearest
+// 10.
+//
+// Returned MultiLineStrings may be invalid due to snapping, even if the input
+// geometry was valid.
+func (m MultiLineString) SnapToGrid(decimalPlaces int) MultiLineString {
+	return m.TransformXY(snapToGridXY(decimalPlaces))
+}

geom/type_multi_point.go

@@ -338,3 +338,15 @@ func (m MultiPoint) Summary() string {
 func (m MultiPoint) String() string {
 	return m.Summary()
 }
+
+// SnapToGrid returns a copy of the MultiPoint with all coordinates snapped to
+// a base 10 grid.
+//
+// The grid spacing is specified by the number of decimal places to round to
+// (with negative decimal places being allowed). E.g., a decimalPlaces value of
+// 2 would cause all coordinates to be rounded to the nearest 0.01, and a
+// decimalPlaces of -1 would cause all coordinates to be rounded to the nearest
+// 10.
+func (m MultiPoint) SnapToGrid(decimalPlaces int) MultiPoint {
+	return m.TransformXY(snapToGridXY(decimalPlaces))
+}

geom/type_multi_polygon.go

@@ -569,3 +569,18 @@ func (m MultiPolygon) Simplify(threshold float64, nv ...NoValidate) (MultiPolygo
 	}
 	return simpl, nil
 }
+
+// SnapToGrid returns a MultiPolygon of the geometry with all coordinates
+// snapped to a base 10 grid.
+//
+// The grid spacing is specified by the number of decimal places to round to
+// (with negative decimal places being allowed). E.g., a decimalPlaces value of
+// 2 would cause all coordinates to be rounded to the nearest 0.01, and a
+// decimalPlaces of -1 would cause all coordinates to be rounded to the nearest
+// 10.
+//
+// Returned MultiPolygons may be invalid due to snapping, even if the input
+// geometry was valid.
+func (m MultiPolygon) SnapToGrid(decimalPlaces int) MultiPolygon {
+	return m.TransformXY(snapToGridXY(decimalPlaces))
+}

geom/type_point.go

@@ -268,3 +268,15 @@ func (p Point) Summary() string {
 func (p Point) String() string {
 	return p.Summary()
 }
+
+// SnapToGrid returns a copy of the Point with all coordinates snapped to a
+// base 10 grid.
+//
+// The grid spacing is specified by the number of decimal places to round to
+// (with negative decimal places being allowed). E.g., a decimalPlaces value of
+// 2 would cause all coordinates to be rounded to the nearest 0.01, and a
+// decimalPlaces of -1 would cause all coordinates to be rounded to the nearest
+// 10.
+func (p Point) SnapToGrid(decimalPlaces int) Point {
+	return p.TransformXY(snapToGridXY(decimalPlaces))
+}

geom/type_polygon.go

@@ -681,3 +681,18 @@ func (p Polygon) Simplify(threshold float64, nv ...NoValidate) (Polygon, error)
 	}
 	return simpl, nil
 }
+
+// SnapToGrid returns a copy of the Polygon with all coordinates snapped to a
+// base 10 grid.
+//
+// The grid spacing is specified by the number of decimal places to round to
+// (with negative decimal places being allowed). E.g., a decimalPlaces value of
+// 2 would cause all coordinates to be rounded to the nearest 0.01, and a
+// decimalPlaces of -1 would cause all coordinates to be rounded to the nearest
+// 10.
+//
+// Returned Polygons may be invalid due to snapping, even if the input geometry
+// was valid.
+func (p Polygon) SnapToGrid(decimalPlaces int) Polygon {
+	return p.TransformXY(snapToGridXY(decimalPlaces))
+}