From e4e7b0b0fd251c1c61436650924c4005fdea81aa Mon Sep 17 00:00:00 2001 From: Peter Stace Date: Sun, 28 Jan 2024 07:28:29 +1100 Subject: [PATCH] Add `SnapToGrid` methods These methods transform the receiver geometry by snapping its control points to a base 10 grid. --- CHANGELOG.md | 8 +++ geom/snap_to_grid.go | 52 ++++++++++++++++ geom/snap_to_grid_test.go | 100 +++++++++++++++++++++++++++++++ geom/type_geometry.go | 32 ++++++++++ geom/type_geometry_collection.go | 15 +++++ geom/type_line_string.go | 15 +++++ geom/type_multi_line_string.go | 15 +++++ geom/type_multi_point.go | 12 ++++ geom/type_multi_polygon.go | 15 +++++ geom/type_point.go | 12 ++++ geom/type_polygon.go | 15 +++++ 11 files changed, 291 insertions(+) create mode 100644 geom/snap_to_grid.go create mode 100644 geom/snap_to_grid_test.go diff --git a/CHANGELOG.md b/CHANGELOG.md index 0f9c47e5..4ef2d422 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/geom/snap_to_grid.go b/geom/snap_to_grid.go new file mode 100644 index 00000000..c2534039 --- /dev/null +++ b/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) + } +} diff --git a/geom/snap_to_grid_test.go b/geom/snap_to_grid_test.go new file mode 100644 index 00000000..4374db8a --- /dev/null +++ b/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) + }) + } +} diff --git a/geom/type_geometry.go b/geom/type_geometry.go index 5d2e0a48..f414d50c 100644 --- a/geom/type_geometry.go +++ b/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()) + } +} diff --git a/geom/type_geometry_collection.go b/geom/type_geometry_collection.go index 56d9fefc..6690cd58 100644 --- a/geom/type_geometry_collection.go +++ b/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)) +} diff --git a/geom/type_line_string.go b/geom/type_line_string.go index c206bc88..88805905 100644 --- a/geom/type_line_string.go +++ b/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)) +} diff --git a/geom/type_multi_line_string.go b/geom/type_multi_line_string.go index 6097cc1b..8318e441 100644 --- a/geom/type_multi_line_string.go +++ b/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)) +} diff --git a/geom/type_multi_point.go b/geom/type_multi_point.go index d3a8d1f9..5383a789 100644 --- a/geom/type_multi_point.go +++ b/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)) +} diff --git a/geom/type_multi_polygon.go b/geom/type_multi_polygon.go index 93fe1e7c..d66899f3 100644 --- a/geom/type_multi_polygon.go +++ b/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)) +} diff --git a/geom/type_point.go b/geom/type_point.go index ceb3c890..2208a8ef 100644 --- a/geom/type_point.go +++ b/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)) +} diff --git a/geom/type_polygon.go b/geom/type_polygon.go index 99178833..da46b2bd 100644 --- a/geom/type_polygon.go +++ b/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)) +}