# Kriging on surface of sphere

I have a few random points on a sphere’s surface and I’d like to interpolate and extrapolate those using kriging. But I’m getting artefacts in the result that look like an issue with the fact that the coordinates are spherical not cartesian.

# The data

``````using Random
Random.seed!(1234);

function rand_on_sphere() # in degrees
longitude = 360rand() .- 180
latitude = acosd.(1 .- 2rand()) .- 90
return longitude, latitude
end

n = 100
lon_lat = [rand_on_sphere() for _ in 1:n]
data = 2rand(n) .- 1 # some value between -1 and 1
``````

which looks like this:

``````using GLMakie, CoordinateTransformations
fig = Figure()
ax = Axis3(fig[1,1], aspect=:data)
mesh!(ax, Sphere(zero(Point3), radius); color = (:gray, 1))
cartesian = CartesianFromSpherical().(spherical)
h = scatter!(ax, cartesian; color=data, colormap=:bluesreds)
Colorbar(fig[1,2], h, label = "data")
``````

# Fitting a Variogram

We now look at the empirical variogram of the data:

``````using GeoStats, Plots
𝒟 = georef((; data), lon_lat)
g = EmpiricalVariogram(𝒟, :data, maxlag=π, distance=Haversine(radius), algo=:full)
Plots.plot(g)
`````` I used a Haversine distance measure since the coordinates are on a sphere.

Note that the distribution of angular distances between randomly spaced points on a surface of a sphere has a distinct parabolic shape:

``````dist = Haversine(radius)
d = map(1:10^6) do _
p1 = rand_on_sphere()
p2 = rand_on_sphere()
end;
hist(d; bins=100, axis=(xlabel="Angular distance (°)",))
``````

We try to auto-fit a Variogram to the data:

``````γ = fit(Variogram, g)
Plots.plot(g)
Plots.plot!(γ)
`````` This resulted in a `SineHoleVariogram`, which doesn’t really look promising.

# Solution

Instead of a `SineHoleVariogram`, I try to solve this with a Gaussian Variogram:

``````dims = (100, 100)
𝒢 = CartesianGrid((-180, -90), (180, 90); dims)
problem = EstimationProblem(𝒟, 𝒢, :data)
solver = Kriging(:data => (variogram=GaussianVariogram(range=30), distance=Haversine(radius)))
solution = solve(problem, solver)
``````

but when I plot the solution, I get these weird artefacts:

``````xyz = [CartesianFromSpherical()(Spherical(radius, longitude, latitude)) for longitude in range(-π, π, dims), latitude in range(-π/2, π/2, dims)];
x = getindex.(xyz, 1);
y = getindex.(xyz, 2);
z = getindex.(xyz, 3);
v = reshape(values(solution).data, dims...);

fig = Figure()
ax = Axis3(fig[1,1], aspect=:data, elevation = 0.4π)
GLMakie.surface!(ax, x, y, z; color=v, shading=false, colormap=:bluesreds)
h = scatter!(ax, cartesian; color=data, colormap=:bluesreds)
Colorbar(fig[1,2], h, label = "data")
``````

# Artefacts

Two main artefacts are visible:

1. There is a “fault line” where the longitude angles start and end
2. The data seems to morph and squeeze into the poles

Both of these artefacts seem to have to do with the fact that the data is “circular”. Something I had hoped that the Haversine distance measure would taken care of that…

2 Likes

Hi @yakir12, on the sphere the methodology is slightly different as there is no “global” variogram model to fit before the interpolation. Here is a MWE where we simulate realizations on the sphere instead:

``````using GeoStats

sphere = Sphere((0.0,0.0,0.0), 1.0)

mesh = triangulate(sphere)

problem = SimulationProblem(mesh, :Z => Float64, 3)

solver = SPDEGS(:Z => (range=0.5, sill=1.0))

solution = solve(problem, solver)
``````
``````3D Ensemble
domain: 20000 SimpleMesh{3,Float64}
variables: Z
N° reals:  3
``````

and here is the visualization:

``````using MeshViz
import GLMakie as Mke

fig = Mke.Figure()

viz(fig[1,1], solution)
viz(fig[1,2], solution)
viz(fig[1,3], solution)

Mke.current_figure()
``````

In order to perform interpolation, we need to make use of the topology of the mesh, and build what is called the Laplace-Beltrami matrix:

``````L = laplacematrix(mesh)
``````
``````10002×10002 SparseArrays.SparseMatrixCSC{Float64, Int64} with 70002 stored entries:
⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢹
⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⢸
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⢸
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⢸
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⢸
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⢸
⣄⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣈⣻⣾
``````

