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

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

radius = 1
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))
spherical = [Spherical(radius, deg2rad(longitude), deg2rad(latitude)) for (longitude, latitude) in lon_lat]
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)


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()
    rad2deg(dist(p1, p2))
hist(d; bins=100, axis=(xlabel="Angular distance (°)",))

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

γ = fit(Variogram, g)


This resulted in a SineHoleVariogram, which doesn’t really look promising.


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[1]), latitude in range(-π/2, π/2, dims[2])];
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")


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…


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[1])
viz(fig[1,2], solution[2])
viz(fig[1,3], solution[3])


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.


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 => (;))
    └─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 => (;))
    └─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.
 [1] 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
 [2] 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
 [3] top-level scope
   @ REPL[280]: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.


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
(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?


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

1 Like