Add EquidistantCylindrical projection (#751)

* Add EquidistantCylindrical projection * Update src/crs.jl Co-authored-by: Júlio Hoffimann <julio.hoffimann@gmail.com> * Update src/crs.jl * Rename file --------- Co-authored-by: Júlio Hoffimann <julio.hoffimann@gmail.com>
JuliaGeometry · Feb 9, 2024 · bc8aae3 · bc8aae3
1 parent 3943bbc
commit bc8aae3
Show file tree

Hide file tree

Showing 8 changed files with 88 additions and 67 deletions.
diff --git a/src/crs.jl b/src/crs.jl
@@ -36,19 +36,20 @@ function typealias end
 # ----
 
 """
-    CRS{ID,Coords,Datum}
+    CRS{ID,Coords,Datum,Params}
 
 A Coordinate Reference System (CRS) with identifier `ID`, coordinates `Coords`
-and given `Datum` can be used to georeference any point in physical space.
+a given `Datum` and parameters `Params` can be used to georeference any point
+in physical space.
 
 The `CRS` type is not intended for end-users. Aliases are provided,
 such as `LatLon` and `WebMercator`, to facilitate coordinate system conversions.
 """
-struct CRS{ID,Coords,Datum}
+struct CRS{ID,Coords,Datum,Params}
   coords::Coords
 end
 
-CRS{ID,Coords,Datum}(args...) where {ID,Coords,Datum} = CRS{ID,Coords,Datum}(Coords(args))
+CRS{ID,Coords,Datum,Params}(args...) where {ID,Coords,Datum,Params} = CRS{ID,Coords,Datum,Params}(Coords(args))
 
 _coords(coords::CRS) = getfield(coords, :coords)
 
@@ -72,7 +73,7 @@ end
 
 Returns the datum of the coordinates `coords`.
 """
-datum(::CRS{ID,Coords,Datum}) where {ID,Coords,Datum} = Datum
+datum(::CRS{ID,Coords,Datum,Params}) where {ID,Coords,Datum,Params} = Datum
 
 """
     ellipsoid(coords)
@@ -130,11 +131,13 @@ const Met{T} = Quantity{T,u"𝐋",typeof(u"m")}
 const Rad{T} = Quantity{T,NoDims,typeof(u"rad")}
 const Deg{T} = Quantity{T,NoDims,typeof(u"°")}
 
+const NoParams = nothing
+
 include("crs/basic.jl")
 include("crs/latlon.jl")
 include("crs/mercator.jl")
 include("crs/webmercator.jl")
-include("crs/platecarree.jl")
+include("crs/eqdistcylindrical.jl")
 
 # ----------
 # FALLBACKS

diff --git a/src/crs/basic.jl b/src/crs/basic.jl
@@ -22,7 +22,7 @@ Cartesian(1.0u"km", 1.0u"km", 1.0u"km")
 * [ISO 80000-2:2019](https://www.iso.org/standard/64973.html)
 * [ISO 80000-3:2019](https://www.iso.org/standard/64974.html)
 """
-const Cartesian{N,L<:Len} = CRS{:Cartesian,NTuple{N,L},NoDatum}
+const Cartesian{N,L<:Len} = CRS{:Cartesian,NTuple{N,L},NoDatum,NoParams}
 
 Cartesian(coords::Vararg{L,N}) where {N,L<:Len} = Cartesian{N,float(L)}(coords)
 Cartesian(coords::Len...) = Cartesian(promote(coords...)...)
@@ -76,7 +76,7 @@ Polar(1.0u"km", (π/4)u"rad")
 * [ISO 80000-2:2019](https://www.iso.org/standard/64973.html)
 * [ISO 80000-3:2019](https://www.iso.org/standard/64974.html)
 """
-const Polar{L<:Len,R<:Rad} = CRS{:Polar,@NamedTuple{ρ::L, ϕ::R},NoDatum}
+const Polar{L<:Len,R<:Rad} = CRS{:Polar,@NamedTuple{ρ::L, ϕ::R},NoDatum,NoParams}
 
 Polar(ρ::L, ϕ::R) where {L<:Len,R<:Rad} = Polar{float(L),float(R)}(ρ, ϕ)
 Polar(ρ::Len, ϕ::Deg) = Polar(ρ, deg2rad(ϕ))
@@ -104,7 +104,7 @@ Cylindrical(1.0u"km", (π/4)u"rad", 1.0u"km")
 * [ISO 80000-2:2019](https://www.iso.org/standard/64973.html)
 * [ISO 80000-3:2019](https://www.iso.org/standard/64974.html)
 """
-const Cylindrical{L<:Len,R<:Rad} = CRS{:Cylindrical,@NamedTuple{ρ::L, ϕ::R, z::L},NoDatum}
+const Cylindrical{L<:Len,R<:Rad} = CRS{:Cylindrical,@NamedTuple{ρ::L, ϕ::R, z::L},NoDatum,NoParams}
 
 Cylindrical(ρ::L, ϕ::R, z::L) where {L<:Len,R<:Rad} = Cylindrical{float(L),float(R)}(ρ, ϕ, z)
 function Cylindrical(ρ::Len, ϕ::Rad, z::Len)
@@ -135,7 +135,7 @@ Spherical(1.0u"km", (π/4)u"rad", (π/4)u"rad")
 * [ISO 80000-2:2019](https://www.iso.org/standard/64973.html)
 * [ISO 80000-3:2019](https://www.iso.org/standard/64974.html)
 """
-const Spherical{L<:Len,R<:Rad} = CRS{:Spherical,@NamedTuple{r::L, θ::R, ϕ::R},NoDatum}
+const Spherical{L<:Len,R<:Rad} = CRS{:Spherical,@NamedTuple{r::L, θ::R, ϕ::R},NoDatum,NoParams}
 
 Spherical(r::L, θ::R, ϕ::R) where {L<:Len,R<:Rad} = Spherical{float(L),float(R)}(r, θ, ϕ)
 Spherical(r::Len, θ::Rad, ϕ::Rad) = Spherical(r, promote(θ, ϕ)...)

diff --git a/src/crs/eqdistcylindrical.jl b/src/crs/eqdistcylindrical.jl
@@ -0,0 +1,65 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    EquidistantCylindrical{ID,latₜₛ}
+
+Equidistant Cylindrical CRS with identifier `ID` and latitude of true scale `latₜₛ` in degrees.
+"""
+const EquidistantCylindrical{ID,latₜₛ,M<:Met} = CRS{ID,@NamedTuple{x::M, y::M},WGS84,latₜₛ}
+
+EquidistantCylindrical{ID,latₜₛ}(x::M, y::M) where {ID,latₜₛ,M<:Met} = EquidistantCylindrical{ID,latₜₛ,float(M)}(x, y)
+EquidistantCylindrical{ID,latₜₛ}(x::Met, y::Met) where {ID,latₜₛ} = EquidistantCylindrical{ID,latₜₛ}(promote(x, y)...)
+EquidistantCylindrical{ID,latₜₛ}(x::Len, y::Len) where {ID,latₜₛ} =
+  EquidistantCylindrical{ID,latₜₛ}(uconvert(u"m", x), uconvert(u"m", y))
+EquidistantCylindrical{ID,latₜₛ}(x::Number, y::Number) where {ID,latₜₛ} =
+  EquidistantCylindrical{ID,latₜₛ}(addunit(x, u"m"), addunit(y, u"m"))
+
+"""
+    PlateCarree(x, y)
+
+Plate Carrée coordinates in length units (default to meter).
+
+## Examples
+
+```julia
+PlateCarree(1, 1) # add default units
+PlateCarree(1u"m", 1u"m") # integers are converted converted to floats
+PlateCarree(1.0u"km", 1.0u"km") # length quantities are converted to meters
+PlateCarree(1.0u"m", 1.0u"m")
+```
+
+See [EPSG:32662](https://epsg.io/32662).
+"""
+const PlateCarree = EquidistantCylindrical{EPSG{32662},0.0u"°"}
+
+typealias(::Type{EPSG{32662}}) = PlateCarree
+
+function Base.convert(::Type{EquidistantCylindrical{ID,latₜₛ}}, coords::LatLon) where {ID,latₜₛ}
+  🌎 = ellipsoid(coords)
+  λ = deg2rad(coords.lon)
+  ϕ = deg2rad(coords.lat)
+  ϕₜₛ = oftype(ϕ, deg2rad(latₜₛ))
+  l = ustrip(λ)
+  o = ustrip(ϕ)
+  a = oftype(l, ustrip(majoraxis(🌎)))
+
+  x = a * l * cos(ϕₜₛ)
+  y = a * o
+
+  EquidistantCylindrical{ID,latₜₛ}(x * u"m", y * u"m")
+end
+
+function Base.convert(::Type{LatLon}, coords::EquidistantCylindrical{ID,latₜₛ}) where {ID,latₜₛ}
+  🌎 = ellipsoid(coords)
+  x = coords.x
+  y = coords.y
+  a = oftype(x, majoraxis(🌎))
+  ϕₜₛ = numconvert(numtype(x), deg2rad(latₜₛ))
+
+  λ = x / (cos(ϕₜₛ) * a)
+  ϕ = y / a
+
+  LatLon(rad2deg(ϕ) * u"°", rad2deg(λ) * u"°")
+end
diff --git a/src/crs/latlon.jl b/src/crs/latlon.jl
@@ -18,7 +18,7 @@ LatLon(45.0u"°", 45.0u"°")
 
 See [EPSG:4326](https://epsg.io/4326).
 """
-const LatLon{D<:Deg} = CRS{EPSG{4326},@NamedTuple{lat::D, lon::D},WGS84}
+const LatLon{D<:Deg} = CRS{EPSG{4326},@NamedTuple{lat::D, lon::D},WGS84,NoParams}
 
 typealias(::Type{EPSG{4326}}) = LatLon
 

diff --git a/src/crs/mercator.jl b/src/crs/mercator.jl
@@ -18,7 +18,7 @@ Mercator(1.0u"m", 1.0u"m")
 
 See [EPSG:3395](https://epsg.io/3395).
 """
-const Mercator{M<:Met} = CRS{EPSG{3395},@NamedTuple{x::M, y::M},WGS84}
+const Mercator{M<:Met} = CRS{EPSG{3395},@NamedTuple{x::M, y::M},WGS84,NoParams}
 
 typealias(::Type{EPSG{3395}}) = Mercator
 

diff --git a/src/crs/platecarree.jl b/src/crs/platecarree.jl
diff --git a/src/crs/webmercator.jl b/src/crs/webmercator.jl
@@ -18,7 +18,7 @@ WebMercator(1.0u"m", 1.0u"m")
 
 See [EPSG:3857](https://epsg.io/3857).
 """
-const WebMercator{M<:Met} = CRS{EPSG{3857},@NamedTuple{x::M, y::M},WGS84}
+const WebMercator{M<:Met} = CRS{EPSG{3857},@NamedTuple{x::M, y::M},WGS84,NoParams}
 
 typealias(::Type{EPSG{3857}}) = WebMercator
 

diff --git a/src/utils.jl b/src/utils.jl
@@ -153,3 +153,10 @@ function atanpos(y, x)
   α = atan(y, x)
   ifelse(α ≥ zero(α), α, α + oftype(α, 2π))
 end
+
+"""
+    numconvert(T, x)
+
+Converts the number type of quantity `x` to `T`.
+"""
+numconvert(::Type{T}, x::Quantity{S,D,U}) where {T,S,D,U} = convert(Quantity{T,D,U}, x)