Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

geo/geogfn: implement ST_Project #49949

Merged
merged 1 commit into from
Jun 8, 2020
Merged

geo/geogfn: implement ST_Project #49949

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions docs/generated/sql/functions.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1110,6 +1110,13 @@ given Geometry.</p>
</span></td></tr>
<tr><td><a name="st_polygonfromwkb"></a><code>st_polygonfromwkb(wkb: <a href="bytes.html">bytes</a>, srid: <a href="int.html">int</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns the Geometry from a WKB representation with an SRID. If the shape underneath is not Polygon, NULL is returned.</p>
</span></td></tr>
<tr><td><a name="st_project"></a><code>st_project(geography: geography, distance: <a href="float.html">float</a>, azimuth: <a href="float.html">float</a>) &rarr; geography</code></td><td><span class="funcdesc"><p>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.</p>
<p>The distance is given in meters. Negative values are supported.</p>
<p>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.</p>
</span></td></tr>
<tr><td><a name="st_relate"></a><code>st_relate(geometry_a: geometry, geometry_b: geometry) &rarr; <a href="string.html">string</a></code></td><td><span class="funcdesc"><p>Returns the DE-9IM spatial relation between geometry_a and geometry_b.</p>
<p>This function utilizes the GEOS module.</p>
</span></td></tr>
Expand Down
69 changes: 69 additions & 0 deletions pkg/geo/geogfn/unary_operators.go
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
)
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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
}
34 changes: 34 additions & 0 deletions pkg/geo/geogfn/unary_operators_test.go
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand Down Expand Up @@ -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)
})
}
}
27 changes: 26 additions & 1 deletion pkg/geo/geographiclib/geographiclib.go
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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))
}
29 changes: 29 additions & 0 deletions pkg/geo/geographiclib/geographiclib_test.go
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,10 @@
package geographiclib

import (
"math"
"testing"

"github.com/golang/geo/s1"
"github.com/golang/geo/s2"
"github.com/stretchr/testify/require"
)
Expand Down Expand Up @@ -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)
})
}
}
5 changes: 5 additions & 0 deletions pkg/sql/logictest/testdata/logic_test/geospatial
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
50 changes: 50 additions & 0 deletions pkg/sql/sem/builtins/geo_builtins.go
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
)
Expand Down Expand Up @@ -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).
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if you wanted extra new lines here to display in HTML format, you'll need an extra new line here.
This new line is otherwise displayed as spaces (see autogenerated code and how the

tags are missing.)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've seen functions.md, but I can't find anything wrong about the format. Can you explain it in more detail?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more. 


<p>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.</p>

if you wanted new line between them, you'd want


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

then you need a new line between each here. if you're ok with just space difference, then this is ok. i'm up for either.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for explain. I want to keep it as it is.

Negative azimuth values and values greater than 2π (360 degrees) are supported.`,
}.String(),
Volatility: tree.VolatilityImmutable,
},
),

//
// Unary functions.
Expand Down