From 5cbdb0dcc42320f8a1e0b0183a33ceea0b00eb05 Mon Sep 17 00:00:00 2001 From: Jaewan Park Date: Tue, 9 Jun 2020 06:54:19 +0900 Subject: [PATCH] geo/geogfn: implement ST_Project Fixes #48402 Release note (sql change): This PR implement adds the ST_Project({geography,float8,float8}) --- docs/generated/sql/functions.md | 7 ++ pkg/geo/geogfn/unary_operators.go | 69 +++++++++++++++++++ pkg/geo/geogfn/unary_operators_test.go | 34 +++++++++ pkg/geo/geographiclib/geographiclib.go | 27 +++++++- pkg/geo/geographiclib/geographiclib_test.go | 29 ++++++++ .../logictest/testdata/logic_test/geospatial | 5 ++ pkg/sql/sem/builtins/geo_builtins.go | 50 ++++++++++++++ 7 files changed, 220 insertions(+), 1 deletion(-) diff --git a/docs/generated/sql/functions.md b/docs/generated/sql/functions.md index 14153d12e99d..4f4e046e7cc8 100644 --- a/docs/generated/sql/functions.md +++ b/docs/generated/sql/functions.md @@ -1110,6 +1110,13 @@ given Geometry.

st_polygonfromwkb(wkb: bytes, srid: int) → geometry

Returns the Geometry from a WKB representation with an SRID. If the shape underneath is not Polygon, NULL is returned.

 +st_project(geography: geography, distance: float, azimuth: float) → geography

Returns a point projected from a start point along a geodesic using a given distance and azimuth (bearing). +This is known as the direct geodesic problem.

+

The distance is given in meters. Negative values are supported.

+

The azimuth (also known as heading or bearing) is given in radians. It is measured clockwise from true north (azimuth zero). +East is azimuth π/2 (90 degrees); south is azimuth π (180 degrees); west is azimuth 3π/2 (270 degrees). +Negative azimuth values and values greater than 2π (360 degrees) are supported.

+ st_relate(geometry_a: geometry, geometry_b: geometry) → string

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

