From 957a46df397df635110ff6ab11b254f687d6b8b5 Mon Sep 17 00:00:00 2001 From: Elias Carvalho <73039601+eliascarv@users.noreply.github.com> Date: Wed, 6 Mar 2024 17:07:56 -0300 Subject: [PATCH] Add tests for projection domain (#17) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * 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 --------- Co-authored-by: Júlio Hoffimann --- src/Cartography.jl | 1 + src/crs/projected.jl | 16 ++- src/crs/projected/mercator.jl | 2 +- src/crs/projected/orthographic.jl | 13 +- src/crs/projected/webmercator.jl | 2 +- test/crs.jl | 200 ++++++++++++++++++++++++++++++ test/runtests.jl | 2 + 7 files changed, 226 insertions(+), 10 deletions(-) diff --git a/src/Cartography.jl b/src/Cartography.jl index e4834f7..bc2ec0c 100644 --- 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 index 674129b..7983a3e 100644 --- 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 index ac956c6..c12d6ab 100644 --- 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 index c56c523..d63b78b 100644 --- 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 index c8b2e24..e312054 100644 --- 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 index 0a41188..283096b 100644 --- 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 index 5f21599..6e0cafa 100644 --- 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"] # --------------------------------