From dd48c375052e9da6655b9cbd76a6181538113d45 Mon Sep 17 00:00:00 2001 From: souma4 Date: Thu, 24 Oct 2024 12:52:59 -0600 Subject: [PATCH 01/38] added export Albers --- src/CoordRefSystems.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/CoordRefSystems.jl b/src/CoordRefSystems.jl index 93eb46c1..681b3da0 100644 --- a/src/CoordRefSystems.jl +++ b/src/CoordRefSystems.jl @@ -89,6 +89,7 @@ export OrthoNorth, OrthoSouth, TransverseMercator, + Albers, datum, indomain, From 70a727a6d9c0c3a6445bfd60f4bc01cf2900b04c Mon Sep 17 00:00:00 2001 From: souma4 Date: Thu, 24 Oct 2024 12:53:45 -0600 Subject: [PATCH 02/38] includes alberts ea --- src/crs/projected.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/crs/projected.jl b/src/crs/projected.jl index 88c59c0c..21f3b018 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/eqareaalbers.jl") # ---------- # FALLBACKS From 75d8679b9de49704d12a16cfcc3bd475105a0974 Mon Sep 17 00:00:00 2001 From: souma4 Date: Thu, 24 Oct 2024 12:54:10 -0600 Subject: [PATCH 03/38] adding albers --- src/crs/projected/eqareaalbers.jl | 153 ++++++++++++++++++++++++++++++ 1 file changed, 153 insertions(+) create mode 100644 src/crs/projected/eqareaalbers.jl diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl new file mode 100644 index 00000000..880ebdbe --- /dev/null +++ b/src/crs/projected/eqareaalbers.jl @@ -0,0 +1,153 @@ +# ------------------------------------------------------------------ +# 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₂,lon₀,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(DomainError("Coordinate transformation outside 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) From 50b92c94f51780549598fc0e150e982aadd2d7a9 Mon Sep 17 00:00:00 2001 From: souma4 Date: Thu, 24 Oct 2024 12:54:47 -0600 Subject: [PATCH 04/38] updated get and tests --- src/get.jl | 3 ++- test/get.jl | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/get.jl b/src/get.jl index 550efb63..e11c678a 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,-96.0,NAD83}) @crscode EPSG{4618} LatLon{SAD69} @crscode EPSG{4674} LatLon{SIRGAS2000} @crscode EPSG{4988} Cartesian3D{shift(ITRF{2000}, 2000.4)} diff --git a/test/get.jl b/test/get.jl index a5b2ba67..20f76909 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}, shift(Albers{23.0,29.5,45.5,-96.0,NAD83})) gettest(EPSG{4618}, LatLon{SAD69}) gettest(EPSG{4674}, LatLon{SIRGAS2000}) gettest(EPSG{4988}, Cartesian{CoordRefSystems.shift(ITRF{2000}, 2000.4),3}) From 23aa433fed1b5baf35e40e38a648911d419cf4e3 Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Fri, 25 Oct 2024 10:20:02 -0600 Subject: [PATCH 05/38] Update src/crs/projected/eqareaalbers.jl Co-authored-by: Elias Carvalho <73039601+eliascarv@users.noreply.github.com> --- src/crs/projected/eqareaalbers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index 880ebdbe..e97bcaef 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -18,7 +18,7 @@ Albers{NAD83}(1.0m, 1.0m) See [EPSG:5070](https://epsg.io/5070). """ -struct Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift,M<:Met} <: Projected{Datum,Shift} +struct Albers{latₒ,lat₁,lat₂,Datum,Shift,M<:Met} <: Projected{Datum,Shift} x::M y::M end From df36644f6612b29b76dec7c4bf5826d161ee5b0b Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Fri, 25 Oct 2024 10:20:17 -0600 Subject: [PATCH 06/38] Update src/crs/projected/eqareaalbers.jl Co-authored-by: Elias Carvalho <73039601+eliascarv@users.noreply.github.com> --- src/crs/projected/eqareaalbers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index e97bcaef..5a4a0daa 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -88,7 +88,7 @@ function formulas(::Type{<:Albers{Datum}}, ::Type{T}) where {lat₀,lat₁,lat n = (m₁^2 - m₂^2) / (α₂ - α₁) C = m₁^2 + n * α₁ - Θ = n * (λ - λ₀) + Θ = n * λ ρ = a * sqrt(C - n * hα(ϕ, e)) / n if ρ < 0 throw(DomainError("Coordinate transformation outside projection domain")) From c8ba857f17c67a1bd0a48ef875330b4ffd87af80 Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Fri, 25 Oct 2024 10:20:23 -0600 Subject: [PATCH 07/38] Update src/crs/projected/eqareaalbers.jl Co-authored-by: Elias Carvalho <73039601+eliascarv@users.noreply.github.com> --- src/crs/projected/eqareaalbers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index 5a4a0daa..9f89bd9a 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -127,7 +127,7 @@ function backward(::Type{<:Albers{Datum}}, x, y) where {lat₀,lat₁,lat₂,lon α′ = (C - (ρ^2 * n^2) / a^2) / n β′ = asin(α′ / (1 - (1 - e) / (2 * e) * log((1 - e) / (1 + e)))) - λ = λ₀ + θ / n + λ = θ / n ϕ = β′ + (e^2 / 3 + 31 * e^4 / 180 + 517 * e^6 / 5040) * sin(2 * β′) + From 533d8c0912904cb11f1b0ae2ff3ed261be051f6b Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Fri, 25 Oct 2024 10:20:35 -0600 Subject: [PATCH 08/38] Update src/crs/projected/eqareaalbers.jl Co-authored-by: Elias Carvalho <73039601+eliascarv@users.noreply.github.com> --- src/crs/projected/eqareaalbers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index 9f89bd9a..bf4dee21 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -91,7 +91,7 @@ function formulas(::Type{<:Albers{Datum}}, ::Type{T}) where {lat₀,lat₁,lat Θ = n * λ ρ = a * sqrt(C - n * hα(ϕ, e)) / n if ρ < 0 - throw(DomainError("Coordinate transformation outside projection domain")) + throw(ArgumentError("coordinates outside of the projection domain")) end ρ₀ = (a * (C - n * α₀))^0.5 / n function fx(λ, ϕ) From 969996a4940daf17a7a5df0bb7c5d9b3f783e79d Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Fri, 25 Oct 2024 10:20:57 -0600 Subject: [PATCH 09/38] Update src/crs/projected/eqareaalbers.jl Co-authored-by: Elias Carvalho <73039601+eliascarv@users.noreply.github.com> --- src/crs/projected/eqareaalbers.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index bf4dee21..0e5c926d 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -51,6 +51,8 @@ lentype(::Type{<:Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift,M}}) where {lat 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 # ------------ From dad2dffc2ca2fcee5c278b9119207a9dcb79b37f Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Fri, 25 Oct 2024 10:23:48 -0600 Subject: [PATCH 10/38] Update get.jl to reflect lon\_o changes --- src/get.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/get.jl b/src/get.jl index e11c678a..8ee1b414 100644 --- a/src/get.jl +++ b/src/get.jl @@ -63,7 +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,-96.0,NAD83}) +@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)} From b1aff8baf64feb6c0220129e786f1bb5d4c0df85 Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Fri, 25 Oct 2024 10:24:34 -0600 Subject: [PATCH 11/38] Update get.jl test to reflect lon\_o changes --- test/get.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/get.jl b/test/get.jl index 20f76909..b5d8101c 100644 --- a/test/get.jl +++ b/test/get.jl @@ -9,7 +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,-96.0,NAD83})) + 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}) From 7ced673bcddbe69ea47dde24f4f1429a5bf89db1 Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Fri, 25 Oct 2024 10:51:03 -0600 Subject: [PATCH 12/38] Change lat\_0 to lat\_o. remove lon\_0 --- src/crs/projected/eqareaalbers.jl | 58 +++++++++++++++---------------- 1 file changed, 28 insertions(+), 30 deletions(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index 0e5c926d..9bca57dd 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -3,8 +3,8 @@ # ------------------------------------------------------------------ """ - Albers(lat₀, lat₁, lat₂, Datum, Shift) - Albers CRS with latitude origin lat₀ standard parallels `lat₁` and `lat₂`, `lon₀`, `Datum` and `Shift`. + Albers(latₒ, lat₁, lat₂, Datum, Shift) + Albers CRS with latitude origin latₒ standard parallels `lat₁` and `lat₂`, `Datum` and `Shift`. ## Examples @@ -23,34 +23,34 @@ struct Albers{latₒ,lat₁,lat₂,Datum,Shift,M<:Met} <: Projected{Datum,Shift} 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₂,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₂,lon₀,Datum}(args...) where {lat₀,lat₁,lat₂,lon₀,Datum} = - Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift()}(args...) +Albers{latₒ,lat₁,lat₂,Datum}(args...) where {latₒ,lat₁,lat₂,Datum} = + Albers{latₒ,lat₁,lat₂,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) + ::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₂,lon₀,Datum,Shift}}) where {lat₀,lat₁,lat₂,lon₀,Datum,Shift} = - Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift} +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₂,lon₀,Datum,Shift,M}}) where {lat₀,lat₁,lat₂,lon₀,Datum,Shift,M} = M +lentype(::Type{<:Albers{latₒ,lat₁,lat₂,Datum,Shift,M}}) where {latₒ,lat₁,lat₂,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 + 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 @@ -63,24 +63,22 @@ lentype(::Type{<:Albers{lat₀,lat₁,lat₂,lon₀,Datum,Shift,M}}) where {lat # 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} = +inbounds(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, λ, ϕ) where {latₒ,lat₁,lat₂,Datum} = -π ≤ λ ≤ π && deg2rad(90) ≤ ϕ ≤ deg2rad(90) -function formulas(::Type{<:Albers{Datum}}, ::Type{T}) where {lat₀,lat₁,lat₂,lon₀,Datum,T} +function formulas(::Type{<:Albers{Datum}}, ::Type{T}) where {latₒ,lat₁,lat₂,Datum,T} # Constants 🌎 = ellipsoid(Datum) e = T(eccentricity(🌎)) a = numconvert(T, majoraxis(🌎)) # Latitude origin - ϕ₀ = T(ustrip(deg2rad(lat₀))) + ϕ₀ = 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) @@ -109,11 +107,11 @@ end # backward projection formulas -function backward(::Type{<:Albers{Datum}}, x, y) where {lat₀,lat₁,lat₂,lon₀,Datum,T} +function backward(::Type{<:Albers{Datum}}, x, y) where {latₒ,lat₁,lat₂,Datum,T} 🌎 = ellipsoid(Datum) e = oftype(x, eccentricity(🌎)) e² = oftype(x, eccentricity²(🌎)) - ϕ₀, ϕ₁, ϕ₂, λ₀ = T.(ustrip.(deg2rad.(lat₀, lat₁, lat₂, lon₀))) + ϕ₀, ϕ₁, ϕ₂ = T.(ustrip.(deg2rad.(latₒ, lat₁, lat₂))) α₀ = hα(ϕ₀, e) ρ₀ = (a * (C - n * α₀))^0.5 / n α₁ = hα(ϕ₁, e) @@ -151,5 +149,5 @@ 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) +Base.convert(::Type{Albers{latₒ,lat₁,lat₂}}, coords::CRS{Datum}) where {latₒ,lat₁,lat₂,Datum} = + convert(Albers{latₒ,lat₁,lat₂,Datum}, coords) From 8aca199f24cd63d84841cadeff6db58ba63bc9d4 Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Fri, 25 Oct 2024 11:37:26 -0600 Subject: [PATCH 13/38] remove redundant type params Co-authored-by: Elias Carvalho <73039601+eliascarv@users.noreply.github.com> --- src/crs/projected/eqareaalbers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index 9bca57dd..8dbebc85 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -63,7 +63,7 @@ lentype(::Type{<:Albers{latₒ,lat₁,lat₂,Datum,Shift,M}}) where {latₒ,lat # 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₂,Datum}}, λ, ϕ) where {latₒ,lat₁,lat₂,Datum} = +inbounds(::Type{<:Albers}, λ, ϕ) = -π ≤ λ ≤ π && deg2rad(90) ≤ ϕ ≤ deg2rad(90) function formulas(::Type{<:Albers{Datum}}, ::Type{T}) where {latₒ,lat₁,lat₂,Datum,T} From 59bdc653bc7b6a56c409adb2e7bba919622fb7ab Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Fri, 25 Oct 2024 11:37:58 -0600 Subject: [PATCH 14/38] add missing type parameters Co-authored-by: Elias Carvalho <73039601+eliascarv@users.noreply.github.com> --- src/crs/projected/eqareaalbers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index 8dbebc85..3ce3c117 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -66,7 +66,7 @@ lentype(::Type{<:Albers{latₒ,lat₁,lat₂,Datum,Shift,M}}) where {latₒ,lat inbounds(::Type{<:Albers}, λ, ϕ) = -π ≤ λ ≤ π && deg2rad(90) ≤ ϕ ≤ deg2rad(90) -function formulas(::Type{<:Albers{Datum}}, ::Type{T}) where {latₒ,lat₁,lat₂,Datum,T} +function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where {latₒ,lat₁,lat₂,Datum,T} # Constants 🌎 = ellipsoid(Datum) e = T(eccentricity(🌎)) From 046536221cd865145af2bfe135ed45f2e83ead2e Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Fri, 25 Oct 2024 11:42:40 -0600 Subject: [PATCH 15/38] Apply suggestions from code review Various code style changes, changing ^0.5 to sqrt, inlined small helper function, split vectorized array operations to convert degrees to radians into individual lines, fixed missing and redundant type parameters. Learning! Co-authored-by: Elias Carvalho <73039601+eliascarv@users.noreply.github.com> --- src/crs/projected/eqareaalbers.jl | 31 ++++++++++++++++--------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index 3ce3c117..f19f04ef 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -79,7 +79,6 @@ function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where ϕ₁ = T(ustrip(deg2rad(lat₁))) ϕ₂ = T(ustrip(deg2rad(lat₂))) - m₁ = hm(ϕ₁, e) m₂ = hm(ϕ₂, e) α₁ = hα(ϕ₁, e) @@ -93,7 +92,7 @@ function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where if ρ < 0 throw(ArgumentError("coordinates outside of the projection domain")) end - ρ₀ = (a * (C - n * α₀))^0.5 / n + ρ₀ = sqrt(a * (C - n * α₀)) / n function fx(λ, ϕ) ρ₀ - ρ * cos(Θ) end @@ -107,13 +106,15 @@ end # backward projection formulas -function backward(::Type{<:Albers{Datum}}, x, y) where {latₒ,lat₁,lat₂,Datum,T} +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²(🌎)) - ϕ₀, ϕ₁, ϕ₂ = T.(ustrip.(deg2rad.(latₒ, lat₁, lat₂))) + ϕ₀ = oftype(x, ustrip(deg2rad(latₒ))) + ϕ₁ = oftype(x, ustrip(deg2rad(lat₁))) + ϕ₂ = oftype(x, ustrip(deg2rad(lat₂))) α₀ = hα(ϕ₀, e) - ρ₀ = (a * (C - n * α₀))^0.5 / n + ρ₀ = sqrt(a * (C - n * α₀)) / n α₁ = hα(ϕ₁, e) α₂ = hα(ϕ₂, e) @@ -133,18 +134,18 @@ function backward(::Type{<:Albers{Datum}}, x, y) where {latₒ,lat₁,lat₂,Dat (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 +# ----------------- +# HELPER FUNCTIONS +# ----------------- + +hm(ϕ, e) = cos(ϕ) / sqrt(1 - e^2 * sin(ϕ)^2) + +hα(ϕ, e) = (1 - e^2) * (sin(ϕ) / (1 - e^2 * sin(ϕ)^2) - (1 / (2 * e)) * log((1 - e * sin(ϕ)) / (1 + e * sin(ϕ)))) + # ---------- # FALLBACKS # ---------- From 79fb6715516d25b654611f831ac591cc9a32d0c2 Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Fri, 25 Oct 2024 11:43:53 -0600 Subject: [PATCH 16/38] newline between docstring definition and description Co-authored-by: Elias Carvalho <73039601+eliascarv@users.noreply.github.com> --- src/crs/projected/eqareaalbers.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index f19f04ef..440874f6 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -3,8 +3,9 @@ # ------------------------------------------------------------------ """ - Albers(latₒ, lat₁, lat₂, Datum, Shift) - Albers CRS with latitude origin latₒ standard parallels `lat₁` and `lat₂`, `Datum` and `Shift`. + Alber{latₒ, lat₁, lat₂, Datum, Shift} + +Albers CRS with latitude origin latₒ standard parallels `lat₁` and `lat₂`, `Datum` and `Shift`. ## Examples From ec50d6b15aedc7eacba981952575ea97470df00c Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Fri, 25 Oct 2024 12:08:53 -0600 Subject: [PATCH 17/38] Forgot to update all `\_0` to `\_o`, so fixed --- src/crs/projected/eqareaalbers.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index 440874f6..a133bb2e 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -74,7 +74,7 @@ function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where a = numconvert(T, majoraxis(🌎)) # Latitude origin - ϕ₀ = T(ustrip(deg2rad(latₒ))) + ϕₒ = T(ustrip(deg2rad(latₒ))) # Standard parallels ϕ₁ = T(ustrip(deg2rad(lat₁))) @@ -84,7 +84,7 @@ function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where m₂ = hm(ϕ₂, e) α₁ = hα(ϕ₁, e) α₂ = hα(ϕ₂, e) - α₀ = hα(ϕ₀, e) + αₒ = hα(ϕₒ, e) n = (m₁^2 - m₂^2) / (α₂ - α₁) C = m₁^2 + n * α₁ @@ -93,9 +93,9 @@ function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where if ρ < 0 throw(ArgumentError("coordinates outside of the projection domain")) end - ρ₀ = sqrt(a * (C - n * α₀)) / n + ρₒ = sqrt(a * (C - n * αₒ)) / n function fx(λ, ϕ) - ρ₀ - ρ * cos(Θ) + ρₒ - ρ * cos(Θ) end function fy(λ, ϕ) @@ -111,11 +111,11 @@ function backward(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, x, y) where {lat 🌎 = 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₁))) ϕ₂ = oftype(x, ustrip(deg2rad(lat₂))) - α₀ = hα(ϕ₀, e) - ρ₀ = sqrt(a * (C - n * α₀)) / n + αₒ = hα(ϕₒ, e) + ρₒ = sqrt(a * (C - n * αₒ)) / n α₁ = hα(ϕ₁, e) α₂ = hα(ϕ₂, e) @@ -124,8 +124,8 @@ function backward(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, x, y) where {lat n = (m₁^2 - m₂^2) / (α₂ - α₁) C = m₁^2 + n * α₁ - θ = atan2(x, ρ₀ - y) - ρ = sqrt(x^2 + (ρ₀ - y)^2) + θ = 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)))) From d741fa0e8bf4deac3e91c44b7e2852e0f5bd0647 Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Fri, 25 Oct 2024 15:04:22 -0600 Subject: [PATCH 18/38] Applied suggestions from code review. Restructured eqareaalbers.jl removed divide by a since this is performed in another step, reordered constants, shifted constants so that all operations on \lambda and \phi occur within locally defined functions Co-authored-by: Elias Carvalho <73039601+eliascarv@users.noreply.github.com> --- src/crs/projected/eqareaalbers.jl | 37 ++++++++++--------------------- 1 file changed, 12 insertions(+), 25 deletions(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index a133bb2e..a0fd2d69 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -68,15 +68,9 @@ inbounds(::Type{<:Albers}, λ, ϕ) = -π ≤ λ ≤ π && deg2rad(90) ≤ ϕ ≤ deg2rad(90) function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where {latₒ,lat₁,lat₂,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₂))) @@ -84,23 +78,16 @@ function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where 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 - ρₒ = sqrt(a * (C - n * αₒ)) / n - function fx(λ, ϕ) - ρₒ - ρ * cos(Θ) - end + Θ(λ) = n * λ + ρ(ϕ) = sqrt(C - n * hα(ϕ, e)) / n + ρₒ = ρ(ϕₒ) - function fy(λ, ϕ) - ρ * sin(Θ) - end + fx(λ, ϕ) = ρ(ϕ) * sin(Θ(λ)) + + fy(λ, ϕ) = ρₒ - ρ(ϕ) * cos(Θ(λ)) fx, fy end @@ -114,19 +101,19 @@ function backward(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, x, y) where {lat ϕₒ = oftype(x, ustrip(deg2rad(latₒ))) ϕ₁ = oftype(x, ustrip(deg2rad(lat₁))) ϕ₂ = oftype(x, ustrip(deg2rad(lat₂))) - αₒ = hα(ϕₒ, e) - ρₒ = sqrt(a * (C - n * αₒ)) / n - α₁ = hα(ϕ₁, e) - α₂ = hα(ϕ₂, e) + ρₒ = sqrt(C - n * αₒ) / n m₁ = hm(ϕ₁, e) m₂ = hm(ϕ₂, e) + αₒ = hα(ϕₒ, e) + α₁ = hα(ϕ₁, e) + α₂ = hα(ϕ₂, 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 + ρ′ = sqrt(x^2 + (ρₒ - y)^2) + α′ = (C - (ρ′^2 * n^2)) / n β′ = asin(α′ / (1 - (1 - e) / (2 * e) * log((1 - e) / (1 + e)))) λ = θ / n From c3082799a1b1cb8e57932c36803a9ce41bf0b63a Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Fri, 25 Oct 2024 17:53:20 -0600 Subject: [PATCH 19/38] Applied suggestions from code review: updated to use e\^2 rather than computing each time, using auth2geod Co-authored-by: Elias Carvalho <73039601+eliascarv@users.noreply.github.com> --- src/crs/projected/eqareaalbers.jl | 31 ++++++++++++++----------------- 1 file changed, 14 insertions(+), 17 deletions(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index a0fd2d69..be45cb9a 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -70,19 +70,20 @@ inbounds(::Type{<:Albers}, λ, ϕ) = 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₂))) - m₁ = hm(ϕ₁, e) - m₂ = hm(ϕ₂, e) - α₁ = hα(ϕ₁, e) - α₂ = hα(ϕ₂, e) + m₁ = hm(ϕ₁, e, e²) + m₂ = hm(ϕ₂, e, e²) + α₁ = hα(ϕ₁, e, e²) + α₂ = hα(ϕ₂, e, e²) n = (m₁^2 - m₂^2) / (α₂ - α₁) C = m₁^2 + n * α₁ Θ(λ) = n * λ - ρ(ϕ) = sqrt(C - n * hα(ϕ, e)) / n + ρ(ϕ) = sqrt(C - n * hα(ϕ, e, e²)) / n ρₒ = ρ(ϕₒ) fx(λ, ϕ) = ρ(ϕ) * sin(Θ(λ)) @@ -103,11 +104,11 @@ function backward(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, x, y) where {lat ϕ₂ = oftype(x, ustrip(deg2rad(lat₂))) ρₒ = sqrt(C - n * αₒ) / n - m₁ = hm(ϕ₁, e) - m₂ = hm(ϕ₂, e) - αₒ = hα(ϕₒ, e) - α₁ = hα(ϕ₁, e) - α₂ = hα(ϕ₂, e) + m₁ = hm(ϕ₁, e, e²) + m₂ = hm(ϕ₂, e, e²) + αₒ = hα(ϕₒ, e, e²) + α₁ = hα(ϕ₁, e, e²) + α₂ = hα(ϕ₂, e, e²) n = (m₁^2 - m₂^2) / (α₂ - α₁) C = m₁^2 + n * α₁ @@ -117,11 +118,7 @@ function backward(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, x, y) where {lat β′ = 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 * β′) + ϕ = auth2geod(β′, e²) λ, ϕ end @@ -130,9 +127,9 @@ end # HELPER FUNCTIONS # ----------------- -hm(ϕ, e) = cos(ϕ) / sqrt(1 - e^2 * sin(ϕ)^2) +hm(ϕ, e, e²) = cos(ϕ) / sqrt(1 - e² * sin(ϕ)^2) -hα(ϕ, e) = (1 - e^2) * (sin(ϕ) / (1 - e^2 * sin(ϕ)^2) - (1 / (2 * e)) * log((1 - e * sin(ϕ)) / (1 + e * sin(ϕ)))) +hα(ϕ, e, e²) = (1 - e²) * (sin(ϕ) / (1 - e² * sin(ϕ)^2) - (1 / (2 * e)) * log((1 - e * sin(ϕ)) / (1 + e * sin(ϕ)))) # ---------- # FALLBACKS From 602f951e2c040235eedb73c636067e158b6170ac Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Fri, 25 Oct 2024 17:58:21 -0600 Subject: [PATCH 20/38] fixed docstring from Alber to Albers --- src/crs/projected/eqareaalbers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index be45cb9a..ea1b8129 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -3,7 +3,7 @@ # ------------------------------------------------------------------ """ - Alber{latₒ, lat₁, lat₂, Datum, Shift} + Albers{latₒ, lat₁, lat₂, Datum, Shift} Albers CRS with latitude origin latₒ standard parallels `lat₁` and `lat₂`, `Datum` and `Shift`. From 4a022d80265258dd2dc35e860855f85ef98deddc Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Fri, 25 Oct 2024 19:08:41 -0600 Subject: [PATCH 21/38] Updated test/get.jl. Added namespace to shift --- test/get.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/get.jl b/test/get.jl index b5d8101c..b720b65b 100644 --- a/test/get.jl +++ b/test/get.jl @@ -9,7 +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{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}) From b06aad4ce8c33b8995d5eca0eccbad456f5e9eae Mon Sep 17 00:00:00 2001 From: souma4 Date: Sat, 26 Oct 2024 16:41:44 -0600 Subject: [PATCH 22/38] src/crs/projected/eqareaalbers.jl updated so that inbounds all longitude, only standard parallels latitude and \lambda properly wraps within the function test/crs/* added tests for Albers in each I COULD NOT TEST AS OF THIS MOMENT. ALL APPROX WASNT ACCEPTING THE METHOD, I BELIEVE IT HAS TO DO WITH THE FACT THE CRS IS SHIFTED I DID DETECT > 1M INACCURACY FROM PROJ. WILL HAVE TO INVESTIGATE Additionally, selected coordinates for tests are not final. --- src/crs/projected/eqareaalbers.jl | 13 ++++++------ test/crs/constructors.jl | 27 +++++++++++++++++++++++++ test/crs/conversions.jl | 33 +++++++++++++++++++++++++++++++ test/crs/crsapi.jl | 15 ++++++++++++++ test/crs/domains.jl | 14 +++++++++++++ test/crs/mactype.jl | 9 +++++++++ test/crs/rand.jl | 3 +++ 7 files changed, 107 insertions(+), 7 deletions(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index ea1b8129..00bb5299 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -4,7 +4,7 @@ """ Albers{latₒ, lat₁, lat₂, Datum, Shift} - + Albers CRS with latitude origin latₒ standard parallels `lat₁` and `lat₂`, `Datum` and `Shift`. ## Examples @@ -33,8 +33,7 @@ Albers{latₒ,lat₁,lat₂,Datum,Shift}(x::Len, y::Len) where {latₒ,lat₁,la 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₂,Datum}(args...) where {latₒ,lat₁,lat₂,Datum} = Albers{latₒ,lat₁,lat₂,Datum,Shift()}(args...) Albers(args...) = Albers{NAD83}(args...) @@ -64,8 +63,7 @@ lentype(::Type{<:Albers{latₒ,lat₁,lat₂,Datum,Shift,M}}) where {latₒ,lat # 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}, λ, ϕ) = - -π ≤ λ ≤ π && deg2rad(90) ≤ ϕ ≤ deg2rad(90) +inbounds(::Type{<:Albers}, λ, ϕ) = -2π ≤ λ ≤ 2π && deg2rad(lat₁) ≤ ϕ ≤ deg2rad(lat₂) function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where {latₒ,lat₁,lat₂,Datum,T} 🌎 = ellipsoid(Datum) @@ -86,9 +84,9 @@ function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where ρ(ϕ) = sqrt(C - n * hα(ϕ, e, e²)) / n ρₒ = ρ(ϕₒ) - fx(λ, ϕ) = ρ(ϕ) * sin(Θ(λ)) + fx(λ, ϕ) = ρ(ϕ) * sin(Θ(hλ(λ))) - fy(λ, ϕ) = ρₒ - ρ(ϕ) * cos(Θ(λ)) + fy(λ, ϕ) = ρₒ - ρ(ϕ) * cos(Θ(hλ(λ))) fx, fy end @@ -131,6 +129,7 @@ hm(ϕ, e, e²) = cos(ϕ) / sqrt(1 - e² * sin(ϕ)^2) hα(ϕ, e, e²) = (1 - e²) * (sin(ϕ) / (1 - e² * sin(ϕ)^2) - (1 / (2 * e)) * log((1 - e * sin(ϕ)) / (1 + e * sin(ϕ)))) +hλ(λ) = λ > π ? λ - 2π : λ < -π ? λ + 2π : λ # ---------- # FALLBACKS # ---------- diff --git a/test/crs/constructors.jl b/test/crs/constructors.jl index aad0479c..19228016 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,NAD83} + @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{NAD83}(x: 1.0 m, y: 1.0 m)" + if T === Float32 + @test sprint(show, MIME("text/plain"), c) == """ + Albers{NAD83} coordinates + ├─ x: 1.0f0 m + └─ y: 1.0f0 m""" + else + @test sprint(show, MIME("text/plain"), c) == """ + Albers{NAD83} 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 057cd776..02f9fcd5 100644 --- a/test/crs/conversions.jl +++ b/test/crs/conversions.jl @@ -1365,6 +1365,39 @@ @inferred convert(LatLon, c2) end + @testset "LatLon <> Albers" begin + AlbersUS = CoordRefSystems.get(EPSG{5070}) + c1 = LatLon(T(45), T(90)) + c2 = convert(AlbersUS, c1) + @test allapprox(c2, AlbersUS(T(-7231430.540202629), T(11853758.709623523))) + c3 = convert(LatLon, c2) + @test allapprox(c3, c1) + + c1 = LatLon(-T(45), T(90)) + c2 = convert(AlbersUS, c1) + @test allapprox(c2, AlbersUS(T(10018754.171394622), -T(5591295.9185533915))) + c3 = convert(LatLon, c2) + @test allapprox(c3, c1) + + c1 = LatLon(T(45), -T(90)) + c2 = convert(AlbersUS, c1) + @test allapprox(c2, AlbersUS(-T(10018754.171394622), T(5591295.9185533915))) + c3 = convert(LatLon, c2) + @test allapprox(c3, c1) + + c1 = LatLon(-T(45), -T(90)) + c2 = convert(AlbersUS, c1) + @test allapprox(c2, AlbersUS(-T(10018754.171394622), -T(5591295.9185533915))) + c3 = convert(LatLon, c2) + @test allapprox(c3, c1) + + # type stability + c1 = LatLon(T(45), T(90)) + c2 = AlbersUS(T(-7231430.540202629), T(11853758.709623523)) + @inferred convert(AlbersUS, c1) + @inferred convert(LatLon, 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 0d3d18e7..1688bec1 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,WGS84Latest}, xₒ=-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 25452b3d..c5bf7ae1 100644 --- a/test/crs/domains.jl +++ b/test/crs/domains.jl @@ -211,6 +211,20 @@ end end + @testset "Albers Forward" begin + Albers = 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(lat, lon) + if indomain(Albers, c1) + c2 = convert(Albers, c1) + @test isfinite(c2.x) + @test isfinite(c2.y) + else + @test_throws ArgumentError convert(Albers, 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 7e2b8bfc..42092ad4 100644 --- a/test/crs/mactype.jl +++ b/test/crs/mactype.jl @@ -141,6 +141,15 @@ c2 = convert(C, c1) @test c2 isa C + Albers = CoordRefSystems.shift(Albers{23.0,29.5,45.5,NAD83}, lonₒ=-96.0°) + C = Albers{Met{T}} + c1 = Albers(1.0, 1.0) + c2 = convert(C, c1) + @test c2 isa C + c1 = Albers(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 295143e3..4d72140d 100644 --- a/test/crs/rand.jl +++ b/test/crs/rand.jl @@ -69,5 +69,8 @@ randtest(OrthoSouth{WGS84Latest}) randtest(OrthoSouth) + + randtest(Albers{NAD83}) + randtest(Albers) end end From 9607f281d5f0e365bca99690825e2874ffe537b6 Mon Sep 17 00:00:00 2001 From: souma4 Date: Sun, 27 Oct 2024 15:04:17 -0600 Subject: [PATCH 23/38] Made adjustments to Albers script and tests changed eqareaalbers to allow latitudes up to \pi and included a helper function to wrap latitudes if greater than abs(\pi/2) Fixed tests to properly construct an Albers object and test it Conversions needs more work. Some of the allapprox tests failed but have submeter inaccuracy. The real issue is in back transforming. Latitudes appear to be back transforming inappropriately, so that will have to be double checked --- src/crs/projected/eqareaalbers.jl | 19 ++++++++++++++----- test/crs/conversions.jl | 26 +++++++++++++------------- test/crs/domains.jl | 8 ++++---- test/crs/mactype.jl | 6 +++--- test/crs/rand.jl | 3 +-- 5 files changed, 35 insertions(+), 27 deletions(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index 00bb5299..7cc3adc8 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -63,7 +63,7 @@ lentype(::Type{<:Albers{latₒ,lat₁,lat₂,Datum,Shift,M}}) where {latₒ,lat # 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}, λ, ϕ) = -2π ≤ λ ≤ 2π && deg2rad(lat₁) ≤ ϕ ≤ deg2rad(lat₂) +inbounds(::Type{<:Albers}, λ, ϕ) = -2π ≤ λ ≤ 2π && -π ≤ ϕ ≤ π function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where {latₒ,lat₁,lat₂,Datum,T} 🌎 = ellipsoid(Datum) @@ -84,9 +84,14 @@ function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where ρ(ϕ) = sqrt(C - n * hα(ϕ, e, e²)) / n ρₒ = ρ(ϕₒ) - fx(λ, ϕ) = ρ(ϕ) * sin(Θ(hλ(λ))) + fx(λ, ϕ) = ρ(hϕ(ϕ)) * sin(Θ(hλ(λ))) - fy(λ, ϕ) = ρₒ - ρ(ϕ) * cos(Θ(hλ(λ))) + function fy(λ, ϕ) + ρ = ρₒ - ρ(hϕ(ϕ)) * cos(Θ(hλ(λ))) + if ρ < 0 + throw(ArgumentError("coordinates outside of the projection domain")) + end + end fx, fy end @@ -101,7 +106,6 @@ function backward(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, x, y) where {lat ϕ₁ = oftype(x, ustrip(deg2rad(lat₁))) ϕ₂ = oftype(x, ustrip(deg2rad(lat₂))) - ρₒ = sqrt(C - n * αₒ) / n m₁ = hm(ϕ₁, e, e²) m₂ = hm(ϕ₂, e, e²) αₒ = hα(ϕₒ, e, e²) @@ -110,7 +114,10 @@ function backward(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, x, y) where {lat n = (m₁^2 - m₂^2) / (α₂ - α₁) C = m₁^2 + n * α₁ - θ = atan2(x, ρₒ - y) + ρ(ϕ) = sqrt(C - n * hα(ϕ, e, e²)) / n + ρₒ = ρ(ϕₒ) + + θ = atan(x, ρₒ - y) ρ′ = sqrt(x^2 + (ρₒ - y)^2) α′ = (C - (ρ′^2 * n^2)) / n β′ = asin(α′ / (1 - (1 - e) / (2 * e) * log((1 - e) / (1 + e)))) @@ -130,6 +137,8 @@ hm(ϕ, e, e²) = cos(ϕ) / sqrt(1 - e² * sin(ϕ)^2) hα(ϕ, e, e²) = (1 - e²) * (sin(ϕ) / (1 - e² * sin(ϕ)^2) - (1 / (2 * e)) * log((1 - e * sin(ϕ)) / (1 + e * sin(ϕ)))) hλ(λ) = λ > π ? λ - 2π : λ < -π ? λ + 2π : λ + +hϕ(ϕ) = ϕ > π / 2 ? ϕ - π : ϕ < -π / 2 ? ϕ + π : ϕ # ---------- # FALLBACKS # ---------- diff --git a/test/crs/conversions.jl b/test/crs/conversions.jl index 02f9fcd5..b62ab6f1 100644 --- a/test/crs/conversions.jl +++ b/test/crs/conversions.jl @@ -1367,35 +1367,35 @@ @testset "LatLon <> Albers" begin AlbersUS = CoordRefSystems.get(EPSG{5070}) - c1 = LatLon(T(45), T(90)) + c1 = LatLon{NAD83}(T(45), T(90)) c2 = convert(AlbersUS, c1) @test allapprox(c2, AlbersUS(T(-7231430.540202629), T(11853758.709623523))) - c3 = convert(LatLon, c2) + c3 = convert(LatLon{NAD83}, c2) @test allapprox(c3, c1) - c1 = LatLon(-T(45), T(90)) + c1 = LatLon{NAD83}(-T(45), T(90)) c2 = convert(AlbersUS, c1) - @test allapprox(c2, AlbersUS(T(10018754.171394622), -T(5591295.9185533915))) - c3 = convert(LatLon, c2) + @test allapprox(c2, AlbersUS(-T(15156419.949174397), T(13963188.032694416))) + c3 = convert(LatLon{NAD83}, c2) @test allapprox(c3, c1) - c1 = LatLon(T(45), -T(90)) + c1 = LatLon{NAD83}(T(45), -T(90)) c2 = convert(AlbersUS, c1) - @test allapprox(c2, AlbersUS(-T(10018754.171394622), T(5591295.9185533915))) - c3 = convert(LatLon, c2) + @test allapprox(c2, AlbersUS(T(472145.12299166346), T(2460630.2617624514))) + c3 = convert(LatLon{NAD83}, c2) @test allapprox(c3, c1) - c1 = LatLon(-T(45), -T(90)) + c1 = LatLon{NAD83}(-T(45), -T(90)) c2 = convert(AlbersUS, c1) - @test allapprox(c2, AlbersUS(-T(10018754.171394622), -T(5591295.9185533915))) - c3 = convert(LatLon, c2) + @test allapprox(c2, AlbersUS(T(8640856.920871278), -T(4596167.748123343))) + c3 = convert(LatLon{NAD83}, c2) @test allapprox(c3, c1) # type stability - c1 = LatLon(T(45), T(90)) + c1 = LatLon{NAD83}(T(45), T(90)) c2 = AlbersUS(T(-7231430.540202629), T(11853758.709623523)) @inferred convert(AlbersUS, c1) - @inferred convert(LatLon, c2) + @inferred convert(LatLon{NAD83}, c2) end @testset "LatLon <> UTM" begin diff --git a/test/crs/domains.jl b/test/crs/domains.jl index c5bf7ae1..44da0b4c 100644 --- a/test/crs/domains.jl +++ b/test/crs/domains.jl @@ -212,15 +212,15 @@ end @testset "Albers Forward" begin - Albers = CoordRefSystems.shift(Albers{23.0,29.5,45.5,NAD83}, lonₒ=-96.0°) + 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(lat, lon) - if indomain(Albers, c1) - c2 = convert(Albers, c1) + if indomain(AlbersUS, c1) + c2 = convert(AlbersUS, c1) @test isfinite(c2.x) @test isfinite(c2.y) else - @test_throws ArgumentError convert(Albers, c1) + @test_throws ArgumentError convert(AlbersUS, c1) end end end diff --git a/test/crs/mactype.jl b/test/crs/mactype.jl index 42092ad4..1afed96e 100644 --- a/test/crs/mactype.jl +++ b/test/crs/mactype.jl @@ -141,12 +141,12 @@ c2 = convert(C, c1) @test c2 isa C - Albers = CoordRefSystems.shift(Albers{23.0,29.5,45.5,NAD83}, lonₒ=-96.0°) + AlbersUS = CoordRefSystems.shift(Albers{23.0,29.5,45.5,NAD83}, lonₒ=-96.0°) C = Albers{Met{T}} - c1 = Albers(1.0, 1.0) + c1 = AlbersUS(1.0, 1.0) c2 = convert(C, c1) @test c2 isa C - c1 = Albers(1.0f0, 1.0f0) + c1 = AlbersUS(1.0f0, 1.0f0) c2 = convert(C, c1) @test c2 isa C diff --git a/test/crs/rand.jl b/test/crs/rand.jl index 4d72140d..3728eeff 100644 --- a/test/crs/rand.jl +++ b/test/crs/rand.jl @@ -70,7 +70,6 @@ randtest(OrthoSouth{WGS84Latest}) randtest(OrthoSouth) - randtest(Albers{NAD83}) - randtest(Albers) + randtest(CoordRefSystems.get(EPSG{5070})) end end From 31d2648bd7e87c904badda9749868f4f3d6be941 Mon Sep 17 00:00:00 2001 From: souma4 Date: Sun, 27 Oct 2024 20:34:12 -0600 Subject: [PATCH 24/38] updated EPSG:5070 and related tests to have proper degree symbols and shifts. tests work on my machine EXCEPT for conversions. Something is wrong with the backward latitude function. The forward function for T(45), -T(90) and -T(45), -T(90) don't pass tests. Will have to investigate further. Reverted change to fy(), \rho probably isn't going to be <0 zero with my wrapping scheme. --- src/crs/projected/eqareaalbers.jl | 7 +------ src/get.jl | 2 +- test/crs/constructors.jl | 2 +- test/crs/conversions.jl | 2 +- test/crs/crsapi.jl | 2 +- test/crs/domains.jl | 2 +- test/crs/mactype.jl | 2 +- test/get.jl | 2 +- 8 files changed, 8 insertions(+), 13 deletions(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index 7cc3adc8..cf83084d 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -86,12 +86,7 @@ function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where fx(λ, ϕ) = ρ(hϕ(ϕ)) * sin(Θ(hλ(λ))) - function fy(λ, ϕ) - ρ = ρₒ - ρ(hϕ(ϕ)) * cos(Θ(hλ(λ))) - if ρ < 0 - throw(ArgumentError("coordinates outside of the projection domain")) - end - end + fy(λ, ϕ) = ρₒ - ρ(hϕ(ϕ)) * cos(Θ(hλ(λ))) fx, fy end diff --git a/src/get.jl b/src/get.jl index 8ee1b414..b45c6b79 100644 --- a/src/get.jl +++ b/src/get.jl @@ -63,7 +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{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/test/crs/constructors.jl b/test/crs/constructors.jl index 19228016..09bb3cc1 100644 --- a/test/crs/constructors.jl +++ b/test/crs/constructors.jl @@ -625,7 +625,7 @@ end @testset "Albers" begin - A = Albers{23.0,29.5,45.5,NAD83} + A = Albers{23.0°,29.5°,45.5°,NAD83} @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) diff --git a/test/crs/conversions.jl b/test/crs/conversions.jl index b62ab6f1..dacd54e4 100644 --- a/test/crs/conversions.jl +++ b/test/crs/conversions.jl @@ -1366,7 +1366,7 @@ end @testset "LatLon <> Albers" begin - AlbersUS = CoordRefSystems.get(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(-7231430.540202629), T(11853758.709623523))) diff --git a/test/crs/crsapi.jl b/test/crs/crsapi.jl index 1688bec1..c87c90dd 100644 --- a/test/crs/crsapi.jl +++ b/test/crs/crsapi.jl @@ -2,7 +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,WGS84Latest}, xₒ=-96.0) + AlbersUS = CoordRefSystems.shift(Albers{23.0°,29.5°,45.0°,NAD83}, lonₒ=-96.0°) @testset "datum" begin c = Cartesian(T(1), T(1)) diff --git a/test/crs/domains.jl b/test/crs/domains.jl index 44da0b4c..11c813a9 100644 --- a/test/crs/domains.jl +++ b/test/crs/domains.jl @@ -212,7 +212,7 @@ end @testset "Albers Forward" begin - AlbersUS = CoordRefSystems.shift(Albers{23.0,29.5,45.5,NAD83}, lonₒ=-96.0°) + 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(lat, lon) if indomain(AlbersUS, c1) diff --git a/test/crs/mactype.jl b/test/crs/mactype.jl index 1afed96e..0cd6f423 100644 --- a/test/crs/mactype.jl +++ b/test/crs/mactype.jl @@ -141,7 +141,7 @@ c2 = convert(C, c1) @test c2 isa C - AlbersUS = CoordRefSystems.shift(Albers{23.0,29.5,45.5,NAD83}, lonₒ=-96.0°) + AlbersUS = CoordRefSystems.shift(Albers{23.0°,29.5°,45.5°,NAD83}, lonₒ=-96.0°) C = Albers{Met{T}} c1 = AlbersUS(1.0, 1.0) c2 = convert(C, c1) diff --git a/test/get.jl b/test/get.jl index b720b65b..153342a8 100644 --- a/test/get.jl +++ b/test/get.jl @@ -9,7 +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{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}) From da6c794cfeed11e9b274ad0fb74098a3bbc71336 Mon Sep 17 00:00:00 2001 From: souma4 Date: Sun, 27 Oct 2024 22:02:30 -0600 Subject: [PATCH 25/38] Fixed \Beta\prime formula so back calculations are correct. All tests pass on my machine except for conversion tests 3 and 4. The forward calculation diverges from PROJ at sub meter --- src/crs/projected/eqareaalbers.jl | 5 ++--- test/crs/conversions.jl | 2 +- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index cf83084d..d29babed 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -103,7 +103,6 @@ function backward(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, x, y) where {lat m₁ = hm(ϕ₁, e, e²) m₂ = hm(ϕ₂, e, e²) - αₒ = hα(ϕₒ, e, e²) α₁ = hα(ϕ₁, e, e²) α₂ = hα(ϕ₂, e, e²) n = (m₁^2 - m₂^2) / (α₂ - α₁) @@ -115,7 +114,7 @@ function backward(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, x, y) where {lat θ = atan(x, ρₒ - y) ρ′ = sqrt(x^2 + (ρₒ - y)^2) α′ = (C - (ρ′^2 * n^2)) / n - β′ = asin(α′ / (1 - (1 - e) / (2 * e) * log((1 - e) / (1 + e)))) + β′ = asin(α′ / (1 - ((1 - e²) / (2 * e)) * log((1 - e) / (1 + e)))) λ = θ / n ϕ = auth2geod(β′, e²) @@ -129,7 +128,7 @@ end hm(ϕ, e, e²) = cos(ϕ) / sqrt(1 - e² * sin(ϕ)^2) -hα(ϕ, e, e²) = (1 - e²) * (sin(ϕ) / (1 - e² * sin(ϕ)^2) - (1 / (2 * e)) * log((1 - e * sin(ϕ)) / (1 + e * sin(ϕ)))) +hα(ϕ, e, e²) = (1 - e²) * ((sin(ϕ) / (1 - e² * sin(ϕ)^2)) - (1 / (2 * e) * log((1 - e * sin(ϕ)) / (1 + e * sin(ϕ))))) hλ(λ) = λ > π ? λ - 2π : λ < -π ? λ + 2π : λ diff --git a/test/crs/conversions.jl b/test/crs/conversions.jl index dacd54e4..c69c434c 100644 --- a/test/crs/conversions.jl +++ b/test/crs/conversions.jl @@ -1394,7 +1394,7 @@ # type stability c1 = LatLon{NAD83}(T(45), T(90)) c2 = AlbersUS(T(-7231430.540202629), T(11853758.709623523)) - @inferred convert(AlbersUS, c1) + @inferred convert(Albers{23°,29.5°,45.5°,NAD83}, c1) @inferred convert(LatLon{NAD83}, c2) end From 3d26c03a4a491147c5b36ef3ada160fab6031d4d Mon Sep 17 00:00:00 2001 From: Jeffrey Chandler Date: Mon, 28 Oct 2024 12:11:24 -0600 Subject: [PATCH 26/38] Applied suggestions from code review, fixed inbounds to check rho (will do more checks later today), and removed the default constructor function Co-authored-by: Elias Carvalho <73039601+eliascarv@users.noreply.github.com> --- src/crs/projected/eqareaalbers.jl | 37 +++++++++++++++++-------------- 1 file changed, 20 insertions(+), 17 deletions(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index d29babed..dd3ba6d2 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -3,21 +3,9 @@ # ------------------------------------------------------------------ """ - Albers{latₒ, lat₁, lat₂, Datum, Shift} + Albers{latₒ,lat₁,lat₂,Datum,Shift} -Albers CRS with latitude origin latₒ standard parallels `lat₁` and `lat₂`, `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). +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 @@ -35,8 +23,6 @@ Albers{latₒ,lat₁,lat₂,Datum,Shift}(x::Number, y::Number) where {latₒ,lat Albers{latₒ,lat₁,lat₂,Datum}(args...) where {latₒ,lat₁,lat₂,Datum} = Albers{latₒ,lat₁,lat₂,Datum,Shift()}(args...) -Albers(args...) = Albers{NAD83}(args...) - Base.convert( ::Type{Albers{latₒ,lat₁,lat₂,Datum,Shift,M}}, coords::Albers{latₒ,lat₁,lat₂,Datum,Shift} @@ -63,7 +49,24 @@ lentype(::Type{<:Albers{latₒ,lat₁,lat₂,Datum,Shift,M}}) where {latₒ,lat # 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}, λ, ϕ) = -2π ≤ λ ≤ 2π && -π ≤ ϕ ≤ π +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₁))) + ϕ₂ = oftype(λ, ustrip(deg2rad(lat₂))) + + m₁ = hm(ϕ₁, e, e²) + m₂ = hm(ϕ₂, e, e²) + α₁ = hα(ϕ₁, e, e²) + α₂ = hα(ϕ₂, e, e²) + n = (m₁^2 - m₂^2) / (α₂ - α₁) + C = m₁^2 + n * α₁ + + ρ = sqrt(C - n * hα(ϕ, e, e²)) / n + ρ ≥ 0 +end function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where {latₒ,lat₁,lat₂,Datum,T} 🌎 = ellipsoid(Datum) From 42b57de797acde6fe28e5cec2b64c3764928f057 Mon Sep 17 00:00:00 2001 From: souma4 Date: Mon, 28 Oct 2024 12:16:00 -0600 Subject: [PATCH 27/38] Changed h functions to _albers to avoid potential namespace conflicts. Updated tests to use the results from Proj.jl and EPSG:4269 (previously used general Proj with EPSG:4326). changed \Theta's to \theta to reflect documentation. updated backwards \theta to switch signs of arguments if n <0 based on latest EPSG guidance note. --- src/crs/projected/eqareaalbers.jl | 37 ++++++++++++++++--------------- test/crs/conversions.jl | 10 ++++----- 2 files changed, 24 insertions(+), 23 deletions(-) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index d29babed..79291478 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -73,20 +73,20 @@ function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where ϕ₁ = T(ustrip(deg2rad(lat₁))) ϕ₂ = T(ustrip(deg2rad(lat₂))) - m₁ = hm(ϕ₁, e, e²) - m₂ = hm(ϕ₂, e, e²) - α₁ = hα(ϕ₁, e, e²) - α₂ = hα(ϕ₂, e, e²) + m₁ = _albersm(ϕ₁, e, e²) + m₂ = _albersm(ϕ₂, e, e²) + α₁ = _albersα(ϕ₁, e, e²) + α₂ = _albersα(ϕ₂, e, e²) n = (m₁^2 - m₂^2) / (α₂ - α₁) C = m₁^2 + n * α₁ - Θ(λ) = n * λ - ρ(ϕ) = sqrt(C - n * hα(ϕ, e, e²)) / n + θ(λ) = n * λ + ρ(ϕ) = sqrt(C - n * _albersα(ϕ, e, e²)) / n ρₒ = ρ(ϕₒ) - fx(λ, ϕ) = ρ(hϕ(ϕ)) * sin(Θ(hλ(λ))) + fx(λ, ϕ) = ρ(_albersϕ(ϕ)) * sin(θ(_albersλ(λ))) - fy(λ, ϕ) = ρₒ - ρ(hϕ(ϕ)) * cos(Θ(hλ(λ))) + fy(λ, ϕ) = ρₒ - ρ(_albersϕ(ϕ)) * cos(θ(_albersλ(λ))) fx, fy end @@ -101,17 +101,17 @@ function backward(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, x, y) where {lat ϕ₁ = oftype(x, ustrip(deg2rad(lat₁))) ϕ₂ = oftype(x, ustrip(deg2rad(lat₂))) - m₁ = hm(ϕ₁, e, e²) - m₂ = hm(ϕ₂, e, e²) - α₁ = hα(ϕ₁, e, e²) - α₂ = hα(ϕ₂, e, e²) + m₁ = _albersm(ϕ₁, e, e²) + m₂ = _albersm(ϕ₂, e, e²) + α₁ = _albersα(ϕ₁, e, e²) + α₂ = _albersα(ϕ₂, e, e²) n = (m₁^2 - m₂^2) / (α₂ - α₁) C = m₁^2 + n * α₁ - ρ(ϕ) = sqrt(C - n * hα(ϕ, e, e²)) / n + ρ(ϕ) = sqrt(C - n * _albersα(ϕ, e, e²)) / n ρₒ = ρ(ϕₒ) - θ = atan(x, ρₒ - y) + θ = 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)))) @@ -126,13 +126,14 @@ end # HELPER FUNCTIONS # ----------------- -hm(ϕ, e, e²) = cos(ϕ) / sqrt(1 - e² * sin(ϕ)^2) +_albersm(ϕ, e, e²) = cos(ϕ) / sqrt(1 - e² * sin(ϕ)^2) -hα(ϕ, e, e²) = (1 - e²) * ((sin(ϕ) / (1 - e² * sin(ϕ)^2)) - (1 / (2 * e) * log((1 - e * sin(ϕ)) / (1 + e * sin(ϕ))))) +_albersα(ϕ, e, e²) = + (1 - e²) * ((sin(ϕ) / (1 - e² * sin(ϕ)^2)) - (1 / (2 * e) * log((1 - e * sin(ϕ)) / (1 + e * sin(ϕ))))) -hλ(λ) = λ > π ? λ - 2π : λ < -π ? λ + 2π : λ +_albersλ(λ) = λ > π ? λ - 2π : λ < -π ? λ + 2π : λ -hϕ(ϕ) = ϕ > π / 2 ? ϕ - π : ϕ < -π / 2 ? ϕ + π : ϕ +_albersϕ(ϕ) = ϕ > π / 2 ? ϕ - π : ϕ < -π / 2 ? ϕ + π : ϕ # ---------- # FALLBACKS # ---------- diff --git a/test/crs/conversions.jl b/test/crs/conversions.jl index c69c434c..0306a015 100644 --- a/test/crs/conversions.jl +++ b/test/crs/conversions.jl @@ -1369,31 +1369,31 @@ 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(-7231430.540202629), T(11853758.709623523))) + @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(15156419.949174397), T(13963188.032694416))) + @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.12299166346), T(2460630.2617624514))) + @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(8640856.920871278), -T(4596167.748123343))) + @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(-7231430.540202629), T(11853758.709623523)) + c2 = AlbersUS(-T(7.231430540202629e6), T(1.1853758709623523e7)) @inferred convert(Albers{23°,29.5°,45.5°,NAD83}, c1) @inferred convert(LatLon{NAD83}, c2) end From ea4d3c607c8fce0b351b719c06999eb8b1c2157e Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Mon, 28 Oct 2024 17:07:52 -0300 Subject: [PATCH 28/38] =?UTF-8?q?Remove=20'=5Falbers=CE=BB'=20and=20'=5Fal?= =?UTF-8?q?bers=CE=BB'=20by=20fixing=20the=20shifted=20longitude?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/crs/projected.jl | 2 +- src/crs/projected/eqareaalbers.jl | 26 ++++++++++---------------- 2 files changed, 11 insertions(+), 17 deletions(-) diff --git a/src/crs/projected.jl b/src/crs/projected.jl index 21f3b018..090b5274 100644 --- a/src/crs/projected.jl +++ b/src/crs/projected.jl @@ -142,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(fixlon(deg2rad(coords.lon - lonₒ))) ϕ = ustrip(deg2rad(coords.lat)) if !inbounds(C, λ, ϕ) throw(ArgumentError("coordinates outside of the projection domain")) diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/eqareaalbers.jl index a4cf0b21..35b608ce 100644 --- a/src/crs/projected/eqareaalbers.jl +++ b/src/crs/projected/eqareaalbers.jl @@ -53,18 +53,17 @@ function inbounds(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, λ, ϕ) where {l 🌎 = ellipsoid(Datum) e = oftype(λ, eccentricity(🌎)) e² = oftype(λ, eccentricity²(🌎)) - ϕₒ = oftype(λ, ustrip(deg2rad(latₒ))) ϕ₁ = oftype(λ, ustrip(deg2rad(lat₁))) ϕ₂ = oftype(λ, ustrip(deg2rad(lat₂))) - - m₁ = hm(ϕ₁, e, e²) - m₂ = hm(ϕ₂, e, e²) - α₁ = hα(ϕ₁, e, e²) - α₂ = hα(ϕ₂, e, e²) + + m₁ = _albersm(ϕ₁, e, e²) + m₂ = _albersm(ϕ₂, e, e²) + α₁ = _albersα(ϕ₁, e, e²) + α₂ = _albersα(ϕ₂, e, e²) n = (m₁^2 - m₂^2) / (α₂ - α₁) C = m₁^2 + n * α₁ - - ρ = sqrt(C - n * hα(ϕ, e, e²)) / n + + ρ = sqrt(C - n * _albersα(ϕ, e, e²)) / n ρ ≥ 0 end @@ -87,15 +86,13 @@ function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where ρ(ϕ) = sqrt(C - n * _albersα(ϕ, e, e²)) / n ρₒ = ρ(ϕₒ) - fx(λ, ϕ) = ρ(_albersϕ(ϕ)) * sin(θ(_albersλ(λ))) + fx(λ, ϕ) = ρ(ϕ) * sin(θ(λ)) - fy(λ, ϕ) = ρₒ - ρ(_albersϕ(ϕ)) * cos(θ(_albersλ(λ))) + fy(λ, ϕ) = ρₒ - ρ(ϕ) * cos(θ(λ)) fx, fy end -# backward projection formulas - function backward(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, x, y) where {latₒ,lat₁,lat₂,Datum} 🌎 = ellipsoid(Datum) e = oftype(x, eccentricity(🌎)) @@ -114,7 +111,7 @@ function backward(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, x, y) where {lat ρ(ϕ) = sqrt(C - n * _albersα(ϕ, e, e²)) / n ρₒ = ρ(ϕₒ) - θ = n >= 0 ? atan(x, ρₒ - y) : atan(-x, y - ρₒ) + θ = 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)))) @@ -134,9 +131,6 @@ _albersm(ϕ, e, 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(ϕ))))) -_albersλ(λ) = λ > π ? λ - 2π : λ < -π ? λ + 2π : λ - -_albersϕ(ϕ) = ϕ > π / 2 ? ϕ - π : ϕ < -π / 2 ? ϕ + π : ϕ # ---------- # FALLBACKS # ---------- From 8e2c9b85ebc0751acf991e6996e24123da075da2 Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Mon, 28 Oct 2024 17:12:33 -0300 Subject: [PATCH 29/38] Rename the file --- src/crs/projected.jl | 2 +- src/crs/projected/{eqareaalbers.jl => albers.jl} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename src/crs/projected/{eqareaalbers.jl => albers.jl} (100%) diff --git a/src/crs/projected.jl b/src/crs/projected.jl index 090b5274..b76e00c4 100644 --- a/src/crs/projected.jl +++ b/src/crs/projected.jl @@ -111,7 +111,7 @@ include("projected/winkeltripel.jl") include("projected/robinson.jl") include("projected/orthographic.jl") include("projected/transversemercator.jl") -include("projected/eqareaalbers.jl") +include("projected/albers.jl") # ---------- # FALLBACKS diff --git a/src/crs/projected/eqareaalbers.jl b/src/crs/projected/albers.jl similarity index 100% rename from src/crs/projected/eqareaalbers.jl rename to src/crs/projected/albers.jl From 26cf1c7b17e2ed762497edfa0c8e5052baade849 Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Mon, 28 Oct 2024 17:28:45 -0300 Subject: [PATCH 30/38] Add more helpers to avoid code repetition --- src/crs/projected/albers.jl | 47 +++++++++++++++++-------------------- test/crs/conversions.jl | 2 +- 2 files changed, 23 insertions(+), 26 deletions(-) diff --git a/src/crs/projected/albers.jl b/src/crs/projected/albers.jl index 35b608ce..84311e48 100644 --- a/src/crs/projected/albers.jl +++ b/src/crs/projected/albers.jl @@ -56,14 +56,8 @@ function inbounds(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, λ, ϕ) where {l ϕ₁ = oftype(λ, ustrip(deg2rad(lat₁))) ϕ₂ = oftype(λ, ustrip(deg2rad(lat₂))) - m₁ = _albersm(ϕ₁, e, e²) - m₂ = _albersm(ϕ₂, e, e²) - α₁ = _albersα(ϕ₁, e, e²) - α₂ = _albersα(ϕ₂, e, e²) - n = (m₁^2 - m₂^2) / (α₂ - α₁) - C = m₁^2 + n * α₁ - - ρ = sqrt(C - n * _albersα(ϕ, e, e²)) / n + C, n = _ambersCn(ϕ₁, ϕ₂, e, e²) + ρ = _ambersρ(ϕ, C, n, e, e²) ρ ≥ 0 end @@ -75,15 +69,10 @@ function formulas(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, ::Type{T}) where ϕ₁ = T(ustrip(deg2rad(lat₁))) ϕ₂ = T(ustrip(deg2rad(lat₂))) - m₁ = _albersm(ϕ₁, e, e²) - m₂ = _albersm(ϕ₂, e, e²) - α₁ = _albersα(ϕ₁, e, e²) - α₂ = _albersα(ϕ₂, e, e²) - n = (m₁^2 - m₂^2) / (α₂ - α₁) - C = m₁^2 + n * α₁ + C, n = _ambersCn(ϕ₁, ϕ₂, e, e²) θ(λ) = n * λ - ρ(ϕ) = sqrt(C - n * _albersα(ϕ, e, e²)) / n + ρ(ϕ) = _ambersρ(ϕ, C, n, e, e²) ρₒ = ρ(ϕₒ) fx(λ, ϕ) = ρ(ϕ) * sin(θ(λ)) @@ -101,15 +90,8 @@ function backward(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, x, y) where {lat ϕ₁ = oftype(x, ustrip(deg2rad(lat₁))) ϕ₂ = oftype(x, ustrip(deg2rad(lat₂))) - m₁ = _albersm(ϕ₁, e, e²) - m₂ = _albersm(ϕ₂, e, e²) - α₁ = _albersα(ϕ₁, e, e²) - α₂ = _albersα(ϕ₂, e, e²) - n = (m₁^2 - m₂^2) / (α₂ - α₁) - C = m₁^2 + n * α₁ - - ρ(ϕ) = sqrt(C - n * _albersα(ϕ, e, e²)) / n - ρₒ = ρ(ϕₒ) + C, n = _ambersCn(ϕ₁, ϕ₂, e, e²) + ρₒ = _ambersρ(ϕₒ, C, n, e, e²) θ = n ≥ 0 ? atan(x, ρₒ - y) : atan(-x, y - ρₒ) ρ′ = sqrt(x^2 + (ρₒ - y)^2) @@ -126,7 +108,22 @@ end # HELPER FUNCTIONS # ----------------- -_albersm(ϕ, e, e²) = cos(ϕ) / sqrt(1 - e² * sin(ϕ)^2) +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(ϕ))))) diff --git a/test/crs/conversions.jl b/test/crs/conversions.jl index 0306a015..d27097ce 100644 --- a/test/crs/conversions.jl +++ b/test/crs/conversions.jl @@ -1394,7 +1394,7 @@ # type stability c1 = LatLon{NAD83}(T(45), T(90)) c2 = AlbersUS(-T(7.231430540202629e6), T(1.1853758709623523e7)) - @inferred convert(Albers{23°,29.5°,45.5°,NAD83}, c1) + @inferred convert(AlbersUS, c1) @inferred convert(LatLon{NAD83}, c2) end From 9bfbbaff388fb628dc9800b4b8bf8f4c0ede460a Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Mon, 28 Oct 2024 17:45:22 -0300 Subject: [PATCH 31/38] Fix code --- src/crs/projected.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/crs/projected.jl b/src/crs/projected.jl index b76e00c4..7f47a255 100644 --- a/src/crs/projected.jl +++ b/src/crs/projected.jl @@ -142,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(fixlon(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")) From b57badeb8bb92b4e47bb842b2a16f65443f116b6 Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Mon, 28 Oct 2024 17:50:01 -0300 Subject: [PATCH 32/38] Fix test --- test/crs/mactype.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/crs/mactype.jl b/test/crs/mactype.jl index 0cd6f423..50ae3228 100644 --- a/test/crs/mactype.jl +++ b/test/crs/mactype.jl @@ -142,7 +142,7 @@ @test c2 isa C AlbersUS = CoordRefSystems.shift(Albers{23.0°,29.5°,45.5°,NAD83}, lonₒ=-96.0°) - C = Albers{Met{T}} + C = AlbersUS{Met{T}} c1 = AlbersUS(1.0, 1.0) c2 = convert(C, c1) @test c2 isa C From 2f866d12256cd2b94e9489cfe2d7c38cf1436c77 Mon Sep 17 00:00:00 2001 From: souma4 Date: Mon, 28 Oct 2024 14:59:12 -0600 Subject: [PATCH 33/38] Added WKT strings to strings.jl --- src/strings.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/strings.jl b/src/strings.jl index 5e0b0613..d841807c 100644 --- a/src/strings.jl +++ b/src/strings.jl @@ -24,7 +24,9 @@ const esriid2code = Dict( "World_Behrmann" => ESRI{54017}, "World_Cylindrical_Equal_Area" => ESRI{54034}, "World_Robinson" => ESRI{54030}, - "World_Winkel_Tripel_NGS" => ESRI{54042} + "World_Winkel_Tripel_NGS" => ESRI{54042}, + "NAD83 / Conus Albers" => EPSG{5070}, + "NAD_1983_Contiguous_USA_Albers" => EPSG{5070} ) """ @@ -50,7 +52,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] From 1b19bf6b9102c0ed9838222fa8698cbe99591713 Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Mon, 28 Oct 2024 18:04:35 -0300 Subject: [PATCH 34/38] Add missing construtor & Update tests --- src/crs/projected/albers.jl | 2 ++ test/crs/constructors.jl | 8 ++++---- test/crs/domains.jl | 2 +- test/crs/rand.jl | 3 ++- 4 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/crs/projected/albers.jl b/src/crs/projected/albers.jl index 84311e48..44c7594d 100644 --- a/src/crs/projected/albers.jl +++ b/src/crs/projected/albers.jl @@ -23,6 +23,8 @@ Albers{latₒ,lat₁,lat₂,Datum,Shift}(x::Number, y::Number) where {latₒ,lat 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} diff --git a/test/crs/constructors.jl b/test/crs/constructors.jl index 09bb3cc1..a03bbe7a 100644 --- a/test/crs/constructors.jl +++ b/test/crs/constructors.jl @@ -625,21 +625,21 @@ end @testset "Albers" begin - A = Albers{23.0°,29.5°,45.5°,NAD83} + 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{NAD83}(x: 1.0 m, y: 1.0 m)" + @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{NAD83} coordinates + Albers{WGS84Latest} coordinates ├─ x: 1.0f0 m └─ y: 1.0f0 m""" else @test sprint(show, MIME("text/plain"), c) == """ - Albers{NAD83} coordinates + Albers{WGS84Latest} coordinates ├─ x: 1.0 m └─ y: 1.0 m""" end diff --git a/test/crs/domains.jl b/test/crs/domains.jl index 11c813a9..a7d90ec4 100644 --- a/test/crs/domains.jl +++ b/test/crs/domains.jl @@ -214,7 +214,7 @@ @testset "Albers Forward" begin 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(lat, lon) + c1 = LatLon{NAD83}(lat, lon) if indomain(AlbersUS, c1) c2 = convert(AlbersUS, c1) @test isfinite(c2.x) diff --git a/test/crs/rand.jl b/test/crs/rand.jl index 3728eeff..aca5698c 100644 --- a/test/crs/rand.jl +++ b/test/crs/rand.jl @@ -70,6 +70,7 @@ randtest(OrthoSouth{WGS84Latest}) randtest(OrthoSouth) - randtest(CoordRefSystems.get(EPSG{5070})) + randtest(Albers{23.0°,29.5°,45.5°,WGS84Latest}) + randtest(Albers{23.0°,29.5°,45.5°}) end end From 29dec859ab5cf0ac2a08671261ea27610ef35282 Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Mon, 28 Oct 2024 18:58:01 -0300 Subject: [PATCH 35/38] Add 'asinclamp' helper function --- src/crs/projected/albers.jl | 2 +- src/crs/projected/orthographic.jl | 5 ++--- src/utils.jl | 7 +++++++ 3 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/crs/projected/albers.jl b/src/crs/projected/albers.jl index 44c7594d..a21e8443 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 88260644..f91343b9 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 be037f62..e2f04cdc 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) From 5edb01380861191b4bff2875ecf5f45b5ad329d7 Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Tue, 29 Oct 2024 08:10:19 -0300 Subject: [PATCH 36/38] Add missing tests --- test/crs/domains.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/test/crs/domains.jl b/test/crs/domains.jl index a7d90ec4..f61fe867 100644 --- a/test/crs/domains.jl +++ b/test/crs/domains.jl @@ -211,7 +211,8 @@ end end - @testset "Albers Forward" begin + @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) @@ -219,6 +220,8 @@ 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 From 16c72889c2c1779640161aac4a807439a16bf931 Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Tue, 29 Oct 2024 08:24:41 -0300 Subject: [PATCH 37/38] Apply suggestions --- src/strings.jl | 5 ++--- test/crs/conversions.jl | 2 ++ 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/strings.jl b/src/strings.jl index d841807c..0b46b390 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}, @@ -24,9 +25,7 @@ const esriid2code = Dict( "World_Behrmann" => ESRI{54017}, "World_Cylindrical_Equal_Area" => ESRI{54034}, "World_Robinson" => ESRI{54030}, - "World_Winkel_Tripel_NGS" => ESRI{54042}, - "NAD83 / Conus Albers" => EPSG{5070}, - "NAD_1983_Contiguous_USA_Albers" => EPSG{5070} + "World_Winkel_Tripel_NGS" => ESRI{54042} ) """ diff --git a/test/crs/conversions.jl b/test/crs/conversions.jl index d27097ce..b5505a71 100644 --- a/test/crs/conversions.jl +++ b/test/crs/conversions.jl @@ -1366,6 +1366,8 @@ end @testset "LatLon <> Albers" begin + # PROJ transformation used: + # 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) From 3357c821fc1781502a527619801111fd908d6327 Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Tue, 29 Oct 2024 08:26:56 -0300 Subject: [PATCH 38/38] Update comments --- test/crs/conversions.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/crs/conversions.jl b/test/crs/conversions.jl index b5505a71..d4add469 100644 --- a/test/crs/conversions.jl +++ b/test/crs/conversions.jl @@ -1366,8 +1366,7 @@ end @testset "LatLon <> Albers" begin - # PROJ transformation used: - # Proj.Transformation("EPSG:4269", "EPSG:5070") + # 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)