Add tests for projection domain (#17)

* Add tests for projection range * Apply suggestions * Fix typo * Update tests * Add tests for 'OrthoNorth' and 'OrthoSouth' * Apply suggestions * Apply suggestions from code review Co-authored-by: Júlio Hoffimann <[email protected]> --------- Co-authored-by: Júlio Hoffimann <[email protected]>
JuliaEarth · Mar 6, 2024 · 957a46d · 957a46d
1 parent 0a648f8
commit 957a46d
Show file tree

Hide file tree

Showing 7 changed files with 226 additions and 10 deletions.
diff --git a/src/Cartography.jl b/src/Cartography.jl
@@ -64,6 +64,7 @@ export
   OrthoNorth,
   OrthoSouth,
   datum,
+  indomain,
 
   # codes
   EPSG,

diff --git a/src/crs/projected.jl b/src/crs/projected.jl
@@ -17,7 +17,19 @@ with `f(λ::T, ϕ::T) -> T` for both functions.
 """
 function formulas end
 
-inrange(::Type{<:Projected}, λ, ϕ) = -π ≤ λ ≤ π && -π / 2 ≤ ϕ ≤ π / 2
+"""
+    inbounds(CRS::Type{<:Projected}, λ, ϕ)
+
+Checks whether `λ` and `ϕ` in radians are within the `CRS` domain.
+"""
+inbounds(::Type{<:Projected}, λ, ϕ) = -π ≤ λ ≤ π && -π / 2 ≤ ϕ ≤ π / 2
+
+"""
+    indomain(CRS::Type{<:Projected}, latlon::LatLon)
+
+Checks whether `latlon` coordinates are within the `CRS` domain.
+"""
+indomain(C::Type{<:Projected}, (; lat, lon)::LatLon) = inbounds(C, ustrip(deg2rad(lon)), ustrip(deg2rad(lat)))
 
 # ----------------
 # IMPLEMENTATIONS
@@ -39,7 +51,7 @@ function Base.convert(::Type{C}, coords::LatLon{Datum}) where {Datum,C<:Projecte
   T = numtype(coords.lon)
   λ = ustrip(deg2rad(coords.lon))
   ϕ = ustrip(deg2rad(coords.lat))
-  if !inrange(C, λ, ϕ)
+  if !inbounds(C, λ, ϕ)
     throw(ArgumentError("coordinates outside of the projection domain"))
   end
   a = numconvert(T, majoraxis(ellipsoid(Datum)))

diff --git a/src/crs/projected/mercator.jl b/src/crs/projected/mercator.jl
@@ -37,7 +37,7 @@ Mercator(args...) = Mercator{WGS84Latest}(args...)
 # CONVERSIONS
 # ------------
 
-inrange(::Type{<:Mercator}, λ, ϕ) = -π ≤ λ ≤ π && -deg2rad(80) ≤ ϕ ≤ deg2rad(84)
+inbounds(::Type{<:Mercator}, λ, ϕ) = -π ≤ λ ≤ π && -deg2rad(80) ≤ ϕ ≤ deg2rad(84)
 
 function formulas(::Type{<:Mercator{Datum}}, ::Type{T}) where {Datum,T}
   e = T(eccentricity(ellipsoid(Datum)))

diff --git a/src/crs/projected/orthographic.jl b/src/crs/projected/orthographic.jl
@@ -74,16 +74,16 @@ const OrthoSouth{Datum} = Orthographic{-90.0u"°",0.0u"°",false,Datum}
 #                     https://mathworld.wolfram.com/OrthographicProjection.html
 #                     https://epsg.org/guidance-notes.html
 
-function inrange(::Type{<:Orthographic{lat₀,lon₀}}, λ, ϕ) where {lat₀,lon₀}
+function inbounds(::Type{<:Orthographic{lat₀,lon₀}}, λ, ϕ) where {lat₀,lon₀}
   λ₀ = ustrip(deg2rad(lon₀))
   ϕ₀ = ustrip(deg2rad(lat₀))
   c = acos(sin(ϕ₀) * sin(ϕ) + cos(ϕ₀) * cos(ϕ) * cos(λ - λ₀))
   -π / 2 < c < π / 2
 end
 
-inrange(::Type{<:OrthoNorth}, λ, ϕ) = -π ≤ λ ≤ π && 0 ≤ ϕ ≤ π / 2
+inbounds(::Type{<:OrthoNorth}, λ, ϕ) = -π ≤ λ ≤ π && 0 ≤ ϕ ≤ π / 2
 
-inrange(::Type{<:OrthoSouth}, λ, ϕ) = -π ≤ λ ≤ π && -π / 2 ≤ ϕ ≤ 0
+inbounds(::Type{<:OrthoSouth}, λ, ϕ) = -π ≤ λ ≤ π && -π / 2 ≤ ϕ ≤ 0
 
 function formulas(::Type{<:Orthographic{lat₀,lon₀,false,Datum}}, ::Type{T}) where {lat₀,lon₀,Datum,T}
   λ₀ = T(ustrip(deg2rad(lon₀)))
@@ -119,17 +119,18 @@ function formulas(::Type{<:Orthographic{lat₀,lon₀,true,Datum}}, ::Type{T}) w
 end
 
 function sphericalinv(x, y, λ₀, ϕ₀)
+  fix(x) = clamp(x, -one(x), one(x))
   ρ = hypot(x, y)
   if ρ < atol(x)
     λ₀, ϕ₀
   else
-    c = asin(ρ)
+    c = asin(fix(ρ))
     sinc = sin(c)
     cosc = cos(c)
     sinϕ₀ = sin(ϕ₀)
     cosϕ₀ = cos(ϕ₀)
-    λ = λ₀ + atan(x * sinc / (ρ * cosϕ₀ * cosc - y * sinϕ₀ * sinc))
-    ϕ = asin(cosc * sinϕ₀ + (y * sinc * cosϕ₀ / ρ))
+    λ = λ₀ + atan(x * sinc, ρ * cosϕ₀ * cosc - y * sinϕ₀ * sinc)
+    ϕ = asin(fix(cosc * sinϕ₀ + (y * sinc * cosϕ₀ / ρ)))
     λ, ϕ
   end
 end

diff --git a/src/crs/projected/webmercator.jl b/src/crs/projected/webmercator.jl
@@ -37,7 +37,7 @@ WebMercator(args...) = WebMercator{WGS84Latest}(args...)
 # CONVERSIONS
 # ------------
 
-function inrange(::Type{<:WebMercator}, λ, ϕ)
+function inbounds(::Type{<:WebMercator}, λ, ϕ)
   θ = deg2rad(85.06)
   -π ≤ λ ≤ π && -θ ≤ ϕ ≤ θ
 end

diff --git a/test/crs.jl b/test/crs.jl
@@ -1277,4 +1277,204 @@
       @inferred convert(WinkelTripel{ITRFLatest}, c2)
     end
   end
+
+  @testset "Projection domain" begin
+    @testset "Mercator" begin
+      for lat in T.(-90:90), lon in T.(-180:180)
+        c1 = LatLon(lat, lon)
+        if indomain(Mercator, c1)
+          c2 = convert(Mercator{WGS84Latest}, c1)
+          @test isnum(c2.x)
+          @test isnum(c2.y)
+          c3 = convert(LatLon{WGS84Latest}, c2)
+          @test c3 ≈ c1
+        else
+          @test_throws ArgumentError convert(Mercator{WGS84Latest}, c1)
+        end
+      end
+    end
+
+    @testset "WebMercator" begin
+      for lat in T.(-90:90), lon in T.(-180:180)
+        c1 = LatLon(lat, lon)
+        if indomain(WebMercator, c1)
+          c2 = convert(WebMercator{WGS84Latest}, c1)
+          @test isnum(c2.x)
+          @test isnum(c2.y)
+          c3 = convert(LatLon{WGS84Latest}, c2)
+          @test c3 ≈ c1
+        else
+          @test_throws ArgumentError convert(WebMercator{WGS84Latest}, c1)
+        end
+      end
+    end
+
+    @testset "PlateCarree" begin
+      for lat in T.(-90:90), lon in T.(-180:180)
+        c1 = LatLon(lat, lon)
+        if indomain(PlateCarree, c1)
+          c2 = convert(PlateCarree{WGS84Latest}, c1)
+          @test isnum(c2.x)
+          @test isnum(c2.y)
+          c3 = convert(LatLon{WGS84Latest}, c2)
+          @test c3 ≈ c1
+        else
+          @test_throws ArgumentError convert(PlateCarree{WGS84Latest}, c1)
+        end
+      end
+    end
+
+    @testset "Lambert" begin
+      atol = T === Float32 ? 1.0f-2u"°" : 1e-4u"°"
+      for lat in T.(-90:90), lon in T.(-180:180)
+        c1 = LatLon(lat, lon)
+        if indomain(Lambert, c1)
+          c2 = convert(Lambert{WGS84Latest}, c1)
+          @test isnum(c2.x)
+          @test isnum(c2.y)
+          c3 = convert(LatLon{WGS84Latest}, c2)
+          @test isapprox(c3, c1; atol)
+        else
+          @test_throws ArgumentError convert(Lambert{WGS84Latest}, c1)
+        end
+      end
+    end
+
+    @testset "Behrmann" begin
+      atol = T === Float32 ? 1.0f-2u"°" : 1e-4u"°"
+      for lat in T.(-90:90), lon in T.(-180:180)
+        c1 = LatLon(lat, lon)
+        if indomain(Behrmann, c1)
+          c2 = convert(Behrmann{WGS84Latest}, c1)
+          @test isnum(c2.x)
+          @test isnum(c2.y)
+          c3 = convert(LatLon{WGS84Latest}, c2)
+          @test isapprox(c3, c1; atol)
+        else
+          @test_throws ArgumentError convert(Behrmann{WGS84Latest}, c1)
+        end
+      end
+    end
+
+    @testset "GallPeters" begin
+      atol = T === Float32 ? 1.0f-2u"°" : 1e-4u"°"
+      for lat in T.(-90:90), lon in T.(-180:180)
+        c1 = LatLon(lat, lon)
+        if indomain(GallPeters, c1)
+          c2 = convert(GallPeters{WGS84Latest}, c1)
+          @test isnum(c2.x)
+          @test isnum(c2.y)
+          c3 = convert(LatLon{WGS84Latest}, c2)
+          @test isapprox(c3, c1; atol)
+        else
+          @test_throws ArgumentError convert(GallPeters{WGS84Latest}, c1)
+        end
+      end
+    end
+
+    @testset "WinkelTripel" begin
+      for lat in T.(-90:90), lon in T.(-180:180)
+        c1 = LatLon(lat, lon)
+        if indomain(WinkelTripel, c1)
+          c2 = convert(WinkelTripel{WGS84Latest}, c1)
+          @test isnum(c2.x)
+          @test isnum(c2.y)
+          c3 = convert(LatLon{WGS84Latest}, c2)
+          @test c3 ≈ c1
+        else
+          @test_throws ArgumentError convert(WinkelTripel{WGS84Latest}, c1)
+        end
+      end
+    end
+
+    @testset "Robinson" begin
+      atol = T(1e-3) * u"°"
+      for lat in T.(-90:90), lon in T.(-180:180)
+        c1 = LatLon(lat, lon)
+        if indomain(Robinson, c1)
+          c2 = convert(Robinson{WGS84Latest}, c1)
+          @test isnum(c2.x)
+          @test isnum(c2.y)
+          c3 = convert(LatLon{WGS84Latest}, c2)
+          @test isapprox(c3, c1; atol)
+        else
+          @test_throws ArgumentError convert(Robinson{WGS84Latest}, c1)
+        end
+      end
+    end
+
+    @testset "OrthoNorth forward" begin
+      for lat in T.(-90:90), lon in T.(-180:180)
+        c1 = LatLon(lat, lon)
+        if indomain(OrthoNorth, c1)
+          c2 = convert(OrthoNorth{WGS84Latest}, c1)
+          @test isnum(c2.x)
+          @test isnum(c2.y)
+        else
+          @test_throws ArgumentError convert(OrthoNorth{WGS84Latest}, c1)
+        end
+      end
+    end
+
+    @testset "OrthoNorth inverse" begin
+      # coordinates at the singularity of the projection (lat ≈ 90) cannot be inverted
+      for lat in T.(1:89), lon in T.(-180:180)
+        c1 = LatLon(lat, lon)
+        if indomain(OrthoNorth, c1)
+          c2 = convert(OrthoNorth{WGS84Latest}, c1)
+          c3 = convert(LatLon{WGS84Latest}, c2)
+          @test c3 ≈ c1
+        end
+      end
+
+      # coordinates at the edge of the projection (lat ≈ 0)
+      # cannot be accurately inverted by numerical problems
+      atol = T(0.5) * u"°"
+      for lon in T.(-180:180)
+        c1 = LatLon(T(0), lon)
+        if indomain(OrthoNorth, c1)
+          c2 = convert(OrthoNorth{WGS84Latest}, c1)
+          c3 = convert(LatLon{WGS84Latest}, c2)
+          @test isapprox(c3, c1; atol)
+        end
+      end
+    end
+
+    @testset "OrthoSouth forward" begin
+      for lat in T.(-90:90), lon in T.(-180:180)
+        c1 = LatLon(lat, lon)
+        if indomain(OrthoSouth, c1)
+          c2 = convert(OrthoSouth{WGS84Latest}, c1)
+          @test isnum(c2.x)
+          @test isnum(c2.y)
+        else
+          @test_throws ArgumentError convert(OrthoSouth{WGS84Latest}, c1)
+        end
+      end
+    end
+
+    @testset "OrthoSouth inverse" begin
+      # coordinates at the singularity of the projection (lat ≈ -90) cannot be inverted
+      for lat in T.(-89:-1), lon in T.(-180:180)
+        c1 = LatLon(lat, lon)
+        if indomain(OrthoSouth, c1)
+          c2 = convert(OrthoSouth{WGS84Latest}, c1)
+          c3 = convert(LatLon{WGS84Latest}, c2)
+          @test c3 ≈ c1
+        end
+      end
+
+      # coordinates at the edge of the projection (lat ≈ 0)
+      # cannot be accurately inverted by numerical problems
+      atol = T(0.5) * u"°"
+      for lon in T.(-180:180)
+        c1 = LatLon(T(0), lon)
+        if indomain(OrthoSouth, c1)
+          c2 = convert(OrthoSouth{WGS84Latest}, c1)
+          c3 = convert(LatLon{WGS84Latest}, c2)
+          @test isapprox(c3, c1; atol)
+        end
+      end
+    end
+  end
 end
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -2,6 +2,8 @@ using Cartography
 using Unitful
 using Test
 
+isnum(x) = !isnan(x) && !isinf(x)
+
 testfiles = ["datums.jl", "crs.jl", "codes.jl"]
 
 # --------------------------------