forked from cockroachdb/cockroach
-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
builtin: Implements ST_Segmentize builtin function
Fixes: cockroachdb#48368 This PR adds ST_Segmentize({geography,float8}) builtin function which allows modify given geography such that no segment longer than the given max_segment_length. Release note (sql change): This PR implements ST_Segmentize builtin function.
- Loading branch information
1 parent
a40921b
commit a36b0f9
Showing
6 changed files
with
367 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,161 @@ | ||
| // Copyright 2020 The Cockroach Authors. | ||
| // | ||
| // Use of this software is governed by the Business Source License | ||
| // included in the file licenses/BSL.txt. | ||
| // | ||
| // As of the Change Date specified in that file, in accordance with | ||
| // the Business Source License, use of this software will be governed | ||
| // by the Apache License, Version 2.0, included in the file | ||
| // licenses/APL.txt. | ||
|
|
||
| package geogfn | ||
|
|
||
| import ( | ||
| "math" | ||
|
|
||
| "github.com/cockroachdb/cockroach/pkg/geo" | ||
| "github.com/cockroachdb/cockroach/pkg/geo/geographiclib" | ||
| "github.com/cockroachdb/errors" | ||
| "github.com/golang/geo/s2" | ||
| "github.com/twpayne/go-geom" | ||
| ) | ||
|
|
||
| // Segmentize return modified Geography having no segment longer | ||
| // that given maximum segment length. | ||
| func Segmentize(geography *geo.Geography, segmentMaxLength float64) (*geo.Geography, error) { | ||
| spheroid := geographiclib.WGS84Spheroid | ||
| // Convert segmentMaxLength to s1.ChordAngle as further calculation | ||
| // is done considering segmentMaxLength as s1.ChordAngle. | ||
| segmentMaxLength = segmentMaxLength / spheroid.SphereRadius | ||
| geometry, err := geography.AsGeomT() | ||
| if err != nil { | ||
| return nil, err | ||
| } | ||
| segGeometry, err := segmentizeGeom(geometry, segmentMaxLength) | ||
| if err != nil { | ||
| return nil, err | ||
| } | ||
| return geo.NewGeographyFromGeom(segGeometry) | ||
| } | ||
|
|
||
| // segmentizeGeom returns a modified geom.T having no segment longer than | ||
| // the given maximum segment length. | ||
| func segmentizeGeom(geometry geom.T, segmentMaxLength float64) (geom.T, error) { | ||
| if segmentMaxLength <= 0 { | ||
| return nil, errors.Newf("maximum segment length must be positive") | ||
| } | ||
| switch geometry := geometry.(type) { | ||
| case *geom.Point, *geom.MultiPoint: | ||
| return geometry, nil | ||
| case *geom.LineString: | ||
| var allFlatCoordinates []float64 | ||
| for pointIdx := 1; pointIdx < geometry.NumCoords(); pointIdx++ { | ||
| allFlatCoordinates = append( | ||
| allFlatCoordinates, | ||
| segmentizeCoords(geometry.Coord(pointIdx-1), geometry.Coord(pointIdx), segmentMaxLength)..., | ||
| ) | ||
| } | ||
| // Appending end point as it wasn't included in the iteration of coordinates. | ||
| allFlatCoordinates = append(allFlatCoordinates, geometry.Coord(geometry.NumCoords()-1)...) | ||
| return geom.NewLineStringFlat(geom.XY, allFlatCoordinates).SetSRID(geometry.SRID()), nil | ||
| case *geom.MultiLineString: | ||
| segMultiLine := geom.NewMultiLineString(geom.XY).SetSRID(geometry.SRID()) | ||
| for lineIdx := 0; lineIdx < geometry.NumLineStrings(); lineIdx++ { | ||
| l, err := segmentizeGeom(geometry.LineString(lineIdx), segmentMaxLength) | ||
| if err != nil { | ||
| return nil, err | ||
| } | ||
| err = segMultiLine.Push(l.(*geom.LineString)) | ||
| if err != nil { | ||
| return nil, err | ||
| } | ||
| } | ||
| return segMultiLine, nil | ||
| case *geom.LinearRing: | ||
| var allFlatCoordinates []float64 | ||
| for pointIdx := 1; pointIdx < geometry.NumCoords(); pointIdx++ { | ||
| allFlatCoordinates = append( | ||
| allFlatCoordinates, | ||
| segmentizeCoords(geometry.Coord(pointIdx-1), geometry.Coord(pointIdx), segmentMaxLength)..., | ||
| ) | ||
| } | ||
| // Appending end point as it wasn't included in the iteration of coordinates. | ||
| allFlatCoordinates = append(allFlatCoordinates, geometry.Coord(geometry.NumCoords()-1)...) | ||
| return geom.NewLinearRingFlat(geom.XY, allFlatCoordinates).SetSRID(geometry.SRID()), nil | ||
| case *geom.Polygon: | ||
| segPolygon := geom.NewPolygon(geom.XY).SetSRID(geometry.SRID()) | ||
| for loopIdx := 0; loopIdx < geometry.NumLinearRings(); loopIdx++ { | ||
| l, err := segmentizeGeom(geometry.LinearRing(loopIdx), segmentMaxLength) | ||
| if err != nil { | ||
| return nil, err | ||
| } | ||
| err = segPolygon.Push(l.(*geom.LinearRing)) | ||
| if err != nil { | ||
| return nil, err | ||
| } | ||
| } | ||
| return segPolygon, nil | ||
| case *geom.MultiPolygon: | ||
| segMultiPolygon := geom.NewMultiPolygon(geom.XY).SetSRID(geometry.SRID()) | ||
| for polygonIdx := 0; polygonIdx < geometry.NumPolygons(); polygonIdx++ { | ||
| p, err := segmentizeGeom(geometry.Polygon(polygonIdx), segmentMaxLength) | ||
| if err != nil { | ||
| return nil, err | ||
| } | ||
| err = segMultiPolygon.Push(p.(*geom.Polygon)) | ||
| if err != nil { | ||
| return nil, err | ||
| } | ||
| } | ||
| return segMultiPolygon, nil | ||
| case *geom.GeometryCollection: | ||
| segGeomCollection := geom.NewGeometryCollection().SetSRID(geometry.SRID()) | ||
| for geoIdx := 0; geoIdx < geometry.NumGeoms(); geoIdx++ { | ||
| g, err := segmentizeGeom(geometry.Geom(geoIdx), segmentMaxLength) | ||
| if err != nil { | ||
| return nil, err | ||
| } | ||
| err = segGeomCollection.Push(g) | ||
| if err != nil { | ||
| return nil, err | ||
| } | ||
| } | ||
| return segGeomCollection, nil | ||
| } | ||
| return nil, errors.Newf("unknown type: %T", geometry) | ||
| } | ||
|
|
||
| // segmentizeCoords inserts multiple points between given two-coordinate and | ||
| // return resultant points as []float64. Such that distance between any two | ||
| // points is less than given maximum segment's length, the total number of | ||
| // segments is the power of 2, and all the segments are of the same length. | ||
| // NOTE: List of points does not consist of end point. | ||
| func segmentizeCoords(a geom.Coord, b geom.Coord, maxSegmentLength float64) []float64 { | ||
| // Converted geom.Coord into s2.Point so we can segmentize the coordinates. | ||
| point1 := s2.PointFromLatLng(s2.LatLngFromDegrees(a.Y(), a.X())) | ||
| point2 := s2.PointFromLatLng(s2.LatLngFromDegrees(b.Y(), b.X())) | ||
|
|
||
| var segmentizedPoints []s2.Point | ||
| distanceBetweenPoints := s2.ChordAngleBetweenPoints(point1, point2).Angle().Radians() | ||
| if maxSegmentLength <= distanceBetweenPoints { | ||
| // This calculation is to determine the total number of segment between given | ||
| // 2 coordinates. For that fraction by segment must be less than or equal to | ||
| // the fraction of max segment length to distance between point, since the | ||
| // total number of segment must be power of 2 therefore we can write as | ||
| // 1 / (2^n)[numberOfSegmentToCreate] <= segmentMaxLength / distanceBetweenPoints < 1 / (2^(n-1)) | ||
| // (2^n)[numberOfSegmentToCreate] >= distanceBetweenPoints / segmentMaxLength > 2^(n-1) | ||
| // therefore n = ceil(log2(segmentMaxLength/distanceBetweenPoints)). Hence | ||
| // numberOfSegmentToCreate = 2^(ceil(log2(segmentMaxLength/distanceBetweenPoints))). | ||
| numberOfSegmentToCreate := int(math.Pow(2, math.Ceil(math.Log2(distanceBetweenPoints/maxSegmentLength)))) | ||
| for pointInserted := 1; pointInserted < numberOfSegmentToCreate; pointInserted++ { | ||
| newPoint := s2.Interpolate(float64(pointInserted)/float64(numberOfSegmentToCreate), point1, point2) | ||
| segmentizedPoints = append(segmentizedPoints, newPoint) | ||
| } | ||
| } | ||
| allSegmentizedCoordinates := a.Clone() | ||
| for _, point := range segmentizedPoints { | ||
| latLng := s2.LatLngFromPoint(point) | ||
| allSegmentizedCoordinates = append(allSegmentizedCoordinates, latLng.Lng.Degrees(), latLng.Lat.Degrees()) | ||
| } | ||
| return allSegmentizedCoordinates | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,142 @@ | ||
| // Copyright 2020 The Cockroach Authors. | ||
| // | ||
| // Use of this software is governed by the Business Source License | ||
| // included in the file licenses/BSL.txt. | ||
| // | ||
| // As of the Change Date specified in that file, in accordance with | ||
| // the Business Source License, use of this software will be governed | ||
| // by the Apache License, Version 2.0, included in the file | ||
| // licenses/APL.txt. | ||
|
|
||
| package geogfn | ||
|
|
||
| import ( | ||
| "fmt" | ||
| "testing" | ||
|
|
||
| "github.com/cockroachdb/cockroach/pkg/geo" | ||
| "github.com/stretchr/testify/require" | ||
| "github.com/twpayne/go-geom" | ||
| ) | ||
|
|
||
| func TestSegmentize(t *testing.T) { | ||
| segmentizeTestCases := []struct { | ||
| wkt string | ||
| maxSegmentLength float64 | ||
| segmentizedGeogWKT string | ||
| }{ | ||
| { | ||
| wkt: "POINT(1.0 1.0)", | ||
| maxSegmentLength: 10000.0, | ||
| segmentizedGeogWKT: "POINT(1.0 1.0)", | ||
| }, | ||
| { | ||
| wkt: "LINESTRING(1.0 1.0, 2.0 2.0, 3.0 3.0)", | ||
| maxSegmentLength: 100000.0, | ||
| segmentizedGeogWKT: "LINESTRING (1 1, 1.4998857365616758 1.5000570914791973, 2 2, 2.4998094835255658 2.500095075163195, 3 3)", | ||
| }, | ||
| { | ||
| wkt: "LINESTRING(1.0 1.0, 2.0 2.0, 3.0 3.0)", | ||
| maxSegmentLength: 500000.0, | ||
| segmentizedGeogWKT: "LINESTRING(1.0 1.0, 2.0 2.0, 3.0 3.0)", | ||
| }, | ||
| { | ||
| wkt: "POLYGON((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 0.0))", | ||
| maxSegmentLength: 100000.0, | ||
| segmentizedGeogWKT: "POLYGON ((0 0, 0.5 0, 1 0, 1 0.4999999999999999, 1 1, 0.4999619199226218 0.5000190382261059, 0 0))", | ||
| }, | ||
| { | ||
| wkt: "POLYGON((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 0.0), (0.1 0.1, 0.2 0.1, 0.2 0.2, 0.1 0.1))", | ||
| maxSegmentLength: 50000.0, | ||
| segmentizedGeogWKT: "POLYGON ((0 0, 0.25 0, 0.5 0, 0.75 0, 1 0, 0.9999999999999998 0.25, 1 0.4999999999999999, 0.9999999999999998 0.7499999999999999, 1 1, 0.7499666792978344 0.7500166590029188, 0.4999619199226218 0.5000190382261059, 0.24997620022356357 0.2500118986534327, 0 0), (0.1 0.1, 0.2 0.1, 0.2 0.2, 0.1 0.1))", | ||
| }, | ||
| { | ||
| wkt: "POLYGON((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 0.0), (0.1 0.1, 0.2 0.1, 0.2 0.2, 0.1 0.1))", | ||
| maxSegmentLength: 1000000.0, | ||
| segmentizedGeogWKT: "POLYGON((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 0.0), (0.1 0.1, 0.2 0.1, 0.2 0.2, 0.1 0.1))", | ||
| }, | ||
| { | ||
| wkt: "MULTIPOINT((1.0 1.0), (2.0 2.0))", | ||
| maxSegmentLength: 100000.0, | ||
| segmentizedGeogWKT: "MULTIPOINT((1.0 1.0), (2.0 2.0))", | ||
| }, | ||
| { | ||
| wkt: "MULTILINESTRING((1.0 1.0, 2.0 2.0, 3.0 3.0), (6.0 6.0, 7.0 6.0))", | ||
| maxSegmentLength: 100000.0, | ||
| segmentizedGeogWKT: "MULTILINESTRING ((1 1, 1.4998857365616758 1.5000570914791973, 2 2, 2.4998094835255658 2.500095075163195, 3 3), (6 6, 6.5 6.000226803574719, 7 6))", | ||
| }, | ||
| { | ||
| wkt: "MULTIPOLYGON(((3.0 3.0, 4.0 3.0, 4.0 4.0, 3.0 3.0)), ((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 0.0), (0.1 0.1, 0.2 0.1, 0.2 0.2, 0.1 0.1)))", | ||
| maxSegmentLength: 100000.0, | ||
| segmentizedGeogWKT: "MULTIPOLYGON (((3 3, 3.500000000000001 3.0001140264716564, 4 3, 3.999999999999999 3.5, 4 4, 3.4997331141752333 3.5001329429928956, 3 3)), ((0 0, 0.5 0, 1 0, 1 0.4999999999999999, 1 1, 0.4999619199226218 0.5000190382261059, 0 0), (0.1 0.1, 0.2 0.1, 0.2 0.2, 0.1 0.1)))", | ||
| }, | ||
| { | ||
| wkt: "GEOMETRYCOLLECTION (POINT (40 10),LINESTRING (10 10, 20 20, 10 40),POLYGON ((40 40, 20 45, 45 30, 40 40)))", | ||
| maxSegmentLength: 1000000.0, | ||
| segmentizedGeogWKT: "GEOMETRYCOLLECTION (POINT (40 10), LINESTRING (10 10, 14.88248913028353 15.054670903122675, 20 20, 17.84783705611024 25.063693381403255, 15.510294720403053 30.09369304154877, 12.922758839587525 35.07789165618019, 10 40), POLYGON ((40 40, 30.404184608679003 42.93646050008845, 20 45, 27.258698904670663 41.785178536245354, 33.78296512799812 38.158856049677816, 39.6629587979805 34.20684491196235, 45 30, 42.6532474552537 35.02554212454631, 40 40)))", | ||
| }, | ||
| } | ||
| for _, test := range segmentizeTestCases { | ||
| t.Run(fmt.Sprintf("%s, maximum segment length: %v", test.wkt, test.maxSegmentLength), func(t *testing.T) { | ||
| geog, err := geo.ParseGeography(test.wkt) | ||
| require.NoError(t, err) | ||
| modifiedGeog, err := Segmentize(geog, test.maxSegmentLength) | ||
| require.NoError(t, err) | ||
| expectedGeog, err := geo.ParseGeography(test.segmentizedGeogWKT) | ||
| require.NoError(t, err) | ||
| require.Equal(t, expectedGeog, modifiedGeog) | ||
| }) | ||
| } | ||
| // Test for segment maximum length as negative. | ||
| t.Run(fmt.Sprintf("%s, maximum segment length: %v", segmentizeTestCases[0].wkt, -1.0), func(t *testing.T) { | ||
| geog, err := geo.ParseGeography(segmentizeTestCases[0].wkt) | ||
| require.NoError(t, err) | ||
| _, err = Segmentize(geog, -1) | ||
| require.EqualError(t, err, "maximum segment length must be positive") | ||
| }) | ||
| } | ||
|
|
||
| func TestSegmentizeCoords(t *testing.T) { | ||
| testCases := []struct { | ||
| desc string | ||
| a geom.Coord | ||
| b geom.Coord | ||
| segmentMaxLength float64 | ||
| resultantCoordinates []float64 | ||
| }{ | ||
| { | ||
| desc: `Coordinate(0, 0) to Coordinate(85, 85), 0.781600222084459`, | ||
| a: geom.Coord{0, 0}, | ||
| b: geom.Coord{85, 85}, | ||
| segmentMaxLength: 0.781600222084459, | ||
| resultantCoordinates: []float64{0, 0, 4.924985039227315, 44.568038920632105}, | ||
| }, | ||
| { | ||
| desc: `Coordinate(0, 0) to Coordinate(85, 85), 0.7816000651234446`, | ||
| a: geom.Coord{0, 0}, | ||
| b: geom.Coord{85, 85}, | ||
| segmentMaxLength: 0.7816000651234446, | ||
| resultantCoordinates: []float64{0, 0, 2.0486953806277866, 22.302074138733936, 4.924985039227315, 44.568038920632105, 11.655816669136822, 66.66485041602017}, | ||
| }, | ||
| { | ||
| desc: `Coordinate(85, 85) to Coordinate(0, 0), 0.29502092024628396`, | ||
| a: geom.Coord{85, 85}, | ||
| b: geom.Coord{0, 0}, | ||
| segmentMaxLength: 0.29502092024628396, | ||
| resultantCoordinates: []float64{85, 85, 22.871662720021178, 77.3609628116894, 11.655816669136824, 66.66485041602017, 7.329091976767396, 55.658764687902504, 4.924985039227315, 44.56803892063211, 3.299940624127866, 33.443216802941045, 2.048695380627787, 22.302074138733943, 0.9845421446758968, 11.1527721155093}, | ||
| }, | ||
| { | ||
| desc: `Coordinate(85, 85) to Coordinate(0, 0), -1`, | ||
| a: geom.Coord{85, 85}, | ||
| b: geom.Coord{0, 0}, | ||
| segmentMaxLength: -1, | ||
| resultantCoordinates: []float64{85, 85}, | ||
| }, | ||
| } | ||
| for _, test := range testCases { | ||
| t.Run(test.desc, func(t *testing.T) { | ||
| convertedPoints := segmentizeCoords(test.a, test.b, test.segmentMaxLength) | ||
| require.Equal(t, test.resultantCoordinates, convertedPoints) | ||
| }) | ||
| } | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters