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.
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")
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.
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")
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…