Add Orthographic projection (#762)

* Add 'Orthographic' projection

* Apply suggestions

* Update code

src/Meshes.jl

@@ -164,6 +164,8 @@ export
   GallPeters,
   WinkelTripel,
   Robinson,
+  OrthoNorth,
+  OrthoSouth,
   datum,
 
   # geometries

src/crs.jl

@@ -146,6 +146,7 @@ include("crs/eqdistcylindrical.jl")
 include("crs/eqareacylindrical.jl")
 include("crs/winkeltripel.jl")
 include("crs/robinson.jl")
+include("crs/orthographic.jl")
 
 # ----------
 # FALLBACKS

src/crs/orthographic.jl

@@ -0,0 +1,97 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    Orthographic{lat₀,lon₀,S}
+
+Orthographic CRS with latitude origin `lat₀` and longitude origin `lon₀` in degrees
+and spherical mode `S` enabled or not.
+"""
+const Orthographic{lat₀,lon₀,S,M<:Met} = CRS{:Orthographic,@NamedTuple{x::M, y::M},WGS84,Tuple{lat₀,lon₀,S}}
+
+Orthographic{lat₀,lon₀,S}(x::M, y::M) where {lat₀,lon₀,S,M<:Met} = Orthographic{lat₀,lon₀,S,float(M)}(x, y)
+Orthographic{lat₀,lon₀,S}(x::Met, y::Met) where {lat₀,lon₀,S} = Orthographic{lat₀,lon₀,S}(promote(x, y)...)
+Orthographic{lat₀,lon₀,S}(x::Len, y::Len) where {lat₀,lon₀,S} =
+  Orthographic{lat₀,lon₀,S}(uconvert(u"m", x), uconvert(u"m", y))
+Orthographic{lat₀,lon₀,S}(x::Number, y::Number) where {lat₀,lon₀,S} =
+  Orthographic{lat₀,lon₀,S}(addunit(x, u"m"), addunit(y, u"m"))
+
+"""
+    OrthoNorth(x, y)
+
+Orthographic North Pole coordinates in length units (default to meter).
+
+## Examples
+
+```julia
+OrthoNorth(1, 1) # add default units
+OrthoNorth(1u"m", 1u"m") # integers are converted converted to floats
+OrthoNorth(1.0u"km", 1.0u"km") # length quantities are converted to meters
+OrthoNorth(1.0u"m", 1.0u"m")
+```
+"""
+const OrthoNorth = Orthographic{90.0u"°",0.0u"°",false}
+
+"""
+    OrthoSouth(x, y)
+
+Orthographic South Pole coordinates in length units (default to meter).
+
+## Examples
+
+```julia
+OrthoSouth(1, 1) # add default units
+OrthoSouth(1u"m", 1u"m") # integers are converted converted to floats
+OrthoSouth(1.0u"km", 1.0u"km") # length quantities are converted to meters
+OrthoSouth(1.0u"m", 1.0u"m")
+```
+"""
+const OrthoSouth = Orthographic{-90.0u"°",0.0u"°",false}
+
+typealias(::Type{ESRI{102035}}) = Orthographic{90.0u"°",0.0u"°",true}
+
+typealias(::Type{ESRI{102037}}) = Orthographic{-90.0u"°",0.0u"°",true}
+
+# ------------
+# CONVERSIONS
+# ------------
+
+function Base.convert(::Type{Orthographic{lat₀,lon₀,false}}, coords::LatLon) where {lat₀,lon₀}
+  🌎 = ellipsoid(Orthographic)
+  λ = deg2rad(coords.lon)
+  ϕ = deg2rad(coords.lat)
+  λ₀ = oftype(λ, deg2rad(lon₀))
+  ϕ₀ = oftype(ϕ, deg2rad(lat₀))
+  l = ustrip(λ)
+  a = oftype(l, ustrip(majoraxis(🌎)))
+  e² = oftype(l, eccentricity²(🌎))
+
+  sinϕ = sin(ϕ)
+  cosϕ = cos(ϕ)
+  sinϕ₀ = sin(ϕ₀)
+  cosϕ₀ = cos(ϕ₀)
+  ν = 1 / sqrt(1 - e² * sinϕ^2)
+  ν₀ = 1 / sqrt(1 - e² * sinϕ₀^2)
+
+  x = a * (ν * cosϕ * sin(λ - λ₀))
+  y = a * (ν * (sinϕ * cosϕ₀ - cosϕ * sinϕ₀ * cos(λ - λ₀)) + e² * (ν₀ * sinϕ₀ - ν * sinϕ) * cosϕ₀)
+
+  Orthographic{lat₀,lon₀,false}(x * u"m", y * u"m")
+end
+
+function Base.convert(::Type{Orthographic{lat₀,lon₀,true}}, coords::LatLon) where {lat₀,lon₀}
+  🌎 = ellipsoid(Orthographic)
+  λ = deg2rad(coords.lon)
+  ϕ = deg2rad(coords.lat)
+  λ₀ = oftype(λ, deg2rad(lon₀))
+  ϕ₀ = oftype(ϕ, deg2rad(lat₀))
+  l = ustrip(λ)
+  a = oftype(l, ustrip(majoraxis(🌎)))
+
+  cosϕ = cos(ϕ)
+  x = a * cosϕ * sin(λ - λ₀)
+  y = a * (sin(ϕ) * cos(ϕ₀) - cosϕ * sin(ϕ₀) * cos(λ - λ₀))
+
+  Orthographic{lat₀,lon₀,true}(x * u"m", y * u"m")
+end

test/crs.jl

@@ -691,5 +691,68 @@
       @inferred convert(Robinson, c1)
       @inferred convert(ESRI{54030}, c1)
     end
+
+    @testset "LatLon <> OrthoNorth" begin
+      c1 = LatLon(T(30), T(60))
+      c2 = convert(OrthoNorth, c1)
+      @test c2 ≈ OrthoNorth(T(4787610.688267582), T(-2764128.319646418))
+
+      c1 = LatLon(T(30), -T(60))
+      c2 = convert(OrthoNorth, c1)
+      @test c2 ≈ OrthoNorth(-T(4787610.688267582), T(-2764128.319646418))
+
+      # type stability
+      c1 = LatLon(T(30), T(60))
+      @inferred convert(OrthoNorth, c1)
+    end
+
+    @testset "LatLon <> OrthoSouth" begin
+      c1 = LatLon(-T(30), T(60))
+      c2 = convert(OrthoSouth, c1)
+      @test c2 ≈ OrthoSouth(T(4787610.688267582), T(2764128.319646418))
+
+      c1 = LatLon(-T(30), -T(60))
+      c2 = convert(OrthoSouth, c1)
+      @test c2 ≈ OrthoSouth(-T(4787610.688267582), T(2764128.319646418))
+
+      # type stability
+      c1 = LatLon(T(30), T(60))
+      @inferred convert(OrthoSouth, c1)
+    end
+
+    @testset "LatLon <> OrthoSpherical" begin
+      OrthoNorthSpherical = Meshes.typealias(ESRI{102035})
+      OrthoSouthSpherical = Meshes.typealias(ESRI{102037})
+
+      c1 = LatLon(T(30), T(60))
+      c2 = convert(OrthoNorthSpherical, c1)
+      @test c2 ≈ OrthoNorthSpherical(T(4783602.75), T(-2761814.335408735))
+
+      c1 = LatLon(T(30), -T(60))
+      c2 = convert(OrthoNorthSpherical, c1)
+      @test c2 ≈ OrthoNorthSpherical(-T(4783602.75), T(-2761814.335408735))
+
+      c1 = LatLon(-T(30), T(60))
+      c2 = convert(OrthoSouthSpherical, c1)
+      @test c2 ≈ OrthoSouthSpherical(T(4783602.75), T(2761814.335408735))
+
+      c1 = LatLon(-T(30), -T(60))
+      c2 = convert(OrthoSouthSpherical, c1)
+      @test c2 ≈ OrthoSouthSpherical(-T(4783602.75), T(2761814.335408735))
+
+      # EPSG/ESRI fallback
+      c1 = LatLon(T(30), T(60))
+      c2 = convert(ESRI{102035}, c1)
+      @test c2 ≈ OrthoNorthSpherical(T(4783602.75), T(-2761814.335408735))
+      c2 = convert(ESRI{102037}, c1)
+      @test c2 ≈ OrthoSouthSpherical(T(4783602.75), T(2761814.335408735))
+
+      # type stability
+      c1 = LatLon(T(30), T(60))
+      @inferred convert(OrthoNorthSpherical, c1)
+      @inferred convert(ESRI{102035}, c1)
+      @inferred convert(OrthoSouthSpherical, c1)
+      @inferred convert(ESRI{102037}, c1)
+    end
   end
 end