The following tutorial gives an introduction on how to instantiate variogram objects with a theoretical model. The examples used were Matern models explained in Chapter 3 of Model-Based Geostatistics (Diggle, Ribeiro, 2007). Spystats can be a Python alternative to the popular R package geoR.
In [2]:
# Load Biospytial modules and etc.
%matplotlib inline
import sys
sys.path.append('/apps')
sys.path.append('/apps/external_plugins/spystats')
import django
django.setup()
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
## Use the ggplot style
plt.style.use('ggplot')
In [3]:
from spystats import tools
from scipy.special import logit,expit
from scipy.stats import bernoulli,poisson
Here we will show how to plot the standard correlation functions. Currently spystats.tools
have implemented the following models:
* Matern Model tools.MaternVariogram()
* Exponential Model tools.ExponentialVariogram()
* Gaussian Model tools.GaussianVariogram()
* Spherical Model tools.SphericalVariogram()
* Whittle Model tools.WhittleVariogram()
We will plot a function similar to the figure 3.2 (Diggle and Ribeiro,2007)
In [3]:
## Create the domain X
x = np.linspace(0,1,101)
In [4]:
m1 = tools.MaternVariogram(sill=1,range_a=0.25,kappa=0.5)
m2 = tools.MaternVariogram(sill=1,range_a=0.16,kappa=1.5)
m3 = tools.MaternVariogram(sill=1,range_a=0.13,kappa=2.5)
In [5]:
plt.plot(x,m1.corr_f(x))
plt.plot(x,m2.corr_f(x))
plt.plot(x,m3.corr_f(x))
plt.legend(['$\kappa$ = 0.5 $\phi$ = 0.25' ,'$\kappa$ = 1.5 $\phi$ = 0.16','$\kappa$ = 2.5 $\phi$ = 0.13'])
Out[5]:
In [6]:
#We need to create a domain now for two dimensions.
grid = tools.createGrid()
## Apply the simualtion to all the models
%time s = map(lambda m : tools.simulateGaussianRandomField(m,grid,random_seed=12345),[m1,m2,m3])
In [7]:
plt.imshow(s[0],interpolation=None)
Out[7]:
In [8]:
plt.imshow(poisson.rvs(np.exp(s[0]),0))
Out[8]:
In [9]:
plt.imshow(bernoulli.rvs(expit(s[0]),1))
Out[9]:
In [10]:
plt.imshow(s[1],interpolation=None)
Out[10]:
In [11]:
plt.imshow(poisson.rvs(np.exp(s[1]),0))
Out[11]:
In [12]:
plt.imshow(bernoulli.rvs(expit(s[1]),0))
Out[12]:
In [13]:
plt.imshow(s[2],interpolation=None)
Out[13]:
In [14]:
plt.imshow(poisson.rvs(np.exp(s[2]),1))
Out[14]:
In [15]:
plt.imshow(bernoulli.rvs(expit(s[2]),0))
Out[15]:
In [16]:
phi_A = (np.pi / 3.0)
#phi_A = 0.0
phi_R = 4.0
a = np.array([[np.cos(phi_A), -np.sin(phi_A)],[np.sin(phi_A),np.cos(phi_A)]])
b = np.array([[1, 0],[0,1.0 / phi_R]])
In [17]:
A = np.matmul(a,b)
In [18]:
new_coords = np.dot(grid[['Lon','Lat']].values,A)
n = len(new_coords)
nc = pd.DataFrame({'nLon':new_coords[:,0],'nLat':new_coords[:,1],'Y': np.zeros(n)})
nc = tools.toGeoDataFrame(nc,'nLon','nLat')
In [19]:
nc.plot()
Out[19]:
In [20]:
vg = tools.Variogram(nc,'Y',model=m2)
In [21]:
Sigma = vg.calculateCovarianceMatrix()
In [22]:
from scipy.stats import multivariate_normal as mvn
np.random.seed(12345)
sim1 = mvn.rvs(mean=nc.Y,cov=Sigma)
n_2 = int(np.sqrt(n))
In [23]:
plt.imshow(sim1.reshape(n_2,n_2))
Out[23]:
Therefore, it has the functions of Geometric Manipulations. (See: http://geopandas.org/geometric_manipulations.html?highlight=affine)
In [24]:
transgrid = grid.geometry.rotate(-phi_A, origin=(0,0),use_radians=True).scale(1,1.0 / phi_R,origin=(0,0))
transgrid.plot()
Out[24]:
If you go back (up) to the manually transformed coordinates you will find out that the plot is the same. Except for the argument phi_R that, for the geodataframe object/method (shapely) it should be -phi_R. Remmember also to use origin=(something). Without this the method will not work.
In [25]:
grid['transcoords'] = transgrid
grid = grid.set_geometry('transcoords')
In [26]:
grid = grid.set_geometry('transcoords')
In [27]:
transfield= tools.simulateGaussianRandomField(grid=grid,variogram_model=m1,random_seed=12345)
plt.imshow(transfield)
Out[27]:
In [28]:
plt.imshow(poisson.rvs(np.exp(transfield),0))
plt.colorbar()
Out[28]:
In [29]:
plt.imshow(bernoulli.rvs(expit(transfield),0))
plt.colorbar()
Out[29]:
In [30]:
transfield= tools.simulateGaussianRandomField(grid=grid,variogram_model=m2,random_seed=12345)
plt.imshow(transfield)
Out[30]:
In [31]:
plt.imshow(poisson.rvs(np.exp(transfield),0))
plt.colorbar()
Out[31]:
In [32]:
plt.imshow(bernoulli.rvs(expit(transfield),0))
plt.colorbar()
Out[32]:
In [33]:
transfield= tools.simulateGaussianRandomField(grid=grid,variogram_model=m3,random_seed=12345)
plt.imshow(transfield)
Out[33]:
In [34]:
plt.imshow(poisson.rvs(np.exp(transfield),0))
plt.colorbar()
Out[34]:
In [35]:
plt.imshow(bernoulli.rvs(expit(transfield),0))
plt.colorbar()
Out[35]:
In [36]:
#grid = grid.set_geometry('geometry')
p = np.pi / 5.0
transgrid = grid.geometry.skew(-p,origin=(0,0),use_radians=True)
transgrid.plot()
Out[36]:
In [37]:
grid['transcoords'] = transgrid
grid = grid.set_geometry('transcoords')
In [38]:
transfield= tools.simulateGaussianRandomField(grid=grid,variogram_model=m1,random_seed=12345)
plt.imshow(transfield)
Out[38]:
In [39]:
plt.imshow(poisson.rvs(np.exp(transfield),0))
plt.colorbar()
Out[39]:
In [40]:
plt.imshow(bernoulli.rvs(expit(transfield),0))
plt.colorbar()
Out[40]:
In [41]:
transfield= tools.simulateGaussianRandomField(grid=grid,variogram_model=m2,random_seed=12345)
plt.imshow(transfield)
Out[41]:
In [42]:
plt.imshow(poisson.rvs(np.exp(transfield),0))
plt.colorbar()
Out[42]:
In [43]:
plt.imshow(bernoulli.rvs(expit(transfield),0))
plt.colorbar()
Out[43]:
In [44]:
transfield= tools.simulateGaussianRandomField(grid=grid,variogram_model=m3,random_seed=12345)
plt.imshow(transfield)
Out[44]:
In [45]:
plt.imshow(poisson.rvs(np.exp(transfield),0))
plt.colorbar()
Out[45]:
In [46]:
plt.imshow(bernoulli.rvs(expit(transfield),0))
plt.colorbar()
Out[46]:
In [59]:
## Simulations with non squared grid
grid = tools.createGrid(grid_sizex=50,minx=-1,maxx=2,miny=-1,maxy=2,grid_sizey=70)
In [60]:
sim1= tools.simulateGaussianRandomField(grid=grid,variogram_model=m3,random_seed=12345,return_matrix=False)
In [61]:
grid = pd.concat([grid,sim1],axis=1)
In [63]:
fig = plt.figure(figsize=(16,16), dpi= 80, facecolor='w', edgecolor='w')
#grid.plot(column='sim')
plt.scatter(grid.Lon,grid.Lat,c=grid.sim.values,s=150)
plt.colorbar()
Out[63]:
In [64]:
a,b,c = tools.simulatedGaussianFieldAsPcolorMesh(m2,grid_sizex=50,minx=-1,maxx=2,miny=-1,maxy=2,grid_sizey=70)
In [97]:
fig = plt.figure(figsize=(16,16), dpi= 80, facecolor='w', edgecolor='w')
#grid.plot(column='sim')
alpha = 0.5
#CS = plt.contour(a,b,c,colors='white',lw=2.0)
#plt.clabel(CS, fontsize=12, inline=1)
plt.pcolormesh(a,b,c)
#plt.scatter(grid.Lon,grid.Lat,c=np.exp(alpha + grid.sim),s=150)
plt.colorbar()#CS = plt.contour(a,b,c,c='k',lw=2.0)
Out[97]:
In [122]:
fig = plt.figure(figsize=(16,16), dpi= 80, facecolor='w', edgecolor='w')
#grid.plot(column='sim')
alpha = 0.5
CS = plt.contour(a,b,np.exp(alpha + c),colors='white',lw=2.0)
plt.clabel(CS, fontsize=12, inline=1)
plt.pcolormesh(a,b,poisson.rvs(np.exp((alpha * np.random.normal() + c ))))
#plt.scatter(grid.Lon,grid.Lat,c=np.exp(alpha + grid.sim),s=150)
plt.colorbar()
Out[122]:
In [132]:
fig = plt.figure(figsize=(16,16), dpi= 80, facecolor='w', edgecolor='w')
#grid.plot(column='sim')
alpha = 7
CS = plt.contour(a,b,np.exp(alpha + c),colors='white',lw=2.0)
plt.clabel(CS, fontsize=12, inline=1)
plt.pcolormesh(a,b,poisson.rvs(np.exp((alpha * np.random.normal() + c ))))
#plt.scatter(grid.Lon,grid.Lat,c=np.exp(alpha + grid.sim),s=150)
plt.colorbar()
Out[132]:
In [170]:
import scipy
import matplotlib
In [180]:
np.arange(-2.5,3,1)
Out[180]:
In [155]:
scipy.stats.binom.rvs(5,p)
Out[155]:
In [193]:
fig = plt.figure(figsize=(16,16), dpi= 80, facecolor='w', edgecolor='w')
#grid.plot(column='sim')
alpha = 7
n_sp = 8
binom_p = expit(c)
binom_v = scipy.stats.binom.rvs(n_sp,binom_p)
cmap=plt.cm.Accent
norm = matplotlib.colors.BoundaryNorm(np.arange(1,n_sp+1), cmap.N)
CS = plt.contour(a,b,np.exp(alpha + c),colors='white',lw=2.0)
plt.clabel(CS, fontsize=12, inline=1)
plt.pcolormesh(a,b,binom_v,cmap=cmap,edgecolor='none',norm=norm)
#plt.scatter(grid.Lon,grid.Lat,c=np.exp(alpha + grid.sim),s=150)
plt.colorbar(ticks=np.linspace(1.5,n_sp+0.5,n_sp))
#cmap=plt.cm.rainbow
#norm = matplotlib.colors.BoundaryNorm(np.arange(-2.5,3,1), cmap.N)
#plt.scatter(x,y,c=z,cmap=cmap,norm=norm,s=100,edgecolor='none')
Out[193]: