[SEDONA-656] Add ST_Project (#1594)

* feat: add ST_Project

* fix: snowflake tests

* fix: snowflake tests 2

* address comments

* chore: format python code

common/src/main/java/org/apache/sedona/common/Functions.java

@@ -32,6 +32,7 @@
 import org.apache.sedona.common.sphere.Spheroid;
 import org.apache.sedona.common.subDivide.GeometrySubDivider;
 import org.apache.sedona.common.utils.*;
+import org.locationtech.jts.algorithm.Angle;
 import org.locationtech.jts.algorithm.MinimumBoundingCircle;
 import org.locationtech.jts.algorithm.Orientation;
 import org.locationtech.jts.algorithm.construct.LargestEmptyCircle;
@@ -1353,6 +1354,45 @@ public static Integer dimension(Geometry geometry) {
     return dimension;
   }
 
+  public static Geometry project(Geometry point, double distance, double azimuth, boolean lenient) {
+    if (!point.getClass().getSimpleName().equals("Point")) {
+      if (lenient) {
+        return point.getFactory().createPoint();
+      } else {
+        throw new IllegalArgumentException(
+            String.format(
+                "Input geometry is %s. It should be a Point type geometry",
+                point.getClass().getSimpleName()));
+      }
+    }
+
+    // Normalize azimuth if it is out of (-360, 360) range
+    // by calculating the number of orbits and subtracting it
+    int orbit = (int) Math.floor(azimuth / Angle.PI_TIMES_2);
+    azimuth -= Angle.PI_TIMES_2 * orbit;
+    // Convert azimuth to conventional slope
+    double slope = Angle.PI_TIMES_2 - azimuth + Angle.PI_OVER_2;
+    if (slope > Angle.PI_TIMES_2) slope -= Angle.PI_TIMES_2;
+    if (slope < -Angle.PI_TIMES_2) slope += Angle.PI_TIMES_2;
+
+    Coordinate projectedCoordinate = Angle.project(point.getCoordinate(), slope, distance);
+
+    if (Functions.hasZ(point)) {
+      projectedCoordinate.setZ(point.getCoordinate().getZ());
+    }
+
+    if (Functions.hasM(point)) {
+      CoordinateXYZM projectedCoordinateM = new CoordinateXYZM(projectedCoordinate);
+      projectedCoordinateM.setM(point.getCoordinate().getM());
+      return point.getFactory().createPoint(projectedCoordinateM);
+    }
+    return point.getFactory().createPoint(projectedCoordinate);
+  }
+
+  public static Geometry project(Geometry point, double distance, double azimuth) {
+    return project(point, distance, azimuth, false);
+  }
+
   /**
    * get the coordinates of a geometry and transform to Google s2 cell id
    *

common/src/test/java/org/apache/sedona/common/FunctionsTest.java

@@ -658,6 +658,41 @@ public void dimensionGeometryEmpty() {
     assertEquals(actualResult, expectedResult);
   }
 
+  @Test
+  public void project() throws ParseException {
+    Geometry point = GEOMETRY_FACTORY.createPoint(new Coordinate(0, 0));
+    String actual = Functions.asWKT(Functions.project(point, 100000, Math.toRadians(45.0)));
+    String expected = "POINT (70710.67811865476 70710.67811865475)";
+    assertEquals(expected, actual);
+
+    actual =
+        Functions.asWKT(Functions.project(Constructors.makeEnvelope(0, 1, 0, 1), 10, 10, true));
+    expected = "POINT EMPTY";
+    assertEquals(expected, actual);
+
+    point = Constructors.geomFromWKT("POINT Z(10 15 12)", 1111);
+    Geometry actualPoint = Functions.project(point, 1000, Math.toRadians(300.0));
+    actual = Functions.asWKT(actualPoint);
+    expected = "POINT Z(-856.0254037844385 515.0000000000003 12)";
+    assertEquals(expected, actual);
+    assertEquals(1111, actualPoint.getSRID());
+
+    point = Constructors.geomFromWKT("POINT M(10 15 12)", 1111);
+    actual = Functions.asWKT(Functions.project(point, 1000, Math.toRadians(300.0)));
+    expected = "POINT M(-856.0254037844385 515.0000000000003 12)";
+    assertEquals(expected, actual);
+
+    point = Constructors.geomFromWKT("POINT ZM(10 15 12 2)", 1111);
+    actual = Functions.asWKT(Functions.project(point, 1000, Math.toRadians(300.0)));
+    expected = "POINT ZM(-856.0254037844385 515.0000000000003 12 2)";
+    assertEquals(expected, actual);
+
+    point = Constructors.geomFromWKT("POINT(2 -1)", 0);
+    actual = Functions.asWKT(Functions.project(point, 100, Math.toRadians(470)));
+    expected = Functions.asWKT(Functions.project(point, 100, Math.toRadians(110)));
+    assertEquals(expected, actual);
+  }
+
   private static boolean intersects(Set<?> s1, Set<?> s2) {
     Set<?> copy = new HashSet<>(s1);
     copy.retainAll(s2);

docs/api/flink/Function.md

@@ -3139,6 +3139,48 @@ Output:
 GEOMETRYCOLLECTION (POLYGON ((0 2, 1 3, 2 4, 2 3, 2 2, 1 2, 0 2)), POLYGON ((2 2, 2 3, 2 4, 3 3, 4 2, 3 2, 2 2)))
 ```
 
+## ST_Project
+
+Introduction: Calculates a new point location given a starting point, distance, and azimuth. The azimuth indicates the direction, expressed in radians, and is measured in a clockwise manner starting from true north. The system can handle azimuth values that are negative or exceed 2π (360 degrees). The optional `lenient` parameter prevents an error if the input geometry is not a Point. Its default value is `false`.
+
+Format:
+
+```
+ST_Project(point: Geometry, distance: Double, azimuth: Double, lenient: Boolean = False)
+```
+
+```
+ST_Project(point: Geometry, distance: Double, Azimuth: Double)
+```
+
+Since: `v1.7.0`
+
+SQL Example:
+
+```sql
+SELECT ST_Project(ST_GeomFromText('POINT (10 15)'), 100, radians(90))
+```
+
+Output:
+
+```
+POINT (110 14.999999999999975)
+```
+
+SQL Example:
+
+```sql
+SELECT ST_Project(
+        ST_GeomFromText('POLYGON ((1 5, 1 1, 3 3, 5 3, 1 5))'),
+        25, radians(270), true)
+```
+
+Output:
+
+```
+POINT EMPTY
+```
+
 ## ST_ReducePrecision
 
 Introduction: Reduce the decimals places in the coordinates of the geometry to the given number of decimal places. The last decimal place will be rounded.

docs/api/snowflake/vector-data/Function.md

@@ -2388,6 +2388,46 @@ Output:
 GEOMETRYCOLLECTION (POLYGON ((0 2, 1 3, 2 4, 2 3, 2 2, 1 2, 0 2)), POLYGON ((2 2, 2 3, 2 4, 3 3, 4 2, 3 2, 2 2)))
 ```
 
+## ST_Project
+
+Introduction: Calculates a new point location given a starting point, distance, and azimuth. The azimuth indicates the direction, expressed in radians, and is measured in a clockwise manner starting from true north. The system can handle azimuth values that are negative or exceed 2π (360 degrees). The optional `lenient` parameter prevents an error if the input geometry is not a Point. Its default value is `false`.
+
+Format:
+
+```
+ST_Project(point: Geometry, distance: Double, azimuth: Double, lenient: Boolean = False)
+```
+
+```
+ST_Project(point: Geometry, distance: Double, Azimuth: Double)
+```
+
+SQL Example:
+
+```sql
+SELECT ST_Project(ST_GeomFromText('POINT (10 15)'), 100, radians(90))
+```
+
+Output:
+
+```
+POINT (110 14.999999999999975)
+```
+
+SQL Example:
+
+```sql
+SELECT ST_Project(
+        ST_GeomFromText('POLYGON ((1 5, 1 1, 3 3, 5 3, 1 5))'),
+        25, radians(270), true)
+```
+
+Output:
+
+```
+POINT EMPTY
+```
+
 ## ST_ReducePrecision
 
 Introduction: Reduce the decimals places in the coordinates of the geometry to the given number of decimal places. The last decimal place will be rounded. This function was called ST_PrecisionReduce in versions prior to v1.5.0.

docs/api/sql/Function.md

@@ -3219,6 +3219,48 @@ Output:
 GEOMETRYCOLLECTION (POLYGON ((0 2, 1 3, 2 4, 2 3, 2 2, 1 2, 0 2)), POLYGON ((2 2, 2 3, 2 4, 3 3, 4 2, 3 2, 2 2)))
 ```
 
+## ST_Project
+
+Introduction: Calculates a new point location given a starting point, distance, and azimuth. The azimuth indicates the direction, expressed in radians, and is measured in a clockwise manner starting from true north. The system can handle azimuth values that are negative or exceed 2π (360 degrees). The optional `lenient` parameter prevents an error if the input geometry is not a Point. Its default value is `false`.
+
+Format:
+
+```
+ST_Project(point: Geometry, distance: Double, azimuth: Double, lenient: Boolean = False)
+```
+
+```
+ST_Project(point: Geometry, distance: Double, Azimuth: Double)
+```
+
+Since: `v1.7.0`
+
+SQL Example:
+
+```sql
+SELECT ST_Project(ST_GeomFromText('POINT (10 15)'), 100, radians(90))
+```
+
+Output:
+
+```
+POINT (110 14.999999999999975)
+```
+
+SQL Example:
+
+```sql
+SELECT ST_Project(
+        ST_GeomFromText('POLYGON ((1 5, 1 1, 3 3, 5 3, 1 5))'),
+        25, radians(270), true)
+```
+
+Output:
+
+```
+POINT EMPTY
+```
+
 ## ST_ReducePrecision
 
 Introduction: Reduce the decimals places in the coordinates of the geometry to the given number of decimal places. The last decimal place will be rounded. This function was called ST_PrecisionReduce in versions prior to v1.5.0.

flink/src/main/java/org/apache/sedona/flink/Catalog.java

@@ -161,6 +161,7 @@ public static UserDefinedFunction[] getFuncs() {
       new Functions.ST_Points(),
       new Functions.ST_Polygon(),
       new Functions.ST_Polygonize(),
+      new Functions.ST_Project(),
       new Functions.ST_MakePolygon(),
       new Functions.ST_MakeValid(),
       new Functions.ST_MaxDistance(),

flink/src/main/java/org/apache/sedona/flink/expressions/Functions.java

@@ -1190,6 +1190,29 @@ public Geometry eval(
     }
   }
 
+  public static class ST_Project extends ScalarFunction {
+    @DataTypeHint(value = "RAW", bridgedTo = org.locationtech.jts.geom.Geometry.class)
+    public Geometry eval(
+        @DataTypeHint(value = "RAW", bridgedTo = org.locationtech.jts.geom.Geometry.class)
+            Object o1,
+        @DataTypeHint(value = "Double") Double distance,
+        @DataTypeHint(value = "Double") Double azimuth,
+        @DataTypeHint("Boolean") Boolean lenient) {
+      Geometry point = (Geometry) o1;
+      return org.apache.sedona.common.Functions.project(point, distance, azimuth, lenient);
+    }
+
+    @DataTypeHint(value = "RAW", bridgedTo = org.locationtech.jts.geom.Geometry.class)
+    public Geometry eval(
+        @DataTypeHint(value = "RAW", bridgedTo = org.locationtech.jts.geom.Geometry.class)
+            Object o1,
+        @DataTypeHint(value = "Double") Double distance,
+        @DataTypeHint(value = "Double") Double azimuth) {
+      Geometry point = (Geometry) o1;
+      return org.apache.sedona.common.Functions.project(point, distance, azimuth);
+    }
+  }
+
   public static class ST_MakeValid extends ScalarFunction {
     @DataTypeHint(value = "RAW", bridgedTo = org.locationtech.jts.geom.Geometry.class)
     public Geometry eval(

flink/src/test/java/org/apache/sedona/flink/FunctionTest.java

@@ -774,6 +774,21 @@ public void testPointOnSurface() {
     assertEquals("POINT (-117.99 32.01)", result.toString());
   }
 
+  @Test
+  public void testProject() {
+    Table pointTable = createPointTable(testDataSize);
+    Table surfaceTable =
+        pointTable.select(
+            call(
+                Functions.ST_Project.class.getSimpleName(),
+                $(pointColNames[0]),
+                100,
+                Math.toRadians(45)));
+    Geometry result = (Geometry) first(surfaceTable).getField(0);
+    String expected = "POINT (70.71067811865476 70.71067811865474)";
+    assertEquals(expected, result.toString());
+  }
+
   @Test
   public void testReducePrecision() {
     Table polygonTable = tableEnv.sqlQuery("SELECT ST_GeomFromText('POINT(0.12 0.23)') AS geom");

python/sedona/sql/st_functions.py

@@ -1210,6 +1210,27 @@ def ST_Polygonize(geometry: ColumnOrName) -> Column:
     return _call_st_function("ST_Polygonize", (geometry))
 
 
+@validate_argument_types
+def ST_Project(
+    geom: ColumnOrName,
+    distance: Union[ColumnOrName, float],
+    azimuth: Union[ColumnOrName, float],
+    lenient: Optional[Union[ColumnOrName, bool]] = None,
+) -> Column:
+    """Calculates a new point location given a starting point, distance, and direction (azimuth).
+
+    @param geom:
+    @param distance:
+    @param azimuth:
+    @param lenient:
+    @return:
+    """
+    args = (geom, distance, azimuth, lenient)
+    if lenient is None:
+        args = (geom, distance, azimuth)
+    return _call_st_function("ST_Project", args)
+
+
 @validate_argument_types
 def ST_MakePolygon(
     line_string: ColumnOrName, holes: Optional[ColumnOrName] = None

python/tests/sql/test_dataframe_api.py

@@ -14,7 +14,7 @@
 #  KIND, either express or implied.  See the License for the
 #  specific language governing permissions and limitations
 #  under the License.
-
+from math import radians
 from typing import Callable, Tuple
 
 from pyspark.sql import functions as f, Row
@@ -755,6 +755,20 @@
         "ST_Normalize(geom)",
         "GEOMETRYCOLLECTION (POLYGON ((0 2, 1 3, 2 4, 2 3, 2 2, 1 2, 0 2)), POLYGON ((2 2, 2 3, 2 4, 3 3, 4 2, 3 2, 2 2)))",
     ),
+    (
+        stf.ST_Project,
+        ("point", 10.0, radians(10)),
+        "point_geom",
+        "",
+        "POINT (1.7364817766693021 10.848077530122081)",
+    ),
+    (
+        stf.ST_Project,
+        ("geom", 10.0, radians(10), True),
+        "triangle_geom",
+        "",
+        "POINT EMPTY",
+    ),
     (
         stf.ST_MakePolygon,
         ("geom",),
@@ -1288,6 +1302,9 @@
     (stf.ST_PointN, ("", None)),
     (stf.ST_PointN, ("", 2.0)),
     (stf.ST_PointOnSurface, (None,)),
+    (stf.ST_Project, (None, "", "", None)),
+    (stf.ST_Project, ("", None, "", None)),
+    (stf.ST_Project, ("", "", None, None)),
     (stf.ST_ReducePrecision, (None, 1)),
     (stf.ST_ReducePrecision, ("", None)),
     (stf.ST_ReducePrecision, ("", 1.0)),

python/tests/sql/test_function.py

@@ -1735,6 +1735,23 @@ def test_st_polygonize(self):
         for actual, expected in result:
             assert actual == expected
 
+    def test_st_project(self):
+        baseDf = self.spark.sql("SELECT ST_GeomFromWKT('POINT(0 0)') as point")
+        actual = baseDf.selectExpr("ST_Project(point, 10, radians(45))").first()[0].wkt
+        expected = "POINT (7.0710678118654755 7.071067811865475)"
+        assert expected == actual
+
+        actual = (
+            self.spark.sql(
+                "SELECT ST_Project(ST_MakeEnvelope(0, 1, 2, 0), 10, radians(50), true)"
+            )
+            .first()[0]
+            .wkt
+        )
+
+        expected = "POINT EMPTY"
+        assert expected == actual
+
     def test_st_make_polygon(self):
         # Given
         geometry_df = self.spark.createDataFrame(