Adding support for EPSG:5070 and Albers projections generally

src/CoordRefSystems.jl

@@ -89,6 +89,7 @@ export
   OrthoNorth,
   OrthoSouth,
   TransverseMercator,
+  Albers,
   datum,
   indomain,
 

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/eqareaalbers.jl")
 
 # ----------
 # FALLBACKS

src/crs/projected/eqareaalbers.jl

@@ -0,0 +1,155 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    Albers(lat₀, lat₁, lat₂, Datum, Shift)
+    Albers CRS with latitude origin lat₀ standard parallels `lat₁` and `lat₂`, `lon₀`, `Datum` and `Shift`.
+
+## Examples
+
+```julia
+Albers(1, 1) # add default units
+Albers(1m, 1m) # integers are converted converted to floats
+Albers(1.0km, 1.0km) # length quantities are converted to meters
+Albers(1.0m, 1.0m)
+Albers{NAD83}(1.0m, 1.0m)
+```
+
+See [EPSG:5070](https://epsg.io/5070).
+"""
+struct Albers{latₒ,lat₁,lat₂,Datum,Shift,M<:Met} <: Projected{Datum,Shift}
+  x::M
+  y::M
+end
+
+Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift}(x::M, y::M) where {lat₀,lat₁,lat₂,lon₀,Datum,Shift,M<:Met} =
+  Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift,float(M)}(x, y)
+Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift}(x::Met, y::Met) where {lat₀,lat₁,lat₂,lon₀,Datum,Shift} =
+  Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift}(promote(x, y)...)
+Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift}(x::Len, y::Len) where {lat₀,lat₁,lat₂,lon₀,Datum,Shift} =
+  Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift}(uconvert(m, x), uconvert(m, y))
+Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift}(x::Number, y::Number) where {lat₀,lat₁,lat₂,lon₀,Datum,Shift} =
+  Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift}(addunit(x, m), addunit(y, m))
+
+Albers{lat₀,lat₁,lat₂,lon₀,Datum}(args...) where {lat₀,lat₁,lat₂,lon₀,Datum} =
+  Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift()}(args...)
+
+Albers(args...) = Albers{NAD83}(args...)
+
+Base.convert(
+  ::Type{Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift,M}},
+  coords::Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift}
+) where {lat₀,lat₁,lat₂,lon₀,Datum,Shift,M} = Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift,M}(coords.x, coords.y)
+
+constructor(::Type{<:Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift}}) where {lat₀,lat₁,lat₂,lon₀,Datum,Shift} =
+  Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift}
+
+lentype(::Type{<:Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift,M}}) where {lat₀,lat₁,lat₂,lon₀,Datum,Shift,M} = M
+
+==(
+  coords₁::Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift},
+  coords₂::Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift}
+) where {lat₀,lat₁,lat₂,lon₀,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
+
+inbounds(::Type{<:Albers{lat₀,lat₁,lat₂,lon₀,Datum}}, λ, ϕ) where {lat₀,lat₁,lat₂,lon₀,Datum} =
+  -π ≤ λ ≤ π && deg2rad(90) ≤ ϕ ≤ deg2rad(90)
+
+function formulas(::Type{<:Albers{Datum}}, ::Type{T}) where {lat₀,lat₁,lat₂,lon₀,Datum,T}
+  # Constants
+  🌎 = ellipsoid(Datum)
+  e = T(eccentricity(🌎))
+  a = numconvert(T, majoraxis(🌎))
+
+  # Latitude origin
+  ϕ₀ = T(ustrip(deg2rad(lat₀)))
+
+  # Standard parallels
+  ϕ₁ = T(ustrip(deg2rad(lat₁)))
+  ϕ₂ = T(ustrip(deg2rad(lat₂)))
+
+  # Longitude origin
+  λ₀ = T(ustrip(deg2rad(lon₀)))
+
+  m₁ = hm(ϕ₁, e)
+  m₂ = hm(ϕ₂, e)
+  α₁ = hα(ϕ₁, e)
+  α₂ = hα(ϕ₂, e)
+  α₀ = hα(ϕ₀, e)
+  n = (m₁^2 - m₂^2) / (α₂ - α₁)
+  C = m₁^2 + n * α₁
+
+  Θ = n * λ
+  ρ = a * sqrt(C - n * hα(ϕ, e)) / n
+  if ρ < 0
+    throw(ArgumentError("coordinates outside of the projection domain"))
+  end
+  ρ₀ = (a * (C - n * α₀))^0.5 / n
+  function fx(λ, ϕ)
+    ρ₀ - ρ * cos(Θ)
+  end
+
+  function fy(λ, ϕ)
+    ρ * sin(Θ)
+  end
+
+  fx, fy
+end
+
+# backward projection formulas
+
+function backward(::Type{<:Albers{Datum}}, x, y) where {lat₀,lat₁,lat₂,lon₀,Datum,T}
+  🌎 = ellipsoid(Datum)
+  e = oftype(x, eccentricity(🌎))
+  e² = oftype(x, eccentricity²(🌎))
+  ϕ₀, ϕ₁, ϕ₂, λ₀ = T.(ustrip.(deg2rad.(lat₀, lat₁, lat₂, lon₀)))
+  α₀ = hα(ϕ₀, e)
+  ρ₀ = (a * (C - n * α₀))^0.5 / n
+  α₁ = hα(ϕ₁, e)
+  α₂ = hα(ϕ₂, e)
+
+  m₁ = hm(ϕ₁, e)
+  m₂ = hm(ϕ₂, e)
+  n = (m₁^2 - m₂^2) / (α₂ - α₁)
+  C = m₁^2 + n * α₁
+
+  θ = atan2(x, ρ₀ - y)
+  ρ = sqrt(x^2 + (ρ₀ - y)^2)
+  α′ = (C - (ρ^2 * n^2) / a^2) / n
+  β′ = asin(α′ / (1 - (1 - e) / (2 * e) * log((1 - e) / (1 + e))))
+
+  λ = θ / n
+  ϕ =
+    β′ +
+    (e^2 / 3 + 31 * e^4 / 180 + 517 * e^6 / 5040) * sin(2 * β′) +
+    (23 * e^4 / 360 + 251 * e^6 / 3780) * sin(4 * β′) +
+    (761 * e^6 / 45360) * sin(6 * β′)
+  return λ, ϕ
+end
+# ----------
+# Helper functions
+# ----------
+function hm(ϕ, e)
+  cos(ϕ) / sqrt(1 - e^2 * sin(ϕ)^2)
+end
+function hα(ϕ, e)
+  (1 - e^2) * (sin(ϕ) / (1 - e^2 * sin(ϕ)^2) - (1 / (2 * e)) * log((1 - e * sin(ϕ)) / (1 + e * sin(ϕ))))
+end
+
+# ----------
+# FALLBACKS
+# ----------
+
+Base.convert(::Type{Albers{lat₀,lat₁,lat₂,lon₀}}, coords::CRS{Datum}) where {lat₀,lat₁,lat₂,lon₀,Datum} =
+  convert(Albers{lat₀,lat₁,lat₂,lon₀,Datum}, coords)

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)}

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}, 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})