Skip to content

Commit

Permalink
Add Orthographic inverse (#9)
Browse files Browse the repository at this point in the history
* Add 'Orthographic' inverse

* Fix code
  • Loading branch information
eliascarv authored Feb 28, 2024
1 parent 160f3c8 commit accb7a8
Show file tree
Hide file tree
Showing 2 changed files with 76 additions and 3 deletions.
49 changes: 49 additions & 0 deletions src/crs/projected/orthographic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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₀)))
Expand Down Expand Up @@ -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
30 changes: 27 additions & 3 deletions test/crs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

0 comments on commit accb7a8

Please sign in to comment.