Add backward methods for Horizontal Grid Shift and Point Motion (#16)

* Add backward methods for Horizontal Grid Shift and Point Motion

* Add utils.jl

* Add tests

* Update comments

* Fix code

* Update tests

* Fix code

* Add comments

* Rename variables

* Update comments

* More adjustments

src/CoordGridTransforms.jl

@@ -22,6 +22,7 @@ function __init__()
   ENV["DATADEPS_ALWAYS_ACCEPT"] = true
 end
 
+include("utils.jl")
 include("grids.jl")
 include("transforms.jl")
 

src/transforms/hgridshift.jl

@@ -19,6 +19,12 @@ macro hgridshift(Datumₛ, Datumₜ)
       latlon′ = hgridshiftfwd(Dₛ, Dₜ, latlon)
       LatLon{Dₜ}(latlon′...)
     end
+
+    function Base.convert(::Type{LatLon{Dₛ}}, coords::LatLon{Dₜ}) where {Dₛ<:$Datumₛ,Dₜ<:$Datumₜ}
+      latlon = (coords.lat, coords.lon)
+      latlon′ = hgridshiftbwd(Dₛ, Dₜ, latlon)
+      LatLon{Dₛ}(latlon′...)
+    end
   end
   esc(expr)
 end
@@ -28,6 +34,37 @@ function hgridshiftfwd(Datumₛ, Datumₜ, (lat, lon))
   lat + latshift, lon + lonshift
 end
 
+# 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/grids.cpp
+
+function hgridshiftbwd(Datumₛ, Datumₜ, (lat, lon))
+  tol = atol(numtype(lon))
+  # initial guess
+  latshift, lonshift = hgridshiftparams(Datumₛ, Datumₜ, lat, lon)
+  latᵢ = lat - latshift
+  lonᵢ = lon - lonshift
+  for _ in 1:MAXITER
+    # to check if the guess is equivalent to the original Datumₛ coordinates,
+    # forward the guess coordinates to compare them with the Datumₜ input coordinates
+    latᵢ′, lonᵢ′ = hgridshiftfwd(Datumₛ, Datumₜ, (latᵢ, lonᵢ))
+    # difference between forward coordinates and input coordinates for comparison
+    latΔ = latᵢ′ - lat
+    lonΔ = lonᵢ′ - lon
+    # if the difference is small, stop the iteration
+    if hypot(latΔ, lonΔ) ≤ tol
+      break
+    end
+    # otherwise, subtract the difference and continue
+    latᵢ -= latΔ
+    lonᵢ -= lonΔ
+  end
+  latᵢ, lonᵢ
+end
+
 function hgridshiftparams(Datumₛ, Datumₜ, lat, lon)
   T = numtype(lon)
   itp = interpolatelatlon(Datumₛ, Datumₜ, lat, lon)

src/transforms/pointmotion.jl

@@ -18,6 +18,12 @@ macro pointmotion(Datumₛ, Datumₜ)
       latlonalt′ = pointmotionfwd(Dₛ, Dₜ, latlonalt)
       LatLonAlt{Dₜ}(latlonalt′...)
     end
+
+    function Base.convert(::Type{LatLonAlt{Dₛ}}, coords::LatLonAlt{Dₜ}) where {Dₛ<:$Datumₛ,Dₜ<:$Datumₜ}
+      latlonalt = (coords.lat, coords.lon, coords.alt)
+      latlonalt′ = pointmotionbwd(Dₛ, Dₜ, latlonalt)
+      LatLonAlt{Dₛ}(latlonalt′...)
+    end
   end
   esc(expr)
 end
@@ -40,6 +46,55 @@ function pointmotionfwd(Datumₛ, Datumₜ, (lat, lon, alt))
   lat′, lon′, alt′
 end
 
+# 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/grids.cpp
+
+function pointmotionbwd(Datumₛ, Datumₜ, (lat, lon, alt))
+  λ = ustrip(deg2rad(lon))
+  ϕ = ustrip(deg2rad(lat))
+  h = ustrip(m, alt)
+  tol = atol(λ)
+
+  # initial guess
+  λₛ, ϕₛ, hₛ = pointmotionparams(Datumₛ, Datumₜ, lat, lon, ϕ, h)
+  λᵢ = λ - λₛ
+  ϕᵢ = ϕ - ϕₛ
+  hᵢ = h - hₛ
+  for _ in 1:MAXITER
+    # to check if the guess is equivalent to the original Datumₛ coordinates,
+    # forward the guess coordinates to compare them with the Datumₜ input coordinates
+    latᵢ = rad2deg(ϕᵢ) * °
+    lonᵢ = rad2deg(λᵢ) * °
+    λᵢₛ, ϕᵢₛ, hᵢₛ = pointmotionparams(Datumₛ, Datumₜ, latᵢ, lonᵢ, ϕᵢ, hᵢ)
+    λᵢ′ = λᵢ + λᵢₛ
+    ϕᵢ′ = ϕᵢ + ϕᵢₛ
+    hᵢ′ = hᵢ + hᵢₛ
+    # difference between forward coordinates and input coordinates for comparison
+    λΔ = λᵢ′ - λ
+    ϕΔ = ϕᵢ′ - ϕ
+    hΔ = hᵢ′ - h
+    # if the difference is small, stop the iteration
+    if hypot(λΔ, ϕΔ, hΔ) ≤ tol
+      break
+    end
+    # otherwise, subtract the difference and continue
+    λᵢ -= λΔ
+    ϕᵢ -= ϕΔ
+    hᵢ -= hΔ
+  end
+
+  # https://github.com/PainterQubits/Unitful.jl/issues/753
+  lonᵢ = rad2deg(λᵢ) * °
+  latᵢ = rad2deg(ϕᵢ) * °
+  altᵢ = uconvert(unit(alt), hᵢ * m)
+
+  latᵢ, lonᵢ, altᵢ
+end
+
 function pointmotionparams(Datumₛ, Datumₜ, lat, lon, ϕ, h)
   🌎 = ellipsoid(Datumₛ)
   T = numtype(lon)

src/utils.jl

@@ -0,0 +1,17 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+# maximum number of iterations in iterative algorithms
+const MAXITER = 10
+
+"""
+    atol(T)
+    atol(x::T)
+
+Absolute tolerance used in source code for approximate
+comparison with numbers of type `T`.
+"""
+atol(x) = atol(typeof(x))
+atol(::Type{Float64}) = 1e-10
+atol(::Type{Float32}) = 1.0f-5