Skip to content

Commit f3dba87

Browse files
asinghvi17asinghvi17
authored
Implement a perimeter method for all manifolds (#334)
2 parents fb4513b + 5095085 commit f3dba87

File tree

9 files changed

+175
-3
lines changed

9 files changed

+175
-3
lines changed

‎Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "GeometryOps"
22
uuid = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab"
33
authors = ["Anshul Singhvi <[email protected]>", "Rafael Schouten <[email protected]>", "Skylar Gering <[email protected]>", "and contributors"]
4-
version = "0.1.28"
4+
version = "0.1.29"
55

66
[deps]
77
AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

‎ext/GeometryOpsProjExt/GeometryOpsProjExt.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,5 +7,6 @@ using GeometryOps, Proj
77

88
include("reproject.jl")
99
include("segmentize.jl")
10+
include("perimeter_area.jl")
1011

1112
end

‎ext/GeometryOpsProjExt/perimeter_area.jl‎

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,68 @@
1+
import GeometryOps: perimeter, area, applyreduce, TraitTarget, WithTrait
2+
import GeoInterface as GI
3+
4+
function perimeter(m::Geodesic, geom, ::Type{T} = Float64; init = zero(T), kwargs...) where T
5+
# Create a Proj geodesic object using the ellipsoid parameters from the Geodesic manifold
6+
proj_geodesic = Ref(Proj.geod_geodesic(m.semimajor_axis, 1/m.inv_flattening))
7+
proj_polygon = Ref(Proj._null(Proj.geod_polygon))
8+
9+
function _perimeter_geodesic_inner(trait, geom)
10+
@assert GI.npoint(geom) >= 2 "Geodesic perimeter requires at least 2 points"
11+
12+
# Initialize the polygon
13+
proj_polygon[] = Proj._null(Proj.geod_polygon)
14+
Proj.geod_polygon_init(proj_polygon, 1)
15+
16+
# Add all points to the polygon
17+
for point in GI.getpoint(trait, geom)
18+
lat, lon = GI.y(point), GI.x(point) # Proj expects lat, lon order
19+
Proj.geod_polygon_addpoint(proj_geodesic, proj_polygon, lat, lon)
20+
end
21+
22+
# Compute the polygon properties
23+
# geod_polygon_compute returns (num_vertices, perimeter, area)
24+
area_result, perimeter_result = Proj.geod_polygon_compute(proj_geodesic[], proj_polygon[], false, true)
25+
26+
return T(perimeter_result)
27+
end
28+
29+
return applyreduce(
30+
WithTrait(_perimeter_geodesic_inner),
31+
+,
32+
TraitTarget(GI.AbstractCurveTrait),
33+
geom; init, kwargs...
34+
)
35+
end
36+
37+
function GeometryOps.area(m::Geodesic, geom, ::Type{T} = Float64; threaded = False(), init = zero(T), kwargs...) where T
38+
# Create a Proj geodesic object using the ellipsoid parameters from the Geodesic manifold
39+
proj_geodesic = Ref(Proj.geod_geodesic(m.semimajor_axis, 1/m.inv_flattening))
40+
proj_polygon = Ref(Proj._null(Proj.geod_polygon))
41+
42+
function _area_geodesic_inner(trait, geom)
43+
@assert GI.npoint(geom) >= 2 "Geodesic perimeter requires at least 2 points"
44+
45+
# Initialize the polygon
46+
proj_polygon[] = Proj._null(Proj.geod_polygon)
47+
Proj.geod_polygon_init(proj_polygon, 0)
48+
49+
# Add all points to the polygon
50+
for point in GI.getpoint(trait, geom)
51+
lat, lon = GI.y(point), GI.x(point) # Proj expects lat, lon order
52+
Proj.geod_polygon_addpoint(proj_geodesic, proj_polygon, lat, lon)
53+
end
54+
55+
# Compute the polygon properties
56+
# geod_polygon_compute returns (num_vertices, perimeter, area)
57+
area_result, perimeter_result = Proj.geod_polygon_compute(proj_geodesic[], proj_polygon[], true, false)
58+
59+
return T(area_result)
60+
end
61+
62+
return applyreduce(
63+
WithTrait(_area_geodesic_inner),
64+
+,
65+
TraitTarget(GI.AbstractCurveTrait),
66+
geom; init, kwargs...
67+
)
68+
end

‎src/GeometryOps.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -64,6 +64,7 @@ include("methods/centroid.jl")
6464
include("methods/convex_hull.jl")
6565
include("methods/distance.jl")
6666
include("methods/equals.jl")
67+
include("methods/perimeter.jl")
6768
include("methods/clipping/predicates.jl")
6869
include("methods/clipping/clipping_processor.jl")
6970
include("methods/clipping/coverage.jl")

‎src/methods/area.jl‎

Lines changed: 9 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -66,8 +66,12 @@ This is computed slightly differently for different geometries:
6666
Result will be of type T, where T is an optional argument with a default value
6767
of Float64.
6868
"""
69-
function area(geom, ::Type{T} = Float64; threaded=false) where T <: AbstractFloat
70-
applyreduce(WithTrait((trait, g) -> _area(T, trait, g)), +, _AREA_TARGETS, geom; threaded, init=zero(T))
69+
function area(geom, ::Type{T} = Float64; threaded=false, kwargs...) where T <: AbstractFloat
70+
area(Planar(), geom, T; threaded, kwargs...)
71+
end
72+
73+
function area(::Planar, geom, ::Type{T} = Float64; threaded=false, kwargs...) where T <: AbstractFloat
74+
applyreduce(WithTrait((trait, g) -> _area(T, trait, g)), +, _AREA_TARGETS, geom; threaded, init=zero(T), kwargs...)
7175
end
7276

7377
"""
@@ -151,3 +155,6 @@ function _signed_area(::Type{T}, geom) where T
151155
area += _area_component(p1, p2)
152156
return T(area / 2)
153157
end
158+
159+
160+
# ## Spherical

‎src/methods/perimeter.jl‎

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
1+
#=
2+
# Perimeter
3+
4+
The perimeter of a geometry is the length of its boundary. In many contexts
5+
this is called `length`; to avoid clashing with `Base.length` in Julia, we
6+
call this `perimeter`.
7+
8+
## Examples
9+
10+
=#
11+
12+
13+
function perimeter(geom, ::Type{T} = Float64; threaded=False(), init = zero(T), kwargs...) where T
14+
perimeter(Planar(), geom, T; threaded, init, kwargs...)
15+
end
16+
17+
#=
18+
The planar implementation is straightforward.
19+
=#
20+
21+
function perimeter(::Planar, geom, ::Type{T} = Float64; init = zero(T), kwargs...) where T
22+
function _perimeter_planar_inner(trait, geom)
23+
@assert GI.npoint(geom) >= 2 "Planar perimeter requires at least 2 points"
24+
distance = zero(T)
25+
for (p1, p2) in eachedge(trait, geom, T)
26+
distance += hypot(GI.x(p2) - GI.x(p1), GI.y(p2) - GI.y(p1))
27+
end
28+
return distance
29+
end
30+
return applyreduce(
31+
WithTrait(_perimeter_planar_inner),
32+
+,
33+
TraitTarget(GI.AbstractCurveTrait),
34+
geom; init, kwargs...
35+
)
36+
end
37+
38+
using .UnitSpherical: UnitSphericalPoint
39+
40+
function perimeter(m::Spherical, geom, ::Type{T} = Float64; init = zero(T), kwargs...) where T
41+
function _perimeter_spherical_inner(trait, geom)
42+
@assert GI.npoint(geom) >= 2 "Spherical perimeter requires at least 2 points"
43+
p1_unknown, rest = Iterators.peel(GI.getpoint(trait, geom))
44+
p1 = UnitSphericalPoint(GI.PointTrait(), p1_unknown)
45+
distance = zero(T)
46+
for p2 in Iterators.map(p -> UnitSphericalPoint(GI.PointTrait(), p), rest)
47+
distance += spherical_distance(p1, p2)
48+
p1 = p2
49+
end
50+
return distance
51+
end
52+
return applyreduce(
53+
WithTrait(_perimeter_spherical_inner),
54+
+,
55+
TraitTarget(GI.AbstractCurveTrait),
56+
geom; init, kwargs...
57+
) * m.radius
58+
end
59+
60+
# The `Geodesic` implementation is in `ext/GeometryOpsProjExt/perimeter.jl`

‎test/methods/area.jl‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -97,3 +97,11 @@ end
9797
@test GO.area($mp1) == a2 + a4
9898
@test GO.area($mp1, Float32) isa Float32
9999
end
100+
101+
102+
highlat_poly = LG.Polygon([[[70., 70.], [70., 80.], [80., 80.], [80., 70.], [70., 70.]]])
103+
104+
@testset_implementations "Spherical/geodesic" begin
105+
@test GO.area(GO.Planar(), $highlat_poly) == 100
106+
@test GO.area(GO.Planar(), $highlat_poly) < GO.area(GO.Geodesic(), $highlat_poly)
107+
end

‎test/methods/perimeter.jl‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
using Test
2+
import GeometryOps as GO
3+
import LibGEOS as LG
4+
5+
# Basic geometries for testing
6+
point = LG.Point([0.0, 0.0])
7+
line = LG.LineString([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0]])
8+
square = LG.Polygon([[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.0, 0.0]]])
9+
10+
@testset "Basic perimeter tests" begin
11+
# Points have zero perimeter
12+
@test GO.perimeter(point) == 0
13+
14+
# Lines have perimeter equal to their length
15+
@test GO.perimeter(line) == 2.0 # 1.0 + 1.0
16+
17+
# Square has perimeter of 4 (each side is 1)
18+
@test GO.perimeter(square) == 4.0
19+
end
20+
21+
@testset "Spherical and geodesic" begin
22+
highlat_poly = LG.Polygon([[[70., 70.], [70., 80.], [80., 80.], [80., 70.], [70., 70.]]])
23+
@test GO.perimeter(GO.Planar(), highlat_poly) == 40
24+
@test GO.perimeter(GO.Planar(), highlat_poly) < GO.perimeter(GO.Spherical(), highlat_poly)
25+
@test GO.perimeter(GO.Spherical(), highlat_poly) < GO.perimeter(GO.Geodesic(), highlat_poly)
26+
end

‎test/runtests.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@ end
2222
# Methods
2323
@safetestset "Angles" begin include("methods/angles.jl") end
2424
@safetestset "Area" begin include("methods/area.jl") end
25+
@safetestset "Perimeter" begin include("methods/perimeter.jl") end
2526
@safetestset "Barycentric coordinate operations" begin include("methods/barycentric.jl") end
2627
@safetestset "Orientation" begin include("methods/orientation.jl") end
2728
@safetestset "Centroid" begin include("methods/centroid.jl") end

0 commit comments

Comments
 (0)