Merge branch 'mk/haurwitz' of https://github.com/SpeedyWeather/Speedy…

…Weather.jl into mk/haurwitz

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Unreleased
 
 - Rossby-Haurwitz wave initial conditions [#591](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/591)
+- Haversine formula and AbstractSphericalDistance [#588](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/588)
 - Array-agnostic SpectralTransform [#583](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/583)
 - Move CUDA dependency into extension [#586](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/586)
 - Stop supporting Julia v1.9 [#585](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/585)

docs/src/orography.md

@@ -171,6 +171,47 @@ as illustrated in the spectral representation or the surface geopotential
 here. This is because scalar fields do not use that last degree,
 see [One more degree for spectral fields](@ref).
 
+## Spherical distance
+
+In the example above we have defined the "Super Hawaii orography" as
+
+```julia
+orography=(λ,φ) -> H*exp((-(λ-λ₀)^2 - (φ-φ₀)^2)/2σ^2)
+```
+
+note however, that this approximates the distance on the sphere with
+Cartesian coordinates which is here not too bad as we are not too far north
+(where longitudinal distances would become considerably shorter) and also as
+we are far away from the prime meridian. If ``\lambda_0 = 0`` in the example
+above then calculating ``\lambda - \lambda_0`` for ``\lambda = 359˚E``
+yields a really far distance even though ``\lambda`` is actually relatively close to
+the prime meridian. To avoid this problem, SpeedyWeather (actually RingGrids)
+defines a function called `spherical_distance` (inputs in degrees, output in meters) to actually
+calculate the [great-circle distance](https://en.wikipedia.org/wiki/Great-circle_distance)
+or spherical distance. Compare
+
+```@example orography
+λ₀ = 0         # move "Super Hawaii" mountain onto the prime meridian
+set!(model, orography=(λ,φ) -> H*exp((-(λ-λ₀)^2 - (φ-φ₀)^2)/2σ^2))
+heatmap(model.orography.orography, title="Mountain on prime meridian, cartesian coordinates")
+save("mountain_cartesian.png", ans) # hide
+nothing # hide
+```
+![Mountain in cartesian coordinates](mountain_cartesian.png)
+
+which clearly shows that the mountain is only on the Eastern hemisphere
+(the signal on the western hemisphere is due to the interpolation during plotting)
+-- probably not what you wanted. Rewrite this as
+
+```@example orography
+λ₀ = 0         # move "Super Hawaii" mountain onto the prime meridian
+set!(model, orography=(λ,φ) -> H*exp(-spherical_distance((λ,φ), (λ₀,φ₀), radius=1)^2/2σ^2))
+heatmap(model.orography.orography, title="Mountain on prime meridian, spherical distance")
+save("mountain_spherical.png", ans) # hide
+nothing # hide
+```
+![Mountain in spherical coordinates](mountain_spherical.png)
+
 ## Defining a new orography type
 
 You can also define a new orography like we defined `ZonalRidge` or `EarthOrography`.

src/RingGrids/RingGrids.jl

@@ -100,6 +100,7 @@ include("general.jl")
 include("full_grids.jl")
 include("reduced_grids.jl")
 include("scaling.jl")
+include("geodesics.jl")
 
 # FULL GRIDS
 include("grids/full_gaussian.jl")

src/RingGrids/geodesics.jl

@@ -0,0 +1,56 @@
+const DEFAULT_RADIUS = 6371000  # Earth mean radius in meters, defined as integer for more type stability 
+
+"""
+    abstract type AbstractSphericalDistance end
+
+Super type of formulas to calculate the spherical distance or great-circle distance.
+To define a `NewFormula`, define `struct NewFormula <: AbstractSphericalDistance end`
+and the actual calculation as a functor
+
+    function NewFormula(lonlat1::Tuple, lonlat2::Tuple; radius=DEFAULT_RADIUS, kwargs...)
+
+assuming inputs in degrees and returning the distance in meters (or radians for radius=1).
+Then use the general interface `spherical_distance(NewFormula, args...; kwargs...)`"""
+abstract type AbstractSphericalDistance end
+struct Haversine <: AbstractSphericalDistance end
+
+"""$(TYPEDSIGNATURES) 
+Haversine formula calculating the great-circle or spherical distance (in meters) on the sphere
+between two tuples of  longitude-latitude points in degrees ˚E, ˚N. Use keyword argument `radius`
+to change the radius of the sphere (default 6371e3 meters, Earth's radius), use `radius=1`
+to return the central angle in radians."""
+function Haversine(             # functor for Haversine struct
+    lonlat1::Tuple,             # point 1 in spherical coordinates
+    lonlat2::Tuple;             # point 2
+    radius = DEFAULT_RADIUS,    # radius of the sphere [m]
+)
+    lon1, lat1 = lonlat1
+    lon2, lat2 = lonlat2    
+
+    φ1 = deg2rad(lat1)
+    φ2 = deg2rad(lat2)
+
+    Δφ = deg2rad(lat2 - lat1)
+    Δλ = deg2rad(lon2 - lon1)
+
+    # Haversine formula, see https://en.wikipedia.org/wiki/Haversine_formula
+    a = sin(Δφ / 2)^2 + cos(φ1) * cos(φ2) * sin(Δλ / 2)^2
+    c = 2 * atan(sqrt(a), sqrt(1 - a))
+    return c*radius             # in meters
+end
+
+# allow for any <:AbstractSphericalDistance also non-tupled arguments lon1, lat1, lon2, lat2
+(F::Type{<:AbstractSphericalDistance})(lon1, lat1, lon2, lat2; kwargs...) = F((lon1, lat1), (lon2, lat2); kwargs...) 
+
+export spherical_distance
+
+"""$(TYPEDSIGNATURES)
+Spherical distance, or great-circle distance, between two points `lonlat1` and `lonlat2`
+using the `Formula` (default `Haversine`)."""
+spherical_distance(Formula::Type{<:AbstractSphericalDistance}, args...; kwargs...) =
+    Formula(args...; kwargs...)
+
+# define Haversine as default
+spherical_distance(args...; kwargs...) = spherical_distance(Haversine, args...; kwargs...)
+
+

src/SpeedyWeather.jl

@@ -60,6 +60,7 @@ export  FullClenshawGrid, FullClenshawArray,
         OctaHEALPixGrid, OctaHEALPixArray,
         eachring, eachgrid, plot
 export  AnvilInterpolator
+export  spherical_distance
 
 include("RingGrids/RingGrids.jl")
 using .RingGrids

test/geodesics.jl

@@ -0,0 +1,88 @@
+@testset "Spherical distance on Earth" begin
+    earth_circumference = 2π*RingGrids.DEFAULT_RADIUS
+    for (point1, point2, expected) in [
+
+        # longitude 
+        ((0, 0), (1, 0), earth_circumference / 360),
+        ((0, 0), (359, 0), earth_circumference / 360),
+        ((0, 0), (180, 0), earth_circumference / 2),
+        ((0, 0), (90, 0), earth_circumference / 4),
+
+        # latitude
+        ((0, 0), (0, 1), earth_circumference / 360),
+        ((0, 0), (0, 10), 10*earth_circumference / 360),
+        ((0, 0), (0, 90), earth_circumference / 4),
+        ((0, 0), (0, -90), earth_circumference / 4),
+    ]
+        @test spherical_distance(point1, point2) ≈ expected
+    end
+end
+
+@testset "Spherical distance degrees -180:180 or 0:360" begin
+    # either -180 to 180 or 0 to 360˚E
+    @test spherical_distance((-90, 0), (-180, 0)) ≈ spherical_distance((270, 0), (180, 0))
+end
+
+@testset "Spherical distance invariants" begin
+
+    # Ensure that function preserves the invariant that shifting in the longitude
+    # direction does not change the distance.
+    rand_lonlat() = (360 * rand(), 180 * rand() - 90)
+    shift((lon, lat), lon_shift) = ((lon + lon_shift) % 360, lat)
+    for _ in 1:100
+        p1, p2 = rand_lonlat(), rand_lonlat()
+        lon_shift = 360 * rand() 
+
+        # Invariance in longitude shift.
+        @test spherical_distance(p1, p2) ≈ spherical_distance(shift(p1, lon_shift), shift(p2, lon_shift))
+        # Commutative.
+        @test spherical_distance(p1, p2) ≈ spherical_distance(p2, p1)
+        # Symmetry.
+        @test spherical_distance(p1, p2) == spherical_distance((-p1[1], -p1[2]), (-p2[1], -p2[2]))
+        # Identity.
+        @test spherical_distance(p1, p1) == 0
+
+        # Ensure that different ways to call the function are identical.
+        @test spherical_distance(p1, p2) == spherical_distance(p1[1], p1[2], p2[1], p2[2])
+        @test spherical_distance(p1, p2) == spherical_distance(p1, p2)
+        @test spherical_distance(p1, p2) == spherical_distance(RingGrids.Haversine, p1, p2)
+    end
+end
+
+@testset "Spherical distance with custom radius" begin
+    custom_radius = 100
+    circumference = 2 * π * custom_radius
+    for (point1, point2, expected) in [
+        ((0, 0), (1, 0), circumference / 360),
+        ((0, 0), (359, 0), circumference / 360),
+        ((0, 0), (180, 0), circumference / 2),
+        ((0, 0), (90, 0), circumference / 4),
+
+        ((0, 0), (0, 1), circumference / 360),
+        ((0, 0), (0, 10), 10*circumference / 360),
+        ((0, 0), (0, 90), circumference / 4),
+        ((0, 0), (0, -90), circumference / 4),
+    ]
+        # There's a small numerical error, hence only approximate equality.
+        @test spherical_distance(point1, point2, radius=custom_radius) ≈ expected
+    end
+end
+
+@testset "Spherical distance type stability" begin
+    for T in (Float16, Float32, Float64)
+        p1, p2 = (rand(T), rand(T)), (rand(T), rand(T))
+        @test typeof(spherical_distance(p1, p2)) == T
+
+        # but if radius provided promote to that type
+        for Tradius in (Float16, Float32, Float64)
+            @test typeof(spherical_distance(p1, p2, radius=Tradius(1))) == promote_type(T, Tradius)
+        end
+    end
+
+    # but integers always go to Float64 due to the trigonometric functions
+    for T in (Int16, Int32, Int64)
+        d = -180:360    # sensible range for degrees
+        p1, p2 = T.((rand(d), rand(d))), T.((rand(d), rand(d)))
+        @test typeof(spherical_distance(p1, p2)) == Float64
+    end
+end

test/runtests.jl

@@ -6,6 +6,7 @@ include("utility_functions.jl")
 include("dates.jl")
 include("lower_triangular_matrix.jl")
 include("grids.jl")
+include("geodesics.jl")
 include("interpolation.jl")
 include("set.jl")