From b343506b756050c9c0f7c6d49f0ca680bac4739c Mon Sep 17 00:00:00 2001
From: Elias Carvalho <73039601+eliascarv@users.noreply.github.com>
Date: Fri, 9 Feb 2024 18:12:29 -0300
Subject: [PATCH] Add 'WinkelTripel' projection (#752)

---
 src/Meshes.jl                |  2 ++
 src/crs.jl                   | 16 ++++++++----
 src/crs/eqdistcylindrical.jl |  2 +-
 src/crs/mercator.jl          |  2 +-
 src/crs/webmercator.jl       |  2 +-
 src/crs/winkeltripel.jl      | 49 ++++++++++++++++++++++++++++++++++++
 src/datums.jl                | 14 +++++++++++
 src/ellipsoids.jl            |  9 +++++++
 test/crs.jl                  | 28 +++++++++++++++++++++
 test/datums.jl               | 14 +++++++++++
 10 files changed, 130 insertions(+), 8 deletions(-)
 create mode 100644 src/crs/winkeltripel.jl

diff --git a/src/Meshes.jl b/src/Meshes.jl
index bc56bdcb0..92ac30bc9 100644
--- a/src/Meshes.jl
+++ b/src/Meshes.jl
@@ -141,6 +141,7 @@ export
   Datum,
   NoDatum,
   WGS84,
+  WIII,
   ellipsoid,
   latitudeₒ,
   longitudeₒ,
@@ -158,6 +159,7 @@ export
   Mercator,
   WebMercator,
   PlateCarree,
+  WinkelTripel,
   datum,
 
   # geometries
diff --git a/src/crs.jl b/src/crs.jl
index c90a24407..29b14c96e 100644
--- a/src/crs.jl
+++ b/src/crs.jl
@@ -73,35 +73,40 @@ end
 
 Returns the datum of the coordinates `coords`.
 """
-datum(::CRS{ID,Coords,Datum,Params}) where {ID,Coords,Datum,Params} = Datum
+datum(coords::CRS) = datum(typeof(coords))
+datum(::Type{<:CRS{<:Any,<:Any,Datum}}) where {Datum} = Datum
 
 """
     ellipsoid(coords)
 
 Returns the ellipsoid of the coordinates `coords`.
 """
-ellipsoid(coords::CRS) = ellipsoid(datum(coords))
+ellipsoid(coords::CRS) = ellipsoid(typeof(coords))
+ellipsoid(T::Type{<:CRS}) = ellipsoid(datum(T))
 
 """
     latitudeₒ(coords)
 
 Returns the latitude origin of the coordinates `coords`.
 """
-latitudeₒ(coords::CRS) = latitudeₒ(datum(coords))
+latitudeₒ(coords::CRS) = latitudeₒ(typeof(coords))
+latitudeₒ(T::Type{<:CRS}) = latitudeₒ(datum(T))
 
 """
     longitudeₒ(coords)
 
 Returns the longitude origin of the coordinates `coords`.
 """
-longitudeₒ(coords::CRS) = longitudeₒ(datum(coords))
+longitudeₒ(coords::CRS) = longitudeₒ(typeof(coords))
+longitudeₒ(T::Type{<:CRS}) = longitudeₒ(datum(T))
 
 """
     altitudeₒ(coords)
 
 Returns the altitude origin of the coordinates `coords`.
 """
-altitudeₒ(coords::CRS) = altitudeₒ(datum(coords))
+altitudeₒ(coords::CRS) = altitudeₒ(typeof(coords))
+altitudeₒ(T::Type{<:CRS}) = altitudeₒ(datum(T))
 
 # -----------
 # IO METHODS
@@ -138,6 +143,7 @@ include("crs/latlon.jl")
 include("crs/mercator.jl")
 include("crs/webmercator.jl")
 include("crs/eqdistcylindrical.jl")
+include("crs/winkeltripel.jl")
 
 # ----------
 # FALLBACKS
diff --git a/src/crs/eqdistcylindrical.jl b/src/crs/eqdistcylindrical.jl
index 85029507f..f6e9b11f7 100644
--- a/src/crs/eqdistcylindrical.jl
+++ b/src/crs/eqdistcylindrical.jl
@@ -37,7 +37,7 @@ 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)
+  🌎 = ellipsoid(EquidistantCylindrical{ID,latₜₛ})
   λ = deg2rad(coords.lon)
   ϕ = deg2rad(coords.lat)
   ϕₜₛ = oftype(ϕ, deg2rad(latₜₛ))
diff --git a/src/crs/mercator.jl b/src/crs/mercator.jl
index 7f57cfb37..fe28c4162 100644
--- a/src/crs/mercator.jl
+++ b/src/crs/mercator.jl
@@ -32,7 +32,7 @@ Mercator(x::Number, y::Number) = Mercator(addunit(x, u"m"), addunit(y, u"m"))
 # ------------
 
 function Base.convert(::Type{Mercator}, coords::LatLon)
-  🌎 = ellipsoid(coords)
+  🌎 = ellipsoid(Mercator)
   λ = deg2rad(coords.lon)
   ϕ = deg2rad(coords.lat)
   l = ustrip(λ)
diff --git a/src/crs/webmercator.jl b/src/crs/webmercator.jl
index beadb96a6..14797d37a 100644
--- a/src/crs/webmercator.jl
+++ b/src/crs/webmercator.jl
@@ -32,7 +32,7 @@ WebMercator(x::Number, y::Number) = WebMercator(addunit(x, u"m"), addunit(y, u"m
 # ------------
 
 function Base.convert(::Type{WebMercator}, coords::LatLon)
-  🌎 = ellipsoid(coords)
+  🌎 = ellipsoid(WebMercator)
   λ = deg2rad(coords.lon)
   ϕ = deg2rad(coords.lat)
   l = ustrip(λ)
diff --git a/src/crs/winkeltripel.jl b/src/crs/winkeltripel.jl
new file mode 100644
index 000000000..d7420d690
--- /dev/null
+++ b/src/crs/winkeltripel.jl
@@ -0,0 +1,49 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+const Winkel{ID,lat₁,M<:Met} = CRS{ID,@NamedTuple{x::M, y::M},WIII,lat₁}
+
+Winkel{ID,lat₁}(x::M, y::M) where {ID,lat₁,M<:Met} = Winkel{ID,lat₁,float(M)}(x, y)
+Winkel{ID,lat₁}(x::Met, y::Met) where {ID,lat₁} = Winkel{ID,lat₁}(promote(x, y)...)
+Winkel{ID,lat₁}(x::Len, y::Len) where {ID,lat₁} = Winkel{ID,lat₁}(uconvert(u"m", x), uconvert(u"m", y))
+Winkel{ID,lat₁}(x::Number, y::Number) where {ID,lat₁} = Winkel{ID,lat₁}(addunit(x, u"m"), addunit(y, u"m"))
+
+"""
+    WinkelTripel(x, y)
+
+Winkel Tripel coordinates in length units (default to meter).
+
+## Examples
+
+```julia
+WinkelTripel(1, 1) # add default units
+WinkelTripel(1u"m", 1u"m") # integers are converted converted to floats
+WinkelTripel(1.0u"km", 1.0u"km") # length quantities are converted to meters
+WinkelTripel(1.0u"m", 1.0u"m")
+```
+
+See [ESRI:3395](https://epsg.io/53042).
+"""
+const WinkelTripel = Winkel{ESRI{53042},50.467u"°"}
+
+typealias(::Type{ESRI{53042}}) = WinkelTripel
+
+# ------------
+# CONVERSIONS
+# ------------
+
+function Base.convert(::Type{Winkel{ID,lat₁}}, coords::LatLon) where {ID,lat₁}
+  🌎 = ellipsoid(Winkel{ID,lat₁})
+  λ = deg2rad(coords.lon)
+  ϕ = deg2rad(coords.lat)
+  ϕ₁ = oftype(ϕ, deg2rad(lat₁))
+  l = ustrip(λ)
+  o = ustrip(ϕ)
+  a = oftype(l, ustrip(majoraxis(🌎)))
+  α = acos(cos(ϕ) * cos(λ / 2))
+  sincα = sin(α) / α # unnormalized sinc
+  x = a / 2 * (l * cos(ϕ₁) + (2cos(ϕ) * sin(λ / 2)) / sincα)
+  y = a / 2 * (o + sin(ϕ) / sincα)
+  Winkel{ID,lat₁}(x * u"m", y * u"m")
+end
diff --git a/src/datums.jl b/src/datums.jl
index 384a292d1..f3948bc61 100644
--- a/src/datums.jl
+++ b/src/datums.jl
@@ -62,3 +62,17 @@ ellipsoid(::Type{WGS84}) = WGS84🌎
 latitudeₒ(::Type{WGS84}) = 0.0u"°"
 longitudeₒ(::Type{WGS84}) = 0.0u"°"
 altitudeₒ(::Type{WGS84}) = 0.0u"m"
+
+"""
+    WGS84
+
+Winkel Tripel datum.
+
+See [ESRI:53042](https://epsg.io/53042)
+"""
+abstract type WIII <: Datum end
+
+ellipsoid(::Type{WIII}) = WIII🌎
+latitudeₒ(::Type{WIII}) = 0.0u"°"
+longitudeₒ(::Type{WIII}) = 0.0u"°"
+altitudeₒ(::Type{WIII}) = 0.0u"m"
diff --git a/src/ellipsoids.jl b/src/ellipsoids.jl
index 854d5e913..76577bccc 100644
--- a/src/ellipsoids.jl
+++ b/src/ellipsoids.jl
@@ -69,3 +69,12 @@ eccentricity(::Type{WGS84🌎}) = _WGS84.e
 eccentricity²(::Type{WGS84🌎}) = _WGS84.e²
 flattening(::Type{WGS84🌎}) = _WGS84.f
 flattening⁻¹(::Type{WGS84🌎}) = _WGS84.f⁻¹
+
+abstract type WIII🌎 <: RevolutionEllipsoid end
+
+majoraxis(::Type{WIII🌎}) = 6371000.0u"m"
+minoraxis(::Type{WIII🌎}) = 6371000.0u"m"
+eccentricity(::Type{WIII🌎}) = 0.0
+eccentricity²(::Type{WIII🌎}) = 0.0
+flattening(::Type{WIII🌎}) = 0.0
+flattening⁻¹(::Type{WIII🌎}) = Inf
diff --git a/test/crs.jl b/test/crs.jl
index 592fd11c6..88b0b4797 100644
--- a/test/crs.jl
+++ b/test/crs.jl
@@ -553,5 +553,33 @@
       @inferred convert(EPSG{32662}, c1)
       @inferred convert(EPSG{4326}, c2)
     end
+
+    @testset "LatLon <> WinkelTripel" begin
+      c1 = LatLon(T(45), T(90))
+      c2 = convert(WinkelTripel, c1)
+      @test c2 ≈ WinkelTripel(T(7036918.714302972), T(5225594.172156317))
+
+      c1 = LatLon(-T(45), T(90))
+      c2 = convert(WinkelTripel, c1)
+      @test c2 ≈ WinkelTripel(T(7036918.714302972), -T(5225594.172156317))
+
+      c1 = LatLon(T(45), -T(90))
+      c2 = convert(WinkelTripel, c1)
+      @test c2 ≈ WinkelTripel(-T(7036918.714302972), T(5225594.172156317))
+
+      c1 = LatLon(-T(45), -T(90))
+      c2 = convert(WinkelTripel, c1)
+      @test c2 ≈ WinkelTripel(-T(7036918.714302972), -T(5225594.172156317))
+
+      # EPSG fallback
+      c1 = LatLon(T(45), T(90))
+      c2 = convert(ESRI{53042}, c1)
+      @test c2 ≈ WinkelTripel(T(7036918.714302972), T(5225594.172156317))
+
+      # type stability
+      c1 = LatLon(T(45), T(90))
+      @inferred convert(WinkelTripel, c1)
+      @inferred convert(ESRI{53042}, c1)
+    end
   end
 end
diff --git a/test/datums.jl b/test/datums.jl
index d17feddb1..94b637022 100644
--- a/test/datums.jl
+++ b/test/datums.jl
@@ -19,4 +19,18 @@
     @test longitudeₒ(WGS84) == 0.0u"°"
     @test altitudeₒ(WGS84) == 0.0u"m"
   end
+
+  @testset "WIII" begin
+    🌎 = ellipsoid(WIII)
+    @test majoraxis(🌎) == 6371000.0u"m"
+    @test minoraxis(🌎) == 6371000.0u"m"
+    @test eccentricity(🌎) == 0.0
+    @test eccentricity²(🌎) == 0.0
+    @test flattening(🌎) == 0.0
+    @test flattening⁻¹(🌎) == Inf
+
+    @test latitudeₒ(WIII) == 0.0u"°"
+    @test longitudeₒ(WIII) == 0.0u"°"
+    @test altitudeₒ(WIII) == 0.0u"m"
+  end
 end