[SEDONA-498] Add Spheroidal ST_Buffer

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

@@ -54,6 +54,7 @@
 import java.util.stream.Collectors;
 
 import static com.google.common.geometry.S2.DBL_EPSILON;
+import static org.apache.sedona.common.FunctionsGeoTools.bufferSpheroid;
 
 
 public class Functions {
@@ -87,25 +88,40 @@ public static Geometry boundary(Geometry geometry) {
     }
 
     public static Geometry buffer(Geometry geometry, double radius) {
-        return buffer(geometry, radius, "");
+        return buffer(geometry, radius, false, "");
     }
 
-    public static Geometry buffer(Geometry geometry, double radius, String params) {
-        if (params.isEmpty()) {
-            return BufferOp.bufferOp(geometry, radius);
-        }
+    public static Geometry buffer(Geometry geometry, double radius, boolean useSpheroid) {
+        return buffer(geometry, radius, useSpheroid, "");
+    }
 
-        BufferParameters bufferParameters = parseBufferParams(params);
+    public static Geometry buffer(Geometry geometry, double radius, boolean useSpheroid, String params) throws IllegalArgumentException {
+        BufferParameters bufferParameters = new BufferParameters();
 
-        // convert the sign to the appropriate direction
-        // left - radius should be positive
-        // right - radius should be negative
-        if (bufferParameters.isSingleSided() &&
-                (params.toLowerCase().contains("left") && radius < 0 || params.toLowerCase().contains("right") && radius > 0)) {
+        // Processing parameters
+        if (!params.isEmpty()) {
+            bufferParameters = parseBufferParams(params);
+
+            // convert the sign to the appropriate direction
+            // left - radius should be positive
+            // right - radius should be negative
+            if (bufferParameters.isSingleSided() &&
+                    (params.toLowerCase().contains("left") && radius < 0 || params.toLowerCase().contains("right") && radius > 0)) {
                 radius = -radius;
+            }
         }
 
-        return BufferOp.bufferOp(geometry, radius, bufferParameters);
+        if (useSpheroid) {
+            // Spheroidal buffering logic
+            try {
+                return bufferSpheroid(geometry, radius, bufferParameters);
+            } catch (RuntimeException e) {
+                throw new RuntimeException("Error processing spheroidal buffer", e);
+            }
+        } else {
+            // Existing planar buffer logic with params handling
+            return BufferOp.bufferOp(geometry, radius, bufferParameters);
+        }
     }
 
     private static BufferParameters parseBufferParams(String params) {
@@ -186,17 +202,26 @@ else if (singleParam[0].equalsIgnoreCase(listBufferParameters[5])) {
         return bufferParameters;
     }
 
-    public static int bestSRID(Geometry geometry) {
+    public static int bestSRID(Geometry geometry) throws IllegalArgumentException {
+        // Shift longitudes if geometry crosses dateline
+        if (crossesDateLine(geometry)) {
+            shiftLongitude(geometry);
+        }
+
+        // Check envelope
         Envelope envelope = geometry.getEnvelopeInternal();
         if (envelope.isNull()) return Spheroid.EPSG_WORLD_MERCATOR; // Fallback EPSG
 
         // Calculate the center of the envelope
-        double centerX =  (envelope.getMinX() + envelope.getMaxX()) / 2.0; // centroid.getX();
+        double centerX = (envelope.getMinX() + envelope.getMaxX()) / 2.0; // centroid.getX();
         double centerY = (envelope.getMinY() + envelope.getMaxY()) / 2.0; // centroid.getY();
 
         // Calculate angular width and height
-        double xwidth = Spheroid.angularWidth(envelope);
-        double ywidth = Spheroid.angularHeight(envelope);
+        Double xwidth = Spheroid.angularWidth(envelope);
+        Double ywidth = Spheroid.angularHeight(envelope);
+        if (xwidth.isNaN() | ywidth.isNaN()) {
+            throw new IllegalArgumentException("Only lon/lat coordinate systems are supported by ST_BestSRID");
+        }
 
         // Prioritize polar regions for Lambert Azimuthal Equal Area projection
         if (centerY >= 70.0 && ywidth < 45.0) return Spheroid.EPSG_NORTH_LAMBERT;
@@ -208,8 +233,8 @@ public static int bestSRID(Geometry geometry) {
             if (centerX == -180.0 || centerX == 180.0) {
                 zone = 59; // UTM zone 60
             } else {
-                zone = (int)Math.floor((centerX + 180.0) / 6.0);
-                zone = Math.min(zone, 59);
+                zone = (int) Math.floor((centerX + 180.0) / 6.0);
+                zone = (zone > 59) ? zone - 60 : zone;
             }
             return (centerY < 0.0) ? Spheroid.EPSG_SOUTH_UTM_START + zone : Spheroid.EPSG_NORTH_UTM_START + zone;
         }
@@ -218,6 +243,41 @@ public static int bestSRID(Geometry geometry) {
         return Spheroid.EPSG_WORLD_MERCATOR;
     }
 
+    /**
+     * Corrects the longitudes of a geometry to be within the -180 to 180 range.
+     * This method modifies the original geometry.
+     *
+     * @param geometry The geometry to be corrected.
+     */
+    public static void normalizeLongitude(Geometry geometry) {
+        if (geometry == null || geometry.isEmpty()) {
+            return;
+        }
+
+        geometry.apply(new CoordinateSequenceFilter() {
+            @Override
+            public void filter(CoordinateSequence seq, int i) {
+                double x = seq.getX(i);
+                // Normalize the longitude to be within -180 to 180 range
+                while (x > 180) x -= 360;
+                while (x < -180) x += 360;
+                seq.setOrdinate(i, CoordinateSequence.X, x);
+            }
+
+            @Override
+            public boolean isDone() {
+                return false; // Continue processing until all coordinates are processed
+            }
+
+            @Override
+            public boolean isGeometryChanged() {
+                return true; // The geometry is changed as we are modifying the coordinates
+            }
+        });
+
+        geometry.geometryChanged(); // Notify the geometry that its coordinates have been changed
+    }
+
     /**
      * Checks if a geometry crosses the International Date Line.
      *
@@ -263,28 +323,7 @@ public boolean isGeometryChanged() {
         return filter.isDone();
     }
 
-    public static Geometry envelope(Geometry geometry) {
-        return geometry.getEnvelope();
-    }
-
-    public static double distance(Geometry left, Geometry right) {
-        return left.distance(right);
-    }
-
-    public static double distance3d(Geometry left, Geometry right) {
-        return new Distance3DOp(left, right).distance();
-    }
-
-    public static double length(Geometry geometry) {
-        return geometry.getLength();
-    }
-
-    public static Geometry normalize(Geometry geometry) {
-        geometry.normalize();
-        return geometry;
-    }
-
-    public static Geometry shiftLongitude(Geometry geometry) {
+    public static Geometry shiftLongitude (Geometry geometry){
         geometry.apply(new CoordinateSequenceFilter() {
             @Override
             public void filter(CoordinateSequence seq, int i) {
@@ -310,6 +349,27 @@ public boolean isGeometryChanged() {
         return geometry;
     }
 
+    public static Geometry envelope (Geometry geometry){
+        return geometry.getEnvelope();
+    }
+
+    public static double distance (Geometry left, Geometry right){
+        return left.distance(right);
+    }
+
+    public static double distance3d (Geometry left, Geometry right){
+        return new Distance3DOp(left, right).distance();
+    }
+
+    public static double length (Geometry geometry){
+        return geometry.getLength();
+    }
+
+    public static Geometry normalize (Geometry geometry){
+        geometry.normalize();
+        return geometry;
+    }
+
     public static Double x(Geometry geometry) {
         if (geometry instanceof Point) {
             return geometry.getCoordinate().x;

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

@@ -20,6 +20,8 @@
 import org.locationtech.jts.geom.Envelope;
 import org.locationtech.jts.geom.Geometry;
 import org.locationtech.jts.geom.GeometryFactory;
+import org.locationtech.jts.operation.buffer.BufferOp;
+import org.locationtech.jts.operation.buffer.BufferParameters;
 import org.locationtech.jts.triangulate.VoronoiDiagramBuilder;
 import org.opengis.referencing.FactoryException;
 import org.opengis.referencing.NoSuchAuthorityCodeException;
@@ -134,4 +136,39 @@ public static Geometry voronoiPolygons(Geometry geom, double tolerance, Geometry
         }
         return builder.getDiagram(FunctionsGeoTools.GEOMETRY_FACTORY);
     }
+
+    public static Geometry bufferSpheroid(Geometry geometry, double radius, BufferParameters params) throws IllegalArgumentException {
+        // Determine the best SRID for spheroidal calculations
+        int bestCRS = Functions.bestSRID(geometry);
+        int originalCRS = geometry.getSRID();
+        final int WGS84CRS = 4326;
+
+        // Shift longitude if geometry crosses dateline
+        if (Functions.crossesDateLine(geometry)) {
+            Functions.shiftLongitude(geometry);
+        }
+//        geometry = (Predicates.crossesDateLine(geometry)) ? Functions.shiftLongitude(geometry) : geometry;
+
+        // If originalCRS is not set, use WGS84 as the originalCRS for transformation
+        String sourceCRSCode = (originalCRS == 0) ? "EPSG:" + WGS84CRS : "EPSG:" + originalCRS;
+        String targetCRSCode = "EPSG:" + bestCRS;
+
+        try {
+            // Transform the geometry to the selected SRID
+            Geometry transformedGeometry = transform(geometry, sourceCRSCode, targetCRSCode);
+            // Apply the buffer operation in the selected SRID
+            Geometry bufferedGeometry = BufferOp.bufferOp(transformedGeometry, radius, params);
+
+            // Transform back to the original SRID or to WGS 84 if original SRID was not set
+            int backTransformCRSCode = (originalCRS == 0) ? WGS84CRS : originalCRS;
+            Geometry bufferedResult = transform(bufferedGeometry, targetCRSCode, "EPSG:" + backTransformCRSCode);
+            bufferedResult.setSRID(backTransformCRSCode);
+
+            // Normalize longitudes between -180 and 180
+            Functions.normalizeLongitude(bufferedResult);
+            return bufferedResult;
+        } catch (FactoryException | TransformException e) {
+            throw new RuntimeException(e);
+        }
+    }
 }

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

@@ -16,8 +16,6 @@
 import org.locationtech.jts.geom.*;
 import org.apache.sedona.common.sphere.Spheroid;
 
-import java.util.concurrent.atomic.AtomicBoolean;
-
 public class Predicates {
     public static boolean contains(Geometry leftGeometry, Geometry rightGeometry) {
         return leftGeometry.contains(rightGeometry);

common/src/main/java/org/apache/sedona/common/sphere/Spheroid.java

@@ -26,6 +26,7 @@
 import org.locationtech.jts.geom.Envelope;
 import org.locationtech.jts.geom.Geometry;
 
+import static java.lang.Float.NaN;
 import static java.lang.Math.abs;
 
 public class Spheroid
@@ -131,7 +132,7 @@ else if (geomType.equals("MultiPolygon") || geomType.equals("GeometryCollection"
         }
     }
 
-    public static double angularWidth(Envelope envelope) {
+    public static Double angularWidth(Envelope envelope) {
         double lon1 = envelope.getMinX();
         double lon2 = envelope.getMaxX();
         double lat = (envelope.getMinY() + envelope.getMaxY()) / 2; // Mid-latitude for width calculation
@@ -141,12 +142,12 @@ public static double angularWidth(Envelope envelope) {
         double distance = g.s12; // Distance in meters
 
         // Convert distance to angular width in degrees
-        double angularWidth = Math.toDegrees(distance / (Geodesic.WGS84.EquatorialRadius() * Math.PI / 180));
+        Double angularWidth = Math.toDegrees(distance / (Geodesic.WGS84.EquatorialRadius() * Math.PI / 180));
 
         return angularWidth;
     }
 
-    public static double angularHeight(Envelope envelope) {
+    public static Double angularHeight(Envelope envelope) {
         double lat1 = envelope.getMinY();
         double lat2 = envelope.getMaxY();
         double lon = (envelope.getMinX() + envelope.getMaxX()) / 2; // Mid-longitude for height calculation
@@ -156,7 +157,7 @@ public static double angularHeight(Envelope envelope) {
         double distance = g.s12; // Distance in meters
 
         // Convert distance to angular height in degrees
-        double angularHeight = Math.toDegrees(distance / (Geodesic.WGS84.EquatorialRadius() * Math.PI / 180));
+        Double angularHeight = Math.toDegrees(distance / (Geodesic.WGS84.EquatorialRadius() * Math.PI / 180));
 
         return angularHeight;
     }