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
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)
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()
rad2deg(dist(p1, p2))
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[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")
Artefacts
Two main artefacts are visible:
- There is a “fault line” where the longitude angles start and end
- 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…
But the interpolation you shared above seems reasonable to me, at least visually.