Merge #49949

49949: geo/geogfn: implement ST_Project r=otan a=hueypark

Fixes #48402

Release note (sql change): This PR implement adds the ST_Project({geography,float8,float8})

Co-authored-by: Jaewan Park <jaewan.huey.park@gmail.com>

docs/generated/sql/functions.md

@@ -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>

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
+}

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)
+		})
+	}
+}

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))
+}

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)
+		})
+	}
+}

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

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.