This function utilizes the GEOS module.

 diff --git a/pkg/geo/geogfn/unary_operators.go b/pkg/geo/geogfn/unary_operators.go index e5ce3ccf7c05..128c7f76a175 100644 --- a/pkg/geo/geogfn/unary_operators.go +++ b/pkg/geo/geogfn/unary_operators.go @@ -11,9 +11,12 @@ 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/s1" "github.com/golang/geo/s2" "github.com/twpayne/go-geom" ) @@ -90,6 +93,43 @@ func Length(g *geo.Geography, useSphereOrSpheroid UseSphereOrSpheroid) (float64, return length(regions, useSphereOrSpheroid) } +// Project returns calculate a projected point given a source point, a distance and a azimuth. +func Project(point *geom.Point, distance float64, azimuth s1.Angle) (*geom.Point, error) { + spheroid := geographiclib.WGS84Spheroid + + // Normalize distance to be positive. + if distance < 0.0 { + distance = -distance + azimuth += math.Pi + } + + // Normalize azimuth + azimuth = azimuth.Normalized() + + // Check the distance validity. + if distance > (math.Pi * spheroid.Radius) { + return nil, errors.Newf("distance must not be greater than %f", math.Pi*spheroid.Radius) + } + + // Convert to ta geodetic point. + x := point.X() + y := point.Y() + + projected := spheroid.Project( + s2.LatLngFromDegrees(x, y), + distance, + azimuth, + ) + + return geom.NewPointFlat( + geom.XY, + []float64{ + float64(projected.Lng.Normalized()) * 180.0 / math.Pi, + normalizeLatitude(float64(projected.Lat)) * 180.0 / math.Pi, + }, + ), nil +} + // length returns the sum of the lengtsh and perimeters in the shapes of the Geography. // In OGC parlance, length returns both LineString lengths _and_ Polygon perimeters. func length(regions []s2.Region, useSphereOrSpheroid UseSphereOrSpheroid) (float64, error) { @@ -128,3 +168,32 @@ func length(regions []s2.Region, useSphereOrSpheroid UseSphereOrSpheroid) (float } return totalLength, nil } + +// normalizeLatitude convert a latitude to the range of -Pi/2, Pi/2. +func normalizeLatitude(lat float64) float64 { + if lat > 2.0*math.Pi { + lat = math.Remainder(lat, 2.0*math.Pi) + } + + if lat < -2.0*math.Pi { + lat = math.Remainder(lat, -2.0*math.Pi) + } + + if lat > math.Pi { + lat = math.Pi - lat + } + + if lat < -1.0*math.Pi { + lat = -1.0*math.Pi - lat + } + + if lat > math.Pi*2 { + lat = math.Pi - lat + } + + if lat < -1.0*math.Pi*2 { + lat = -1.0*math.Pi - lat + } + + return lat +} diff --git a/pkg/geo/geogfn/unary_operators_test.go b/pkg/geo/geogfn/unary_operators_test.go index 6149ad6762cb..a6918407ce32 100644 --- a/pkg/geo/geogfn/unary_operators_test.go +++ b/pkg/geo/geogfn/unary_operators_test.go @@ -15,7 +15,9 @@ import ( "testing" "github.com/cockroachdb/cockroach/pkg/geo" + "github.com/golang/geo/s1" "github.com/stretchr/testify/require" + "github.com/twpayne/go-geom" ) type unaryOperatorExpectedResult struct { @@ -219,3 +221,35 @@ func TestLength(t *testing.T) { }) } } + +func TestProject(t *testing.T) { + var testCases = []struct { + desc string + point *geom.Point + distance float64 + azimuth float64 + projected *geom.Point + }{ + { + "POINT(0 0), 100000, radians(45)", + geom.NewPointFlat(geom.XY, []float64{0, 0}), + 100000, + 45 * math.Pi / 180.0, + geom.NewPointFlat(geom.XY, []float64{0.6352310291255374, 0.6394723347291977}), + }, + } + + for _, tc := range testCases { + t.Run(tc.desc, func(t *testing.T) { + projected, err := Project(tc.point, tc.distance, s1.Angle(tc.azimuth)) + require.NoError(t, err) + require.Equalf( + t, + tc.projected, + projected, + "expected %f, found %f", + &tc.projected, + projected) + }) + } +} diff --git a/pkg/geo/geographiclib/geographiclib.go b/pkg/geo/geographiclib/geographiclib.go index e1fc0a67eb34..c63ba65747c8 100644 --- a/pkg/geo/geographiclib/geographiclib.go +++ b/pkg/geo/geographiclib/geographiclib.go @@ -18,7 +18,12 @@ package geographiclib // #include "geographiclib.h" import "C" -import "github.com/golang/geo/s2" +import ( + "math" + + "github.com/golang/geo/s1" + "github.com/golang/geo/s2" +) var ( // WGS84Spheroid represents the default WGS84 ellipsoid. @@ -109,3 +114,23 @@ func (s *Spheroid) AreaAndPerimeter(points []s2.Point) (area float64, perimeter ) return float64(areaDouble), float64(perimeterDouble) } + +// Project returns computes the location of the projected point. +// +// Using the direct geodesic problem from GeographicLib (Karney 2013). +func (s *Spheroid) Project(point s2.LatLng, distance float64, azimuth s1.Angle) s2.LatLng { + var lat, lng C.double + + C.geod_direct( + &s.cRepr, + C.double(point.Lat.Degrees()), + C.double(point.Lng.Degrees()), + C.double(azimuth*180.0/math.Pi), + C.double(distance), + &lat, + &lng, + nil, + ) + + return s2.LatLngFromDegrees(float64(lat), float64(lng)) +} diff --git a/pkg/geo/geographiclib/geographiclib_test.go b/pkg/geo/geographiclib/geographiclib_test.go index 0d051b86365f..b383db8921e7 100644 --- a/pkg/geo/geographiclib/geographiclib_test.go +++ b/pkg/geo/geographiclib/geographiclib_test.go @@ -11,8 +11,10 @@ package geographiclib import ( + "math" "testing" + "github.com/golang/geo/s1" "github.com/golang/geo/s2" "github.com/stretchr/testify/require" ) @@ -101,3 +103,30 @@ func TestAreaAndPerimeter(t *testing.T) { }) } } + +func TestProject(t *testing.T) { + testCases := []struct { + desc string + spheroid Spheroid + point s2.LatLng + distance float64 + azimuth float64 + project s2.LatLng + }{ + { + desc: "{0,0} project to 100000, radians(45.0) on WGS84Spheroid", + spheroid: WGS84Spheroid, + point: s2.LatLng{Lat: 0, Lng: 0}, + distance: 100000, + azimuth: 45 * math.Pi / 180.0, + project: s2.LatLng{Lat: 0.011160897716439782, Lng: 0.011086872969072624}, + }, + } + + for _, tc := range testCases { + t.Run(tc.desc, func(t *testing.T) { + project := tc.spheroid.Project(tc.point, tc.distance, s1.Angle(tc.azimuth)) + require.Equal(t, tc.project, project) + }) + } +} diff --git a/pkg/sql/logictest/testdata/logic_test/geospatial b/pkg/sql/logictest/testdata/logic_test/geospatial index 75b829a5cb78..84923a4b83df 100644 --- a/pkg/sql/logictest/testdata/logic_test/geospatial +++ b/pkg/sql/logictest/testdata/logic_test/geospatial @@ -98,6 +98,11 @@ SELECT ST_AsText(p) FROM (VALUES POINT (1 2) POINT (3 4) +query T +SELECT ST_AsText(ST_Project('POINT(0 0)'::geography, 100000, radians(45.0))) +---- +POINT (0.6352310291255374 0.6394723347291977) + subtest cast_test query T diff --git a/pkg/sql/sem/builtins/geo_builtins.go b/pkg/sql/sem/builtins/geo_builtins.go index 6ed97ff883f1..f3e552952975 100644 --- a/pkg/sql/sem/builtins/geo_builtins.go +++ b/pkg/sql/sem/builtins/geo_builtins.go @@ -27,6 +27,7 @@ import ( "github.com/cockroachdb/cockroach/pkg/util/errorutil/unimplemented" "github.com/cockroachdb/cockroach/pkg/util/json" "github.com/cockroachdb/errors" + "github.com/golang/geo/s1" "github.com/twpayne/go-geom" "github.com/twpayne/go-geom/encoding/ewkb" ) @@ -854,6 +855,55 @@ var geoBuiltins = map[string]builtinDefinition{ tree.VolatilityImmutable, ), ), + "st_project": makeBuiltin( + defProps(), + tree.Overload{ + Types: tree.ArgTypes{ + {"geography", types.Geography}, + {"distance", types.Float}, + {"azimuth", types.Float}, + }, + ReturnType: tree.FixedReturnType(types.Geography), + Fn: func(ctx *tree.EvalContext, args tree.Datums) (tree.Datum, error) { + g := args[0].(*tree.DGeography) + distance := float64(*args[1].(*tree.DFloat)) + azimuth := float64(*args[2].(*tree.DFloat)) + + geomT, err := g.AsGeomT() + if err != nil { + return nil, err + } + + point, ok := geomT.(*geom.Point) + if !ok { + return nil, errors.Newf("ST_Project(geography) is only valid for point inputs") + } + + projected, err := geogfn.Project(point, distance, s1.Angle(azimuth)) + if err != nil { + return nil, err + } + + geog, err := geo.NewGeographyFromGeom(projected) + if err != nil { + return nil, err + } + + return &tree.DGeography{Geography: geog}, nil + }, + Info: infoBuilder{ + info: `Returns a point projected from a start point along a geodesic using a given distance and azimuth (bearing). +This is known as the direct geodesic problem. + +The distance is given in meters. Negative values are supported. + +The azimuth (also known as heading or bearing) is given in radians. It is measured clockwise from true north (azimuth zero). +East is azimuth π/2 (90 degrees); south is azimuth π (180 degrees); west is azimuth 3π/2 (270 degrees). +Negative azimuth values and values greater than 2π (360 degrees) are supported.`, + }.String(), + Volatility: tree.VolatilityImmutable, + }, + ), // // Unary functions.