Add method for measure of Segment with Ellipsoid manifold (#1069)

* Add method for 'measure' of 'Segment' with Ellipsoid manifold

* Apply suggestions from code review

Co-authored-by: Júlio Hoffimann <[email protected]>

---------

Co-authored-by: Júlio Hoffimann <[email protected]>

Author: eliascarv

src/discretization/maxlength.jl

@@ -20,7 +20,7 @@ function discretize(box::Box, method::MaxLengthDiscretization)
 end
 
 function discretize(segment::Segment, method::MaxLengthDiscretization)
-  size = ceil(Int, _measure(segment) / method.length)
+  size = ceil(Int, measure(segment) / method.length)
   discretize(segment, RegularDiscretization(size))
 end
 
@@ -41,18 +41,5 @@ function _sides(box::Box{<:🌐})
 
   AP = Segment(A, P)
   PB = Segment(P, B)
-  (_measure(AP), _measure(PB))
-end
-
-_measure(segment::Segment{<:𝔼}) = measure(segment)
-
-# TODO: Haversine returns the shortest distance between two points
-# this is not always equal to the distance between two directed points
-function _measure(segment::Segment{<:🌐})
-  T = numtype(lentype(segment))
-  🌎 = ellipsoid(datum(crs(segment)))
-  r = numconvert(T, majoraxis(🌎))
-
-  A, B = extrema(segment)
-  evaluate(Haversine(r), A, B)
+  (measure(AP), measure(PB))
 end

src/geometries/primitives/box.jl

@@ -51,9 +51,9 @@ Base.maximum(b::Box) = b.max
 
 Base.extrema(b::Box) = b.min, b.max
 
-diagonal(b::Box) = norm(b.max - b.min)
+diagonal(b::Box{<:𝔼}) = norm(b.max - b.min)
 
-sides(b::Box) = Tuple(b.max - b.min)
+sides(b::Box{<:𝔼}) = Tuple(b.max - b.min)
 
 ==(b₁::Box, b₂::Box) = b₁.min == b₂.min && b₁.max == b₂.max
 

src/measures.jl

@@ -75,7 +75,17 @@ function measure(t::Torus)
   4 * T(π)^2 * R * r
 end
 
-measure(s::Segment) = norm(maximum(s) - minimum(s))
+measure(s::Segment{<:𝔼}) = norm(maximum(s) - minimum(s))
+
+# TODO: replace Haversine by an appropriate geodesic distance
+# that considers the west-east orientation of segments
+function measure(s::Segment{<:🌐})
+  T = numtype(lentype(s))
+  🌎 = ellipsoid(datum(crs(s)))
+  r = numconvert(T, majoraxis(🌎))
+
+  evaluate(Haversine(r), extrema(s)...)
+end
 
 function measure(t::Triangle)
   A, B, C = vertices(t)

test/polytopes.jl

@@ -71,6 +71,17 @@
   s = Segment(merc(0, 0), merc(1, 1))
   @test crs(s(T(0))) === crs(s)
 
+  # measure
+  s = Segment(merc(0, 0), merc(1, 1))
+  @test measure(s) ≈ T(√2) * u"m"
+  s = Segment(latlon(0, 45), latlon(0, 135))
+  r = majoraxis(ellipsoid(datum(crs(s))))
+  C = 2 * T(π) * ℳ(r)
+  @test measure(s) ≈ C / 4
+  # TODO: fix measure of segments on the globe
+  s = Segment(latlon(0, 135), latlon(0, 45))
+  @test_broken measure(s) ≈ 3C / 4
+
   s = Segment(cart(0, 0), cart(1, 1))
   @test sprint(show, s) == "Segment((x: 0.0 m, y: 0.0 m), (x: 1.0 m, y: 1.0 m))"
   if T === Float32