SPDE simulation on a sphere¶

In [1]:
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.document as gdoc
import gstlearn.plot3D as gop
import numpy as np 
import matplotlib.pyplot as plt
import plotly.graph_objects as go
import IPython
from IPython.display import Markdown
gdoc.setNoScroll()
In [2]:
Markdown(gdoc.loadDoc("SPDE_Sphere.md"))
Out[2]:

SPDE on the Sphere

The aim of this tutorial is to show how to use gstlearn to simulate the solution of

$$(\kappa^2-\Delta_{\mathcal{S}_R})^{\alpha/2}Z = \sigma\mathcal{W}$$

on the sphere $\mathcal{S}_R$ of radius $R$.

  • $\Delta_{{\mathcal{S}_R}}$ is the Laplace-Beltrami operator, i.e, it acts on each point of the sphere as the usual Laplacian on the tangent plane at this point.

  • $\kappa$ is the inverse of the scale parameter

  • $\alpha \geq 2$ is an integer describing the smoothness of the solution.

  • $\mathcal{W}$ is a Gaussian white-noise suitably normalized such as $\sigma^2$ is the variance of the solution.

In this notebook, we will define the covariance of Matérn on the sphere, as the covariance of the solution of this SPDE (other extensions of the Matérn function are possible). By analogy with the Euclidian case, its smoothness parameter will be defined by $\nu = \alpha -1$. To compute the covariance function with respect on the geodetic distance, one have to use a decomposition on the Legendre polynomial (see below).

We also treat the more general case $$P^{1/2}(-\Delta_{\mathcal{S}_R})Z = \sigma\mathcal{W}$$

where $P$ is a polynom positive on $\mathbb{R}^+$
In [3]:
# Markov (if False, Matérn covariance will be used)
Markov = True

Parametrization¶

In [4]:
#Sphere radius
R = gl.EARTH_RADIUS
gl.defineDefaultSpace(gl.ESpaceType.SN,param=R)

#Scale parameter (for convenience, it is defined as the proportion of the radius)
ratioRange = 0.2
scale = R * ratioRange

# sill 
sill = 2. 

# Smoothness parameter (for Matérn case)
nu = 2

# Markov coefficients (for Markov case)
coeffs = [1,-1,.5]

Meshing¶

In [5]:
mesh = gl.MeshSphericalExt()
err = mesh.resetFromDb(None,None,triswitch = "-r5",verbose=False)

Sampling Db creation

In [6]:
nsample = 4000
#sub-sampling to reduce computational burden
np.random.seed(123)
ind = np.random.choice(mesh.getNApices(),size=nsample,replace=False)

#Creation of the db
X = mesh.getCoordinates(0)
Y = mesh.getCoordinates(1)

dbdat = gl.Db.create()
dbdat["x"] = np.array(X)[ind]
dbdat["y"] = np.array(Y)[ind]
dbdat.setLocators(["x","y"],gl.ELoc.X)

Covariance model¶

The covariance model is Markov or Matérn.

In [7]:
if Markov : 
    model = gl.Model.createFromParam(type=gl.ECov.MARKOV,
                                 range = scale,
                                 sill = sill,
                                 flagRange= False)
    model.setMarkovCoeffs(0, coeffs)
    
else :

    model = gl.Model.createFromParam(type=gl.ECov.MATERN,
                                 range = scale,
                                 sill = sill,
                                 param=nu,
                                 flagRange= False)

Precision matrix¶

In [8]:
Q = gl.PrecisionOp(mesh,model.getCova(0))

Simulation¶

In [9]:
result = np.array(Q.simulateOne())

Display the realization¶

In [10]:
surface = gop.SurfaceOnMesh(mesh, result,opacity=1)
fig = go.Figure(data=[ surface ])
fig.update_scenes(xaxis_visible=False, yaxis_visible=False,zaxis_visible=False )
f = fig.show()