Implement geodesic area

Author: asinghvi17

ext/GeometryOpsProjExt/GeometryOpsProjExt.jl

@@ -7,6 +7,6 @@ using GeometryOps, Proj
 
 include("reproject.jl")
 include("segmentize.jl")
-include("perimeter.jl")
+include("perimeter_area.jl")
 
 end

ext/GeometryOpsProjExt/perimeter.jl

@@ -0,0 +1,68 @@
+import GeometryOps: perimeter, area, applyreduce, TraitTarget, WithTrait
+import GeoInterface as GI
+
+function perimeter(m::Geodesic, geom, ::Type{T} = Float64; init = zero(T), kwargs...) where T
+    # Create a Proj geodesic object using the ellipsoid parameters from the Geodesic manifold
+    proj_geodesic = Ref(Proj.geod_geodesic(m.semimajor_axis, 1/m.inv_flattening))
+    proj_polygon = Ref(Proj._null(Proj.geod_polygon))
+    
+    function _perimeter_geodesic_inner(trait, geom)
+        @assert GI.npoint(geom) >= 2 "Geodesic perimeter requires at least 2 points"
+        
+        # Initialize the polygon
+        proj_polygon[] = Proj._null(Proj.geod_polygon)
+        Proj.geod_polygon_init(proj_polygon, 1)
+        
+        # Add all points to the polygon
+        for point in GI.getpoint(trait, geom)
+            lat, lon = GI.y(point), GI.x(point)  # Proj expects lat, lon order
+            Proj.geod_polygon_addpoint(proj_geodesic, proj_polygon, lat, lon)
+        end
+        
+        # Compute the polygon properties
+        # geod_polygon_compute returns (num_vertices, perimeter, area)
+        area_result, perimeter_result = Proj.geod_polygon_compute(proj_geodesic[], proj_polygon[], false, true)
+        
+        return T(perimeter_result)
+    end
+    
+    return applyreduce(
+        WithTrait(_perimeter_geodesic_inner), 
+        +, 
+        TraitTarget(GI.AbstractCurveTrait), 
+        geom; init, kwargs...
+    )
+end
+
+function GeometryOps.area(m::Geodesic, geom, ::Type{T} = Float64; threaded = False(), init = zero(T), kwargs...) where T
+     # Create a Proj geodesic object using the ellipsoid parameters from the Geodesic manifold
+     proj_geodesic = Ref(Proj.geod_geodesic(m.semimajor_axis, 1/m.inv_flattening))
+     proj_polygon = Ref(Proj._null(Proj.geod_polygon))
+     
+     function _area_geodesic_inner(trait, geom)
+         @assert GI.npoint(geom) >= 2 "Geodesic perimeter requires at least 2 points"
+         
+         # Initialize the polygon
+         proj_polygon[] = Proj._null(Proj.geod_polygon)
+         Proj.geod_polygon_init(proj_polygon, 0)
+         
+         # Add all points to the polygon
+         for point in GI.getpoint(trait, geom)
+             lat, lon = GI.y(point), GI.x(point)  # Proj expects lat, lon order
+             Proj.geod_polygon_addpoint(proj_geodesic, proj_polygon, lat, lon)
+         end
+         
+         # Compute the polygon properties
+         # geod_polygon_compute returns (num_vertices, perimeter, area)
+         area_result, perimeter_result = Proj.geod_polygon_compute(proj_geodesic[], proj_polygon[], true, false)
+         
+         return T(area_result)
+     end
+     
+     return applyreduce(
+         WithTrait(_area_geodesic_inner), 
+         +, 
+         TraitTarget(GI.AbstractCurveTrait), 
+         geom; init, kwargs...
+     )
+end

ext/GeometryOpsProjExt/perimeter_area.jl

@@ -66,8 +66,12 @@ This is computed slightly differently for different geometries:
 Result will be of type T, where T is an optional argument with a default value
 of Float64.
 """
-function area(geom, ::Type{T} = Float64; threaded=false) where T <: AbstractFloat
-    applyreduce(WithTrait((trait, g) -> _area(T, trait, g)), +, _AREA_TARGETS, geom; threaded, init=zero(T))
+function area(geom, ::Type{T} = Float64; threaded=false, kwargs...) where T <: AbstractFloat
+    area(Planar(), geom, T; threaded, kwargs...)
+end
+
+function area(::Planar, geom, ::Type{T} = Float64; threaded=false, kwargs...) where T <: AbstractFloat
+    applyreduce(WithTrait((trait, g) -> _area(T, trait, g)), +, _AREA_TARGETS, geom; threaded, init=zero(T), kwargs...)
 end
 
 """
@@ -151,3 +155,6 @@ function _signed_area(::Type{T}, geom) where T
     area += _area_component(p1, p2)
     return T(area / 2)
 end
+
+
+# ## Spherical