From accb7a8877020f7b20f962105afeeac07bde1ff3 Mon Sep 17 00:00:00 2001
From: Elias Carvalho <73039601+eliascarv@users.noreply.github.com>
Date: Wed, 28 Feb 2024 13:59:46 -0300
Subject: [PATCH] Add `Orthographic` inverse (#9)

* Add 'Orthographic' inverse

* Fix code
---
 src/crs/projected/orthographic.jl | 49 +++++++++++++++++++++++++++++++
 test/crs.jl                       | 30 +++++++++++++++++--
 2 files changed, 76 insertions(+), 3 deletions(-)

diff --git a/src/crs/projected/orthographic.jl b/src/crs/projected/orthographic.jl
index bce7763..9bb38f8 100644
--- a/src/crs/projected/orthographic.jl
+++ b/src/crs/projected/orthographic.jl
@@ -65,6 +65,15 @@ const OrthoSouth{Datum} = Orthographic{-90.0u"°",0.0u"°",false,Datum}
 # 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.
+# reference code: https://github.com/OSGeo/PROJ/blob/master/src/projections/robin.cpp
+# reference formulas: https://neacsu.net/docs/geodesy/snyder/5-azimuthal/sect_20/
+#                     https://mathworld.wolfram.com/OrthographicProjection.html
+#                     https://epsg.org/guidance-notes.html
+
 function formulas(::Type{<:Orthographic{lat₀,lon₀,false,Datum}}, ::Type{T}) where {lat₀,lon₀,Datum,T}
   λ₀ = T(ustrip(deg2rad(lon₀)))
   ϕ₀ = T(ustrip(deg2rad(lat₀)))
@@ -97,3 +106,43 @@ function formulas(::Type{<:Orthographic{lat₀,lon₀,true,Datum}}, ::Type{T}) w
 
   fx, fy
 end
+
+function sphericalinv(x, y, λ₀, ϕ₀)
+  ρ = sqrt(x^2 + y^2)
+  c = asin(ρ)
+  if ρ < atol(x)
+    λ₀, ϕ₀
+  else
+    sinc = sin(c)
+    cosc = cos(c)
+    sinϕ₀ = sin(ϕ₀)
+    cosϕ₀ = cos(ϕ₀)
+    λ = λ₀ + atan(x * sinc / (ρ * cosϕ₀ * cosc - y * sinϕ₀ * sinc))
+    ϕ = asin(cosc * sinϕ₀ + (y * sinc * cosϕ₀ / ρ))
+    λ, ϕ
+  end
+end
+
+function Base.convert(::Type{LatLon{Datum}}, coords::C) where {lat₀,lon₀,Datum,C<:Orthographic{lat₀,lon₀,true,Datum}}
+  T = numtype(coords.x)
+  a = numconvert(T, majoraxis(ellipsoid(Datum)))
+  x = coords.x / a
+  y = coords.y / a
+  λ₀ = T(ustrip(deg2rad(lon₀)))
+  ϕ₀ = T(ustrip(deg2rad(lat₀)))
+  λ, ϕ = sphericalinv(x, y, λ₀, ϕ₀)
+  LatLon{Datum}(rad2deg(ϕ) * u"°", rad2deg(λ) * u"°")
+end
+
+function Base.convert(::Type{LatLon{Datum}}, coords::C) where {lat₀,lon₀,S,Datum,C<:Orthographic{lat₀,lon₀,S,Datum}}
+  T = numtype(coords.x)
+  a = numconvert(T, majoraxis(ellipsoid(Datum)))
+  x = coords.x / a
+  y = coords.y / a
+  λ₀ = T(ustrip(deg2rad(lon₀)))
+  ϕ₀ = T(ustrip(deg2rad(lat₀)))
+  λₛ, ϕₛ = sphericalinv(x, y, λ₀, ϕ₀)
+  fx, fy = formulas(C, T)
+  λ, ϕ = projinv(fx, fy, x, y, λₛ, ϕₛ)
+  LatLon{Datum}(rad2deg(ϕ) * u"°", rad2deg(λ) * u"°")
+end
diff --git a/test/crs.jl b/test/crs.jl
index 10f5ce5..7d59b19 100644
--- a/test/crs.jl
+++ b/test/crs.jl
@@ -1006,54 +1006,78 @@
       c1 = LatLon(T(30), T(60))
       c2 = convert(OrthoNorth{WGS84}, c1)
       @test c2 ≈ OrthoNorth(T(4787610.688267582), T(-2764128.319646418))
+      c3 = convert(LatLon{WGS84}, c2)
+      @test c3 ≈ c1
 
       c1 = LatLon(T(30), -T(60))
       c2 = convert(OrthoNorth{WGS84}, c1)
       @test c2 ≈ OrthoNorth(-T(4787610.688267582), T(-2764128.319646418))
+      c3 = convert(LatLon{WGS84}, c2)
+      @test c3 ≈ c1
 
       # type stability
       c1 = LatLon(T(30), T(60))
+      c2 = OrthoNorth(T(4787610.688267582), T(-2764128.319646418))
       @inferred convert(OrthoNorth{WGS84}, c1)
+      @inferred convert(LatLon{WGS84}, c2)
     end
 
     @testset "LatLon <> OrthoSouth" begin
       c1 = LatLon(-T(30), T(60))
       c2 = convert(OrthoSouth{WGS84}, c1)
       @test c2 ≈ OrthoSouth(T(4787610.688267582), T(2764128.319646418))
+      c3 = convert(LatLon{WGS84}, c2)
+      @test c3 ≈ c1
 
       c1 = LatLon(-T(30), -T(60))
       c2 = convert(OrthoSouth{WGS84}, c1)
       @test c2 ≈ OrthoSouth(-T(4787610.688267582), T(2764128.319646418))
+      c3 = convert(LatLon{WGS84}, c2)
+      @test c3 ≈ c1
 
       # type stability
-      c1 = LatLon(T(30), T(60))
+      c1 = LatLon(-T(30), T(60))
+      c2 = OrthoSouth(T(4787610.688267582), T(2764128.319646418))
       @inferred convert(OrthoSouth{WGS84}, c1)
+      @inferred convert(LatLon{WGS84}, c2)
     end
 
     @testset "LatLon <> OrthoSpherical" begin
-      OrthoNorthSpherical = Cartography.Orthographic{90.0u"°",0.0u"°",true,WGS84}
-      OrthoSouthSpherical = Cartography.Orthographic{-90.0u"°",0.0u"°",true,WGS84}
+      OrthoNorthSpherical = Cartography.crs(ESRI{102035})
+      OrthoSouthSpherical = Cartography.crs(ESRI{102037})
 
       c1 = LatLon(T(30), T(60))
       c2 = convert(OrthoNorthSpherical, c1)
       @test c2 ≈ OrthoNorthSpherical(T(4783602.75), T(-2761814.335408735))
+      c3 = convert(LatLon{WGS84}, c2)
+      @test c3 ≈ c1
 
       c1 = LatLon(T(30), -T(60))
       c2 = convert(OrthoNorthSpherical, c1)
       @test c2 ≈ OrthoNorthSpherical(-T(4783602.75), T(-2761814.335408735))
+      c3 = convert(LatLon{WGS84}, c2)
+      @test c3 ≈ c1
 
       c1 = LatLon(-T(30), T(60))
       c2 = convert(OrthoSouthSpherical, c1)
       @test c2 ≈ OrthoSouthSpherical(T(4783602.75), T(2761814.335408735))
+      c3 = convert(LatLon{WGS84}, c2)
+      @test c3 ≈ c1
 
       c1 = LatLon(-T(30), -T(60))
       c2 = convert(OrthoSouthSpherical, c1)
       @test c2 ≈ OrthoSouthSpherical(-T(4783602.75), T(2761814.335408735))
+      c3 = convert(LatLon{WGS84}, c2)
+      @test c3 ≈ c1
 
       # type stability
       c1 = LatLon(T(30), T(60))
+      c2 = OrthoNorthSpherical(T(4783602.75), T(-2761814.335408735))
+      c3 = OrthoSouthSpherical(T(4783602.75), T(2761814.335408735))
       @inferred convert(OrthoNorthSpherical, c1)
       @inferred convert(OrthoSouthSpherical, c1)
+      @inferred convert(LatLon{WGS84}, c2)
+      @inferred convert(LatLon{WGS84}, c3)
     end
   end
 end