diff --git a/src/crs/projected/albers.jl b/src/crs/projected/albers.jl index 44c7594..a21e844 100644 --- a/src/crs/projected/albers.jl +++ b/src/crs/projected/albers.jl @@ -98,7 +98,7 @@ function backward(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, x, y) where {lat θ = n ≥ 0 ? atan(x, ρₒ - y) : atan(-x, y - ρₒ) ρ′ = sqrt(x^2 + (ρₒ - y)^2) α′ = (C - (ρ′^2 * n^2)) / n - β′ = asin(α′ / (1 - ((1 - e²) / (2 * e)) * log((1 - e) / (1 + e)))) + β′ = asinclamp(α′ / (1 - ((1 - e²) / (2 * e)) * log((1 - e) / (1 + e)))) λ = θ / n ϕ = auth2geod(β′, e²) diff --git a/src/crs/projected/orthographic.jl b/src/crs/projected/orthographic.jl index 8826064..f91343b 100644 --- a/src/crs/projected/orthographic.jl +++ b/src/crs/projected/orthographic.jl @@ -160,18 +160,17 @@ Base.convert(::Type{Orthographic{Mode,latₒ}}, coords::CRS{Datum}) where {Mode, # ----------------- function sphericalinv(x, y, λₒ, ϕₒ) - fix(x) = clamp(x, -one(x), one(x)) ρ = hypot(x, y) if ρ < atol(x) λₒ, ϕₒ else - c = asin(fix(ρ)) + c = asinclamp(ρ) sinc = sin(c) cosc = cos(c) sinϕₒ = sin(ϕₒ) cosϕₒ = cos(ϕₒ) λ = atan(x * sinc, ρ * cosϕₒ * cosc - y * sinϕₒ * sinc) - ϕ = asin(fix(cosc * sinϕₒ + (y * sinc * cosϕₒ / ρ))) + ϕ = asinclamp(cosc * sinϕₒ + (y * sinc * cosϕₒ / ρ)) λ, ϕ end end diff --git a/src/utils.jl b/src/utils.jl index be037f6..e2f04cd 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -64,6 +64,13 @@ function atanpos(y, x) ifelse(α ≥ zero(α), α, α + oftype(α, 2π)) end +""" + asinclamp(x) + +Adjust `x` to the range `[-1,1]` to handle floating-point errors. +""" +asinclamp(x) = asin(clamp(x, -one(x), one(x))) + """ checklat(lat)