diff --git a/src/CoordRefSystems.jl b/src/CoordRefSystems.jl index 93eb46c..681b3da 100644 --- a/src/CoordRefSystems.jl +++ b/src/CoordRefSystems.jl @@ -89,6 +89,7 @@ export OrthoNorth, OrthoSouth, TransverseMercator, + Albers, datum, indomain, diff --git a/src/crs/projected.jl b/src/crs/projected.jl index 88c59c0..7f47a25 100644 --- a/src/crs/projected.jl +++ b/src/crs/projected.jl @@ -111,6 +111,7 @@ include("projected/winkeltripel.jl") include("projected/robinson.jl") include("projected/orthographic.jl") include("projected/transversemercator.jl") +include("projected/albers.jl") # ---------- # FALLBACKS @@ -141,7 +142,7 @@ function Base.convert(::Type{C}, coords::LatLon{Datum}) where {Datum,C<:Projecte xₒ = numconvert(T, S.xₒ) yₒ = numconvert(T, S.yₒ) lonₒ = numconvert(T, S.lonₒ) - λ = ustrip(deg2rad(coords.lon - lonₒ)) + λ = ustrip(deg2rad(fixlon(coords.lon - lonₒ))) ϕ = ustrip(deg2rad(coords.lat)) if !inbounds(C, λ, ϕ) throw(ArgumentError("coordinates outside of the projection domain")) diff --git a/src/crs/projected/albers.jl b/src/crs/projected/albers.jl new file mode 100644 index 0000000..a21e844 --- /dev/null +++ b/src/crs/projected/albers.jl @@ -0,0 +1,138 @@ +# ------------------------------------------------------------------ +# Licensed under the MIT License. See LICENSE in the project root. +# ------------------------------------------------------------------ + +""" + Albers{latₒ,lat₁,lat₂,Datum,Shift} + +Albers Conic Equal Area CRS with latitude origin `latₒ` standard parallels `lat₁` and `lat₂`, `Datum` and `Shift`. +""" +struct Albers{latₒ,lat₁,lat₂,Datum,Shift,M<:Met} <: Projected{Datum,Shift} + x::M + y::M +end + +Albers{latₒ,lat₁,lat₂,Datum,Shift}(x::M, y::M) where {latₒ,lat₁,lat₂,Datum,Shift,M<:Met} = + Albers{latₒ,lat₁,lat₂,Datum,Shift,float(M)}(x, y) +Albers{latₒ,lat₁,lat₂,Datum,Shift}(x::Met, y::Met) where {latₒ,lat₁,lat₂,Datum,Shift} = + Albers{latₒ,lat₁,lat₂,Datum,Shift}(promote(x, y)...) +Albers{latₒ,lat₁,lat₂,Datum,Shift}(x::Len, y::Len) where {latₒ,lat₁,lat₂,Datum,Shift} = + Albers{latₒ,lat₁,lat₂,Datum,Shift}(uconvert(m, x), uconvert(m, y)) +Albers{latₒ,lat₁,lat₂,Datum,Shift}(x::Number, y::Number) where {latₒ,lat₁,lat₂,Datum,Shift} = + Albers{latₒ,lat₁,lat₂,Datum,Shift}(addunit(x, m), addunit(y, m)) + +Albers{latₒ,lat₁,lat₂,Datum}(args...) where {latₒ,lat₁,lat₂,Datum} = Albers{latₒ,lat₁,lat₂,Datum,Shift()}(args...) + +Albers{latₒ,lat₁,lat₂}(args...) where {latₒ,lat₁,lat₂} = Albers{latₒ,lat₁,lat₂,WGS84Latest}(args...) + +Base.convert( + ::Type{Albers{latₒ,lat₁,lat₂,Datum,Shift,M}}, + coords::Albers{latₒ,lat₁,lat₂,Datum,Shift} +) where {latₒ,lat₁,lat₂,Datum,Shift,M} = Albers{latₒ,lat₁,lat₂,Datum,Shift,M}(coords.x, coords.y) + +constructor(::Type{<:Albers{latₒ,lat₁,lat₂,Datum,Shift}}) where {latₒ,lat₁,lat₂,Datum,Shift} = + Albers{latₒ,lat₁,lat₂,Datum,Shift} + +lentype(::Type{<:Albers{latₒ,lat₁,lat₂,Datum,Shift,M}}) where {latₒ,lat₁,lat₂,Datum,Shift,M} = M + +==( + coords₁::Albers{latₒ,lat₁,lat₂,Datum,Shift}, + coords₂::Albers{latₒ,lat₁,lat₂,Datum,Shift} +) where {latₒ,lat₁,lat₂,Datum,Shift} = coords₁.x == coords₂.x && coords₁.y == coords₂.y + +# ------------ +# CONVERSIONS +# ------------ + +# Adapted from PROJ coordinate transformation software +# Initial PROJ 4.3 public domain code was put as Frank Warmerdam as copyright +# holder, but he didn't mean to imply he did the work. Essentially all work was +# done by Gerald Evenden. +# Authors of the original algorithm: Gerald Evenden and Thomas Knudsen +# reference code: https://github.com/OSGeo/PROJ/blob/master/src/projections/aea.cpp + +function inbounds(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, λ, ϕ) where {latₒ,lat₁,lat₂,Datum} + 🌎 = ellipsoid(Datum) + e = oftype(λ, eccentricity(🌎)) + e² = oftype(λ, eccentricity²(🌎)) + ϕ₁ = oftype(λ, ustrip(deg2rad(lat₁))) + ϕ₂ = oftype(λ, ustrip(deg2rad(lat₂))) + + C, n = _ambersCn(ϕ₁, ϕ₂, e, e²) + ρ = _ambersρ(ϕ, C, n, e, e²) + ρ ≥ 0 +end + +function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where {latₒ,lat₁,lat₂,Datum,T} + 🌎 = ellipsoid(Datum) + e = T(eccentricity(🌎)) + e² = T(eccentricity²(🌎)) + ϕₒ = T(ustrip(deg2rad(latₒ))) + ϕ₁ = T(ustrip(deg2rad(lat₁))) + ϕ₂ = T(ustrip(deg2rad(lat₂))) + + C, n = _ambersCn(ϕ₁, ϕ₂, e, e²) + + θ(λ) = n * λ + ρ(ϕ) = _ambersρ(ϕ, C, n, e, e²) + ρₒ = ρ(ϕₒ) + + fx(λ, ϕ) = ρ(ϕ) * sin(θ(λ)) + + fy(λ, ϕ) = ρₒ - ρ(ϕ) * cos(θ(λ)) + + fx, fy +end + +function backward(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, x, y) where {latₒ,lat₁,lat₂,Datum} + 🌎 = ellipsoid(Datum) + e = oftype(x, eccentricity(🌎)) + e² = oftype(x, eccentricity²(🌎)) + ϕₒ = oftype(x, ustrip(deg2rad(latₒ))) + ϕ₁ = oftype(x, ustrip(deg2rad(lat₁))) + ϕ₂ = oftype(x, ustrip(deg2rad(lat₂))) + + C, n = _ambersCn(ϕ₁, ϕ₂, e, e²) + ρₒ = _ambersρ(ϕₒ, C, n, e, e²) + + θ = n ≥ 0 ? atan(x, ρₒ - y) : atan(-x, y - ρₒ) + ρ′ = sqrt(x^2 + (ρₒ - y)^2) + α′ = (C - (ρ′^2 * n^2)) / n + β′ = asinclamp(α′ / (1 - ((1 - e²) / (2 * e)) * log((1 - e) / (1 + e)))) + + λ = θ / n + ϕ = auth2geod(β′, e²) + + λ, ϕ +end + +# ----------------- +# HELPER FUNCTIONS +# ----------------- + +function _ambersCn(ϕ₁, ϕ₂, e, e²) + m₁ = _albersm(ϕ₁, e²) + m₂ = _albersm(ϕ₂, e²) + α₁ = _albersα(ϕ₁, e, e²) + α₂ = _albersα(ϕ₂, e, e²) + n = (m₁^2 - m₂^2) / (α₂ - α₁) + C = m₁^2 + n * α₁ + C, n +end + +function _ambersρ(ϕ, C, n, e, e²) + α = _albersα(ϕ, e, e²) + sqrt(C - n * α) / n +end + +_albersm(ϕ, e²) = cos(ϕ) / sqrt(1 - e² * sin(ϕ)^2) + +_albersα(ϕ, e, e²) = + (1 - e²) * ((sin(ϕ) / (1 - e² * sin(ϕ)^2)) - (1 / (2 * e) * log((1 - e * sin(ϕ)) / (1 + e * sin(ϕ))))) + +# ---------- +# FALLBACKS +# ---------- + +Base.convert(::Type{Albers{latₒ,lat₁,lat₂}}, coords::CRS{Datum}) where {latₒ,lat₁,lat₂,Datum} = + convert(Albers{latₒ,lat₁,lat₂,Datum}, coords) 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/get.jl b/src/get.jl index 550efb6..b45c6b7 100644 --- a/src/get.jl +++ b/src/get.jl @@ -18,7 +18,7 @@ For the inverse operation, see the [`CoordRefSystems.code`](@ref) function. """ function get(code::Type{<:CRSCode}) throw(ArgumentError(""" - The provided code $code is not mapped to a CRS type yet. + The provided code $code is not mapped to a CRS type yet. Please check https://github.com/JuliaEarth/CoordRefSystems.jl/blob/main/src/get.jl If you know the CRS type of a given code, please submit a pull request. See https://discourse.julialang.org/t/esri-code-for-british-national-grid-not-known-by-geoio/117641 @@ -63,6 +63,7 @@ end @crscode EPSG{4208} LatLon{Aratu} @crscode EPSG{4269} LatLon{NAD83} @crscode EPSG{4326} LatLon{WGS84Latest} +@crscode EPSG{5070} shift(Albers{23.0°,29.5°,45.5°,NAD83}, lonₒ=-96.0°) @crscode EPSG{4618} LatLon{SAD69} @crscode EPSG{4674} LatLon{SIRGAS2000} @crscode EPSG{4988} Cartesian3D{shift(ITRF{2000}, 2000.4)} diff --git a/src/strings.jl b/src/strings.jl index 5e0b061..0b46b39 100644 --- a/src/strings.jl +++ b/src/strings.jl @@ -14,6 +14,7 @@ const esriid2code = Dict( "GCS_South_American_1969" => EPSG{4618}, "GCS_WGS_1984" => EPSG{4326}, "IRENET95_Irish_Transverse_Mercator" => EPSG{2157}, + "NAD_1983_Contiguous_USA_Albers" => EPSG{5070}, "North_Pole_Orthographic" => ESRI{102035}, "South_Pole_Orthographic" => ESRI{102037}, "TM75_Irish_Grid" => EPSG{29903}, @@ -50,7 +51,7 @@ function string2code(crsstr) type == "EPSG" ? EPSG{code} : ESRI{code} elseif endswith(keyword, "CS") # WKT1 # use the datum name to differentiate ESRI WKT1 from OGC WKT1 - # if the datum name starts with "D_", the format is ESRI WKT1 + # if the datum name starts with "D_", the format is ESRI WKT1 datumregex = r"DATUM\[\"(.*?)\"" datummatch = match(datumregex, content) |> checkmatch datumname = datummatch.captures[1] 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) diff --git a/test/crs/constructors.jl b/test/crs/constructors.jl index aad0479..a03bbe7 100644 --- a/test/crs/constructors.jl +++ b/test/crs/constructors.jl @@ -623,5 +623,32 @@ @test_throws ArgumentError TM(T(1) * m, T(1) * s) @test_throws ArgumentError TM(T(1) * s, T(1) * s) end + + @testset "Albers" begin + A = Albers{23.0°,29.5°,45.5°} + @test A(T(1), T(1)) == A(T(1) * m, T(1) * m) + @test A(T(1) * m, 1 * m) == A(T(1) * m, T(1) * m) + @test A(T(1) * km, T(1) * km) == A(T(1000) * m, T(1000) * m) + + c = A(T(1), T(1)) + @test sprint(show, c) == "Albers{WGS84Latest}(x: 1.0 m, y: 1.0 m)" + if T === Float32 + @test sprint(show, MIME("text/plain"), c) == """ + Albers{WGS84Latest} coordinates + ├─ x: 1.0f0 m + └─ y: 1.0f0 m""" + else + @test sprint(show, MIME("text/plain"), c) == """ + Albers{WGS84Latest} coordinates + ├─ x: 1.0 m + └─ y: 1.0 m""" + end + + # error: invalid units for coordinates + @test_throws ArgumentError A(T(1), T(1) * m) + @test_throws ArgumentError A(T(1) * s, T(1) * m) + @test_throws ArgumentError A(T(1) * m, T(1) * s) + @test_throws ArgumentError A(T(1) * s, T(1) * s) + end end end diff --git a/test/crs/conversions.jl b/test/crs/conversions.jl index 057cd77..d4add46 100644 --- a/test/crs/conversions.jl +++ b/test/crs/conversions.jl @@ -1365,6 +1365,40 @@ @inferred convert(LatLon, c2) end + @testset "LatLon <> Albers" begin + # tested aganist Proj.Transformation("EPSG:4269", "EPSG:5070") + AlbersUS = CoordRefSystems.shift(Albers{23.0°,29.5°,45.5°,NAD83}, lonₒ=-96°) + c1 = LatLon{NAD83}(T(45), T(90)) + c2 = convert(AlbersUS, c1) + @test allapprox(c2, AlbersUS(-T(7.231430540202629e6), T(1.1853758709623523e7))) + c3 = convert(LatLon{NAD83}, c2) + @test allapprox(c3, c1) + + c1 = LatLon{NAD83}(-T(45), T(90)) + c2 = convert(AlbersUS, c1) + @test allapprox(c2, AlbersUS(-T(1.5156419949174397e7), T(1.3963188032694416e7))) + c3 = convert(LatLon{NAD83}, c2) + @test allapprox(c3, c1) + + c1 = LatLon{NAD83}(T(45), -T(90)) + c2 = convert(AlbersUS, c1) + @test allapprox(c2, AlbersUS(T(472145.2567221438), T(2.4606302944375253e6))) + c3 = convert(LatLon{NAD83}, c2) + @test allapprox(c3, c1) + + c1 = LatLon{NAD83}(-T(45), -T(90)) + c2 = convert(AlbersUS, c1) + @test allapprox(c2, AlbersUS(T(989573.4665648951), -T(5.723953827495912e6))) + c3 = convert(LatLon{NAD83}, c2) + @test allapprox(c3, c1) + + # type stability + c1 = LatLon{NAD83}(T(45), T(90)) + c2 = AlbersUS(-T(7.231430540202629e6), T(1.1853758709623523e7)) + @inferred convert(AlbersUS, c1) + @inferred convert(LatLon{NAD83}, c2) + end + @testset "LatLon <> UTM" begin UTMNorth32 = utmnorth(32) UTMSouth59 = utmsouth(59) diff --git a/test/crs/crsapi.jl b/test/crs/crsapi.jl index 0d3d18e..c87c90d 100644 --- a/test/crs/crsapi.jl +++ b/test/crs/crsapi.jl @@ -2,6 +2,7 @@ ShiftedMercator = CoordRefSystems.shift(Mercator{WGS84Latest}, lonₒ=15.0°, xₒ=200.0m, yₒ=200.0m) ShiftedTM = CoordRefSystems.shift(TransverseMercator{0.9996,15.0°,WGS84Latest}, lonₒ=45.0°) UTMNorth32 = utmnorth(32) + AlbersUS = CoordRefSystems.shift(Albers{23.0°,29.5°,45.0°,NAD83}, lonₒ=-96.0°) @testset "datum" begin c = Cartesian(T(1), T(1)) @@ -163,6 +164,8 @@ @test CoordRefSystems.raw(c) == (T(1), T(2)) c = ShiftedMercator(T(1), T(2)) @test CoordRefSystems.raw(c) == (T(1), T(2)) + c = AlbersUS(T(1), T(2)) + @test CoordRefSystems.raw(c) == (T(1), T(2)) end @testset "units" begin @@ -206,6 +209,8 @@ @test CoordRefSystems.units(c) == (m, m) c = ShiftedMercator(T(1), T(2)) @test CoordRefSystems.units(c) == (m, m) + c = AlbersUS(T(1), T(2)) + @test CoordRefSystems.units(c) == (m, m) end @testset "constructor" begin @@ -247,6 +252,8 @@ @test CoordRefSystems.constructor(c) === UTMNorth32 c = ShiftedMercator(T(1), T(2)) @test CoordRefSystems.constructor(c) === ShiftedMercator + c = AlbersUS(T(1), T(2)) + @test CoordRefSystems.constructor(c) === AlbersUS end @testset "reconstruct" begin @@ -310,6 +317,9 @@ c = ShiftedMercator(T(1), T(2)) rv = CoordRefSystems.raw(c) @test CoordRefSystems.reconstruct(typeof(c), rv) == c + c = AlbersUS(T(1), T(2)) + rv = CoordRefSystems.raw(c) + @test CoordRefSystems.reconstruct(typeof(c), rv) == c end @testset "lentype" begin @@ -351,6 +361,8 @@ @test CoordRefSystems.lentype(c) == typeof(T(1) * m) c = ShiftedMercator(T(1), T(2)) @test CoordRefSystems.lentype(c) == typeof(T(1) * m) + c = AlbersUS(T(1), T(2)) + @test CoordRefSystems.lentype(c) == typeof(T(1) * m) end @testset "mactype" begin @@ -392,6 +404,8 @@ @test CoordRefSystems.mactype(c) == T c = ShiftedMercator(T(1), T(2)) @test CoordRefSystems.mactype(c) == T + c = AlbersUS(T(1), T(2)) + @test CoordRefSystems.mactype(c) == T end @testset "equality" begin @@ -415,6 +429,7 @@ equaltest(UTMNorth32) equaltest(ShiftedTM) equaltest(ShiftedMercator) + equaltest(AlbersUS) end @testset "isapprox" begin diff --git a/test/crs/domains.jl b/test/crs/domains.jl index 25452b3..f61fe86 100644 --- a/test/crs/domains.jl +++ b/test/crs/domains.jl @@ -211,6 +211,23 @@ end end + @testset "Albers" begin + atol = T === Float32 ? 1.0f-1° : 1e-7° + AlbersUS = CoordRefSystems.shift(Albers{23.0°,29.5°,45.5°,NAD83}, lonₒ=-96.0°) + for lat in T.(-90:90), lon in T.(-180:180) + c1 = LatLon{NAD83}(lat, lon) + if indomain(AlbersUS, c1) + c2 = convert(AlbersUS, c1) + @test isfinite(c2.x) + @test isfinite(c2.y) + c3 = convert(LatLon{NAD83}, c2) + @test allapprox(c3, c1; atol) + else + @test_throws ArgumentError convert(AlbersUS, c1) + end + end + end + @testset "UTMNorth forward" begin UTMNorth32 = utmnorth(32) for lat in T.(-90:90), lon in T.(-180:180) diff --git a/test/crs/mactype.jl b/test/crs/mactype.jl index 7e2b8bf..50ae322 100644 --- a/test/crs/mactype.jl +++ b/test/crs/mactype.jl @@ -141,6 +141,15 @@ c2 = convert(C, c1) @test c2 isa C + AlbersUS = CoordRefSystems.shift(Albers{23.0°,29.5°,45.5°,NAD83}, lonₒ=-96.0°) + C = AlbersUS{Met{T}} + c1 = AlbersUS(1.0, 1.0) + c2 = convert(C, c1) + @test c2 isa C + c1 = AlbersUS(1.0f0, 1.0f0) + c2 = convert(C, c1) + @test c2 isa C + UTMNorth32 = utmnorth(32) C = UTMNorth32{Met{T}} c1 = UTMNorth32(1.0, 1.0) diff --git a/test/crs/rand.jl b/test/crs/rand.jl index 295143e..aca5698 100644 --- a/test/crs/rand.jl +++ b/test/crs/rand.jl @@ -69,5 +69,8 @@ randtest(OrthoSouth{WGS84Latest}) randtest(OrthoSouth) + + randtest(Albers{23.0°,29.5°,45.5°,WGS84Latest}) + randtest(Albers{23.0°,29.5°,45.5°}) end end diff --git a/test/get.jl b/test/get.jl index a5b2ba6..153342a 100644 --- a/test/get.jl +++ b/test/get.jl @@ -9,6 +9,7 @@ gettest(EPSG{4208}, LatLon{Aratu}) gettest(EPSG{4269}, LatLon{NAD83}) gettest(EPSG{4326}, LatLon{WGS84Latest}) + gettest(EPSG{5070}, CoordRefSystems.shift(Albers{23.0°,29.5°,45.5°,NAD83}, lonₒ=-96.0°)) gettest(EPSG{4618}, LatLon{SAD69}) gettest(EPSG{4674}, LatLon{SIRGAS2000}) gettest(EPSG{4988}, Cartesian{CoordRefSystems.shift(ITRF{2000}, 2000.4),3})