diff --git a/docs/generated/sql/functions.md b/docs/generated/sql/functions.md index 0a1750091957..d5b660b8febd 100644 --- a/docs/generated/sql/functions.md +++ b/docs/generated/sql/functions.md @@ -976,6 +976,11 @@ has no relationship with the commit order of concurrent transactions.

st_relate(geometry_a: geometry, geometry_b: geometry, pattern: string) → bool

Returns whether the DE-9IM spatial relation between geometry_a and geometry_b matches the DE-9IM pattern.

This function utilizes the GEOS module.

 +st_segmentize(geography: geography, max_segment_length: float) → geography

turns a modified Geography having no segment longer than the given max_segment_length meters. +The calculations are done on a sphere.

+

This function utilizes the S2 library for spherical calculations.

+

This function will automatically use any available index.

+ st_srid(geography: geography) → int

Returns the Spatial Reference Identifier (SRID) for the ST_Geography as defined in spatial_ref_sys table.

 st_srid(geometry: geometry) → int

Returns the Spatial Reference Identifier (SRID) for the ST_Geometry as defined in spatial_ref_sys table.

diff --git a/pkg/geo/geo.go b/pkg/geo/geo.go index a11f005ee075..041fac966dbb 100644 --- a/pkg/geo/geo.go +++ b/pkg/geo/geo.go @@ -188,6 +188,15 @@ func NewGeography(spatialObject geopb.SpatialObject) *Geography { return &Geography{SpatialObject: spatialObject} } +// NewGeographyFromGeom Geography Object from geom.T. +func NewGeographyFromGeom(g geom.T) (*Geography, error) { + spatialObject, err := spatialObjectFromGeom(g) + if err != nil { + return nil, err + } + return NewGeography(spatialObject), nil +} + // ParseGeography parses a Geography from a given text. func ParseGeography(str string) (*Geography, error) { spatialObject, err := parseAmbiguousText(str, geopb.DefaultGeographySRID) diff --git a/pkg/geo/geogfn/segmentize.go b/pkg/geo/geogfn/segmentize.go new file mode 100644 index 000000000000..372925d47c35 --- /dev/null +++ b/pkg/geo/geogfn/segmentize.go @@ -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 +} diff --git a/pkg/geo/geogfn/segmentize_test.go b/pkg/geo/geogfn/segmentize_test.go new file mode 100644 index 000000000000..b57ce5014eea --- /dev/null +++ b/pkg/geo/geogfn/segmentize_test.go @@ -0,0 +1,143 @@ +// 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" +) + +var 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)))", + }, +} + +func TestSegmentize(t *testing.T) { + 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.resultantCoordinate, convertedPoints) + }) + } +} diff --git a/pkg/sql/logictest/testdata/logic_test/geospatial b/pkg/sql/logictest/testdata/logic_test/geospatial index ffe7af75b22b..56258d855a1f 100644 --- a/pkg/sql/logictest/testdata/logic_test/geospatial +++ b/pkg/sql/logictest/testdata/logic_test/geospatial @@ -805,6 +805,26 @@ Square overlapping left and right square Square (left) Square overlapping left and right square Square (right) true true true Square overlapping left and right square Square overlapping left and right square true true true +# ST_Segmentize +query TTT +SELECT + dsc, + ST_AsText(ST_Segmentize(geog, 100000)), + ST_AsText(ST_Segmentize(geog, 50000)) +FROM geog_operators_test +ORDER BY dsc +---- +Faraway point POINT (5 5) POINT (5 5) +Line going through left and right square LINESTRING (-0.5 0.5, -0.00000000000000009939611878359099 0.5000190382262164, 0.5 0.5) LINESTRING (-0.5 0.5, -0.25000000036247944 0.500014278647005, -0.00000000000000009939611878359099 0.5000190382262164, 0.2500000003624792 0.5000142786470051, 0.5 0.5) +NULL NULL NULL +Point middle of Left Square POINT (-0.5 0.5) POINT (-0.5 0.5) +Point middle of Right Square POINT (0.5 0.5) POINT (0.5 0.5) +Square (left) POLYGON ((-1 0, -0.5000000000000001 0, 0 0, 0 0.5, 0 1, -0.4999999999999998 1.0000380706528733, -1 1, -0.9999999999999998 0.5000000000000001, -1 0)) POLYGON ((-1 0, -0.7499999999999998 0, -0.5000000000000001 0, -0.2499999999999997 0, 0 0, 0 0.25, 0 0.5, 0 0.75, 0 1, -0.2499999985501929 1.0000285529443267, -0.4999999999999998 1.0000380706528733, -0.7500000014498067 1.0000285529443265, -1 1, -1 0.7499999999999998, -0.9999999999999998 0.5000000000000001, -0.9999999999999998 0.25, -1 0)) +Square (right) POLYGON ((0 0, 0.5 0, 1 0, 1 0.4999999999999999, 1 1, 0.5 1.0000380706528733, 0 1, 0 0.5000000000000001, 0 0)) 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.750000001449807 1.0000285529443267, 0.5 1.0000380706528733, 0.2499999985501931 1.0000285529443267, 0 1, 0 0.7499999999999998, 0 0.5000000000000001, 0 0.2499999999999997, 0 0)) +Square overlapping left and right square POLYGON ((-0.1 0, 0.44999999999999996 0, 1 0, 1 0.4999999999999999, 1 1, 0.44999999999999996 1.0000460657968335, -0.1 1, -0.1 0.5000000000000001, -0.1 0)) POLYGON ((-0.1 0, 0.17500000000000007 0, 0.44999999999999996 0, 0.7249999999999999 0, 1 0, 0.9999999999999998 0.25, 1 0.4999999999999999, 0.9999999999999998 0.7499999999999999, 1 1, 0.7250000019297163 1.0000345492812595, 0.44999999999999996 1.0000460657968335, 0.17499999807028374 1.0000345492812592, -0.1 1, -0.1 0.75, -0.1 0.5000000000000001, -0.10000000000000002 0.2499999999999997, -0.1 0)) + +statement error st_segmentize\(\): maximum segment length must be positive +SELECT ST_Segmentize('point(0 0)'::geography, 0); subtest pg_extension diff --git a/pkg/sql/sem/builtins/geo_builtins.go b/pkg/sql/sem/builtins/geo_builtins.go index 4775cf4e2cf1..5a774da33018 100644 --- a/pkg/sql/sem/builtins/geo_builtins.go +++ b/pkg/sql/sem/builtins/geo_builtins.go @@ -1129,6 +1129,36 @@ Note ST_Perimeter is only valid for Polygon - use ST_Length for LineString.`, Volatility: tree.VolatilityImmutable, }, ), + + // + // Geography modification + // + + "st_segmentize": makeBuiltin( + defProps(), + tree.Overload{ + Types: tree.ArgTypes{ + {"geography", types.Geography}, + {"max_segment_length", types.Float}, + }, + ReturnType: tree.FixedReturnType(types.Geography), + Fn: func(_ *tree.EvalContext, args tree.Datums) (tree.Datum, error) { + g := args[0].(*tree.DGeography) + segmentMaxLength := float64(*args[1].(*tree.DFloat)) + segGeometry, err := geogfn.Segmentize(g.Geography, segmentMaxLength) + if err != nil { + return nil, err + } + return tree.NewDGeography(segGeometry), nil + }, + Info: infoBuilder{ + info: `turns a modified Geography having no segment longer than the given max_segment_length meters. + The calculations are done on a sphere.`, + libraryUsage: usesS2, + canUseIndex: true, + }.String(), + Volatility: tree.VolatilityImmutable, + }), } // geometryOverload1 hides the boilerplate for builtins operating on one geometry.