# Kriging on surface of sphere

**URL:** <https://discourse.julialang.org/t/kriging-on-surface-of-sphere/85058>\
**Category:** Geo\
**Tags:** interpolations\
**Created:** [July 31, 2022, 1:28pm UTC](https://discourse.julialang.org/t/kriging-on-surface-of-sphere/85058 "2022-07-31T13:28:18Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [July 31, 2022, 1:28pm UTC](https://discourse.julialang.org/t/kriging-on-surface-of-sphere/85058/1 "2022-07-31T13:28:18Z")

</div>

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

```julia-auto
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:

```julia-auto
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")

```

 ![tmp](https://global.discourse-cdn.com/julialang/original/3X/6/5/658fc2f22fe5d9963c62354f1b816e4360cf442e.png)

# Fitting a Variogram

We now look at the empirical variogram of the data:

```julia-auto
using GeoStats, Plots
𝒟 = georef((; data), lon_lat)
g = EmpiricalVariogram(𝒟, :data, maxlag=π, distance=Haversine(radius), algo=:full)
Plots.plot(g)

```

 ![tmp](https://global.discourse-cdn.com/julialang/original/3X/4/5/458a1a880285b3348a8b3bcfde198657b12ed45a.png)

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:

```julia-auto
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 (°)",))

```

 ![tmp](https://global.discourse-cdn.com/julialang/original/3X/6/9/6972a039b6ce2840976c5fde3b815119b8d4499e.png)

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

```julia-auto
γ = fit(Variogram, g)
Plots.plot(g)
Plots.plot!(γ)

```

 ![tmp](https://global.discourse-cdn.com/julialang/original/3X/6/2/6214406cbe50609fc451bb3d2eb80c85b1ef03f2.png)

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:

```julia-auto
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:

```julia-auto
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")

```

 ![tmp](https://global.discourse-cdn.com/julialang/original/3X/5/d/5d457c94381d1e13b9eb46fad09f0c5153382c9f.png)

# 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…

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [July 31, 2022, 2:14pm UTC](https://discourse.julialang.org/t/kriging-on-surface-of-sphere/85058/2 "2022-07-31T14:14:16Z")

</div>

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:

```julia
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)

```

```julia
3D Ensemble
  domain: 20000 SimpleMesh{3,Float64}
  variables: Z
  N° reals: 3

```

and here is the visualization:

```julia
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])

Mke.current_figure()

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/6/1635811b219325838582c376278e12f0ea0faadc.jpeg)

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:

```julia
L = laplacematrix(mesh)

```

```julia
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.

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [August 1, 2022, 8:32am UTC](https://discourse.julialang.org/t/kriging-on-surface-of-sphere/85058/3 "2022-08-01T08:32:00Z")

</div>

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:

```julia
# 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:

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

```

I then tried your sphere domain:

```julia
𝒟 = 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
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:
 [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

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

```

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [August 1, 2022, 11:04am UTC](https://discourse.julialang.org/t/kriging-on-surface-of-sphere/85058/4 "2022-08-01T11:04:32Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [August 1, 2022, 12:41pm UTC](https://discourse.julialang.org/t/kriging-on-surface-of-sphere/85058/5 "2022-08-01T12:41:06Z")

</div>

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

```julia
# 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:

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

```

and used that in order to solve the problem:

```julia
# 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)

```

 ![tmp](https://global.discourse-cdn.com/julialang/original/3X/e/6/e63a5ef8e0e8477e39e4c994e1ccc3e2f9ca0445.png)

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

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [August 1, 2022, 1:47pm UTC](https://discourse.julialang.org/t/kriging-on-surface-of-sphere/85058/6 "2022-08-01T13:47:39Z")

</div>

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