geo/geogfn: implement ST_Azimuth({geometry,geometry})

Fixes cockroachdb#48887

Release note (sql change): This PR implement adds the ST_Azimuth({geometry,geometry})

docs/generated/sql/functions.md

@@ -776,6 +776,9 @@ has no relationship with the commit order of concurrent transactions.</p>
 </span></td></tr>
 <tr><td><a name="st_astext"></a><code>st_astext(geometry: geometry, maximum_decimal_digits: <a href="int.html">int</a>) &rarr; <a href="string.html">string</a></code></td><td><span class="funcdesc"><p>Returns the WKT representation of a given Geometry. The maximum_decimal_digits parameter controls the maximum decimal digits to print after the <code>.</code>. Use -1 to print as many digits as possible.</p>
 </span></td></tr>
+<tr><td><a name="st_azimuth"></a><code>st_azimuth(geometry_a: geometry, geometry_b: geometry) &rarr; <a href="float.html">float</a></code></td><td><span class="funcdesc"><p>Returns the azimuth in radians of the segment defined by the given point geometries, or NULL if the two points are coincident.</p>
+<p>The azimuth is angle is referenced from north, and is positive clockwise: North = 0; East = π/2; South = π; West = 3π/2.</p>
+</span></td></tr>
 <tr><td><a name="st_buffer"></a><code>st_buffer(geometry: geometry, distance: <a href="float.html">float</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns a Geometry that represents all points whose distance is less than or equal to the given distance
 from the given Geometry.</p>
 <p>This function utilizes the GEOS module.</p>

pkg/geo/geomfn/binary_predicates.go

@@ -11,10 +11,25 @@
 package geomfn
 
 import (
+	"fmt"
+	"math"
+
 	"github.com/cockroachdb/cockroach/pkg/geo"
 	"github.com/cockroachdb/cockroach/pkg/geo/geos"
+	"github.com/twpayne/go-geom"
 )
 
+// Azimuth returns the azimuth in radians of the segment defined by the given point geometries.
+// The azimuth is angle is referenced from north, and is positive clockwise.
+// North = 0; East = π/2; South = π; West = 3π/2.
+func Azimuth(a *geom.Point, b *geom.Point) (float64, error) {
+	if a.X() == b.X() && a.Y() == b.Y() {
+		return 0, fmt.Errorf("points are the same")
+	}
+
+	return math.Mod(2*math.Pi+math.Pi/2-math.Atan2(b.Y()-a.Y(), b.X()-a.X()), 2*math.Pi), nil
+}
+
 // Covers returns whether geometry A covers geometry B.
 func Covers(a *geo.Geometry, b *geo.Geometry) (bool, error) {
 	if a.SRID() != b.SRID() {

pkg/geo/geomfn/binary_predicates_test.go

@@ -16,6 +16,7 @@ import (
 
 	"github.com/cockroachdb/cockroach/pkg/geo"
 	"github.com/stretchr/testify/require"
+	"github.com/twpayne/go-geom"
 )
 
 var (
@@ -29,6 +30,54 @@ var (
 	middleLine           = geo.MustParseGeometry("LINESTRING(-0.5 0.5, 0.5 0.5)")
 )
 
+func TestAzimuth(t *testing.T) {
+	testCases := []struct {
+		a           *geom.Point
+		b           *geom.Point
+		expected    float64
+		expectedErr error
+	}{
+		{
+			geom.NewPointFlat(geom.XY, []float64{0, 0}),
+			geom.NewPointFlat(geom.XY, []float64{0, 0}),
+			0,
+			fmt.Errorf("points are the same"),
+		},
+		{
+			geom.NewPointFlat(geom.XY, []float64{0, 0}),
+			geom.NewPointFlat(geom.XY, []float64{1, 1}),
+			0.7853981633974483, // math.Pi * 1 / 4
+			nil,
+		},
+		{
+			geom.NewPointFlat(geom.XY, []float64{0, 0}),
+			geom.NewPointFlat(geom.XY, []float64{1, 0}),
+			1.5707963267948966, // math.Pi * 2 / 4
+			nil,
+		},
+		{
+			geom.NewPointFlat(geom.XY, []float64{0, 0}),
+			geom.NewPointFlat(geom.XY, []float64{1, -1}),
+			2.356194490192344, // almost math.Pi * 3 / 4
+			nil,
+		},
+		{
+			geom.NewPointFlat(geom.XY, []float64{0, 0}),
+			geom.NewPointFlat(geom.XY, []float64{0, 1}),
+			0,
+			nil,
+		},
+	}
+
+	for i, tc := range testCases {
+		t.Run(fmt.Sprintf("tc:%d", i), func(t *testing.T) {
+			r, err := Azimuth(tc.a, tc.b)
+			require.Equal(t, tc.expected, r)
+			require.Equal(t, tc.expectedErr, err)
+		})
+	}
+}
+
 func TestCovers(t *testing.T) {
 	testCases := []struct {
 		a        *geo.Geometry

pkg/sql/logictest/testdata/logic_test/geospatial

@@ -112,6 +112,16 @@ SELECT ST_AsText(ST_Project('POINT(0 0)'::geography, 100000, radians(45.0)))
 ----
 POINT (0.635231029125537 0.639472334729198)
 
+statement error Argument must be POINT geometries
+SELECT ST_Azimuth('POLYGON((0 0, 0 0, 0 0, 0 0))'::geometry, 'POLYGON((0 0, 0 0, 0 0, 0 0))'::geometry)
+
+query RR
+SELECT
+	degrees(ST_Azimuth(ST_Point(25, 45), ST_Point(75, 100))) AS degA_B,
+	degrees(ST_Azimuth(ST_Point(75, 100), ST_Point(25, 45))) AS degB_A
+----
+42.2736890060937  222.273689006094
+
 subtest cast_test
 
 query T

pkg/sql/sem/builtins/geo_builtins.go

@@ -1769,6 +1769,46 @@ Note If the result has zero or one points, it will be returned as a POINT. If it
 	//
 	// Binary functions
 	//
+	"st_azimuth": makeBuiltin(
+		defProps(),
+		geometryOverload2(
+			func(ctx *tree.EvalContext, a, b *tree.DGeometry) (tree.Datum, error) {
+				aGeomT, err := a.AsGeomT()
+				if err != nil {
+					return nil, err
+				}
+
+				aPoint, ok := aGeomT.(*geom.Point)
+				if !ok {
+					return nil, errors.Newf("Argument must be POINT geometries")
+				}
+
+				bGeomT, err := b.AsGeomT()
+				if err != nil {
+					return nil, err
+				}
+
+				bPoint, ok := bGeomT.(*geom.Point)
+				if !ok {
+					return nil, errors.Newf("Argument must be POINT geometries")
+				}
+
+				azimuth, err := geomfn.Azimuth(aPoint, bPoint)
+				if err != nil {
+					return nil, err
+				}
+
+				return tree.NewDFloat(tree.DFloat(azimuth)), nil
+			},
+			types.Float,
+			infoBuilder{
+				info: `Returns the azimuth in radians of the segment defined by the given point geometries, or NULL if the two points are coincident.
+
+The azimuth is angle is referenced from north, and is positive clockwise: North = 0; East = π/2; South = π; West = 3π/2.`,
+			},
+			tree.VolatilityImmutable,
+		),
+	),
 	"st_distance": makeBuiltin(
 		defProps(),
 		geometryOverload2(