From 2d06e027f78e40b29355d9a059990cce85c25104 Mon Sep 17 00:00:00 2001 From: Tim Reichelt Date: Fri, 18 Oct 2024 19:53:14 +0100 Subject: [PATCH] Add Haversine and general `spherical_distance` and `AbstractSphericalDistance` functionality (#588) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Add haversine distance function * introduce AbstractSphericalDistance * changelog updated * functor method for type not instance * define DEFAULT_RADIUS as integer * Add tests * Add additional tests * reexport spherical distance * add type stability test * test spherical distance -180:180 or 0:360 * Docs for spherical distance * restrict rand to -180:360 integers * IntN.(::Tuple) for spherical distance test * (()) typo * clarification in docs --------- Co-authored-by: Milan Klöwer --- CHANGELOG.md | 1 + docs/src/orography.md | 41 ++++++++++++++++++ src/RingGrids/RingGrids.jl | 1 + src/RingGrids/geodesics.jl | 56 ++++++++++++++++++++++++ src/SpeedyWeather.jl | 1 + test/geodesics.jl | 88 ++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + 7 files changed, 189 insertions(+) create mode 100644 src/RingGrids/geodesics.jl create mode 100644 test/geodesics.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index f18b8c5b9..ce58efd02 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Unreleased +- 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) diff --git a/docs/src/orography.md b/docs/src/orography.md index 68b63ab12..c3341c226 100644 --- a/docs/src/orography.md +++ b/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`. diff --git a/src/RingGrids/RingGrids.jl b/src/RingGrids/RingGrids.jl index 2f90a70f5..b1ccb3046 100644 --- a/src/RingGrids/RingGrids.jl +++ b/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") diff --git a/src/RingGrids/geodesics.jl b/src/RingGrids/geodesics.jl new file mode 100644 index 000000000..be51b95f3 --- /dev/null +++ b/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...) + + diff --git a/src/SpeedyWeather.jl b/src/SpeedyWeather.jl index 3765758a8..b709b3742 100644 --- a/src/SpeedyWeather.jl +++ b/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 diff --git a/test/geodesics.jl b/test/geodesics.jl new file mode 100644 index 000000000..78d4db6f0 --- /dev/null +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 556bf37a0..ae0bac7a0 100644 --- a/test/runtests.jl +++ b/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")