I didn’t have time to wrap these low-level concepts into the high-level solvers yet, but could try to work on it over the week after I finish a deck of slides that is more urgent.

4 Likes

Thanks Júlio!
I managed to only partly adapt your simulation example to my use-case.

Consider the simpler case where the data is just 6 points, one at each of the intersection points of the sphere with the axes:

``````# longitude 0°..360° and latitude 0°..180°
lon_lat = [(0,0), (0,180), (0,90), (90, 90), (180, 90), (270,90)]
``````

and the values are zero at the equator and ±1 at the poles:

``````data = [1, -1, 0, 0, 0, 0]
``````

I then tried your sphere domain:

``````𝒟 = georef((; data), lon_lat)
sphere = GeoStats.Sphere((0,0,0), 1)
𝒢 = triangulate(sphere)
problem = EstimationProblem(𝒟, 𝒢, :data)
``````

So far so good (I hope).

But when I try to solve this using one of the `GeoEstimation` solvers that are not `Kriging`, i.e. `IDW` or `LWR` I run into problems:

``````julia> solver = IDW(:data => (;))
IDW
└─data
└─distance ⇨ Euclidean(0.0)
└─power ⇨ 1

julia> solution = solve(problem, solver)
20000 MeshData{3,Float64}
variables (rank 2)
└─data (Float64)
└─data_distance (Float64)
domain: 20000 SimpleMesh{3,Float64}

julia> all(isnan, values(solution).data)
true # so it's just an array of NaNs...

julia> solver = LWR(:data => (;))
LWR
└─data
└─distance ⇨ Euclidean(0.0)
└─weightfun ⇨ #20

julia> solution = solve(problem, solver)
ERROR: DimensionMismatch: first array has length 3 which does not match the length of the second, 4.
Stacktrace:
 dot(x::Vector{Float64}, y::StaticArraysCore.SVector{4, Float64})
@ LinearAlgebra ~/.julia/juliaup/julia-1.8.0-rc3+0.aarch64/share/julia/stdlib/v1.8/LinearAlgebra/src/generic.jl:876
 solve(problem::EstimationProblem{MeshData{PointSet{2, Float64}, Dict{Int64, TypedTables.Table{NamedTuple{(:data,), Tuple{Int64}}, 1, NamedTuple{(:data,), Tuple{Vector{Int64}}}}}}, SimpleMesh{3, Float64, Vector{Meshes.Point3}, FullTopology{Connectivity{Triangle{Dim, T} where {Dim, T}, 3}}}, 1}, solver::LWR)
@ GeoEstimation ~/.julia/packages/GeoEstimation/aVH3Q/src/lwr.jl:95
 top-level scope
@ REPL:1
``````

Any attempt at setting the distance to `Haversine` results in

``````ERROR: ArgumentError: expected both inputs to have length 2 in Haversine{Int64}(1) distance
``````

Hi @yakir12 , the geospatial data and domain must currently share the same coordinate system. I think you are mixing lat/lon coordinates for the measurements and x/y/z coordinates for the domain (sphere). Ultimately, we want to generalize this interface so that users can mix coordinate systems, but time is too limited and there are too few contributors.

2 Likes

OK, I switched from the spherical coordinate system to a cartesian one:

``````# the data in cartesian
data = Float64[1, -1, 0, 0, 0, 0]
xyz = NTuple{3, Float64}[(0,0,1), (0,0,-1), (0,1,0), (0,-1,0), (1,0,0), (-1,0,0)]

𝒟 = georef((; data), xyz)
sphere = GeoStats.Sphere((0.0,0.0,0.0), 1.0)
𝒢 = triangulate(sphere)
problem = EstimationProblem(𝒟, 𝒢, :data)
``````

But now I need to use a different distance measure, one that can calculate the angular difference between two cartesian coordinates on the surface of a sphere. To that end I defined a new distance measure:

``````using LinearAlgebra
import Distances.Metric
struct CartesianHaversine{T<:Number} <: Metric
end
(dist::CartesianHaversine)(x, y) = acos((x ⋅ y)/dist.radius)
``````

and used that in order to solve the problem:

``````# all three work
solver = LWR(:data => (distance=CartesianHaversine(1), ))
# solver = IDW(:data => (distance=CartesianHaversine(1), ))
# solver = Kriging(:data => (variogram=GaussianVariogram(range=30, ), distance=CartesianHaversine(1), ))
solution = solve(problem, solver)
``````

So I think problem solved…? Or would you disagree?

2 Likes

It depends on your ultimate goals But the interpolation you shared above seems reasonable to me, at least visually.

1 Like