In [1]:
import gstlearn as gl
import numpy as np
import matplotlib.pyplot as plt
In [2]:
cov = gl.CovDiffusionAdvection()

vel = [0.,0]
scaleTime = 0.2
sigma2 = 3
p1 = 1
p2 = 1
scale = 2
ctxt = gl.CovContext()
cova1 = gl.CovAniso.createIsotropic(type=gl.ECov.BESSEL_K,
                                    range = scale, 
                                    param =p1, flagRange= False,
                                    ctxt=ctxt)
cova2 = gl.CovAniso.createIsotropic(type=gl.ECov.BESSEL_K,
                                    range = scale, 
                                    param =p2, 
                                    flagRange= False,ctxt=ctxt)

cova1 = None
   
cov = gl.CovDiffusionAdvection.create(cova1,cova2,scaleTime,vel,sigma2)
covprod = cov.getSpatialTrace()
covprod.display()

Markov (Third Parameter = 1)
- Sill         =      3.000
- Range        =     12.000
- Theo. Range  =      2.000
In [3]:
cov.getSpatialTrace().getMarkovCoeffs()
Out[3]:
array([1., 2., 1.])
In [4]:
cov.getGlobalCorrec()
Out[4]:
0.24989756768264187

Spatial Trace¶

In [5]:
mH1 = 50
mH2 = 50
N=2**8
result = np.array(cov.evalCovFFT([mH1,mH2],0,N).getValues()).reshape(N,N)
trace= np.array(covprod.evalCovFFT([mH1,mH2],N).getValues()).reshape(N,N)
In [6]:
plt.imshow(result,extent = [-mH1,mH1,-mH2,mH2])
plt.colorbar()
plt.show()
plt.imshow(trace,extent = [-mH1,mH1,-mH2,mH2])
c = plt.colorbar()
In [7]:
X1 = np.linspace(-mH1,mH1,N) 
X2 = np.linspace(-mH2,mH2,N) 

plt.plot(X1,trace[int(N/2),:],label = "Trace")
plt.plot(X1,result[int(N/2),:],"--",label="FFT")
l = plt.legend()
In [8]:
N=2**7

ntimes = 40
restime = np.zeros(shape=ntimes)
result = np.zeros(shape = (N*N,ntimes))
times = np.arange(ntimes) * 0.5
for it,time in enumerate(times):
    result[:,it] = np.array(cov.evalCovFFT([mH1,mH2],time,N).getValues())
    restime[it] = result[:,it].reshape(N,N)[int(N/2),int(N/2)]
In [9]:
fig, axs = plt.subplots(nrows=4,ncols=5, figsize=(10,5))
k = 0
for i in range(4) : 
    for j in range(5):
        pt = axs[i,j].imshow(result[:,k].reshape(N,N),origin = "lower",extent = [-mH1,mH1,-mH2,mH2])
        pt.set_clim([result.min(),result.max()])
        k+=1

Temporal covariance¶

In [10]:
plt.plot(times,restime,label = "FFT")
if cova1 is None and np.sum(vel) == 0:
    plt.plot(times,sigma2 * np.exp(-times*scaleTime),"--",label = "Theoretical")

#plt.ylim((0,restime.max()))
ax = plt.legend()