Add 'WinkelTripel' projection (#752)

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

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

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ₜₛ))

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(λ)

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(λ)

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

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"

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

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

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