Introduction
Presentation of the 2 days workshop
The outline of the course is presented together with a brief history of Geostatistics.
Slides of the presentation
First overview on RGeostats
The main classes and concepts of RGeostats are presented.
Slides of the presentation
Scotland_Temperatures.csv
Geostatistics course
I- Variography
The variogram and its extension (variogram cloud, variogram map) are presented.
Slides of the presentation
RGeostats scripts
II- Estimation
Several linear predictors are compared.
Stationarity assumptions and main variogram models are presented.
Kriging is introduced (simple, ordinary, intrinsic)
+ kriging weights behavior, cross-validation, neighborhood
Slides of the presentation
RGeostats scripts
III- Multivariate Geostatistics
Reminders on bivariate statistics (covariance, linear model, ...)
Multivariate spatial tools (cross-variogram, cross-covariance)
Linear model of coregionalisation
Co-kriging (simple, ordinary)
Simplifications of co-kriging
Collocated co-kriging
Factorial kriging analysis.
Slides of the presentation
RGeostats scripts
IV- Simulations
Why simulations?
Spatial law and Gaussian random function.
Anamorphosis.
Conditional laws.
Gibbs and Sequential Gaussian Simulation (SGS) algorithms.
Turning bands algorithm.
Slides of the presentation
RGeostats scripts
V- Plurigaussian simulations
Categorical simulations
Principles: a single underlying gaussian random function (GRF)
Inference of the parameters of the spatial characteristics
Using one GRF for more than two facies
Using several GRFs
Conditioning
Slides of the presentation
Practice
Illustration of the geostatistical concepts with RGeostats.
Slides
(The scripts written in the slides could be out of date. Please refer to the scripts below which are updated at each new version)
-------------------------------------------------------------------------------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------------------------------------------------------------------------------
Scripts
Require RGeostats > = 10.0.4
Below, you will find the scripts of RGeostats used during the course.
These scripts will be modified in order to stay compatible with the latest version of RGeostats (or to obtain the same results in a more elegant manner).
I) Variography
Return to the table of content
I-1) Preamble
Load the data and change the name (to have a shorter name) and make the variable 4 (Elevation) inactive
- Code: Select all
data(Exdemo_Scotland_Temperatures)
dat=Exdemo_Scotland_Temperatures
dat=db.locate(dat,4)
Note that the last line is equivalent to
- Code: Select all
dat=db.locate(dat,4,"NA")
or
- Code: Select all
dat=db.locate(dat,"Elevation","NA")
I-2) Exploratory data analysis
Display information on the data base
- Code: Select all
dat
Display the data
- Code: Select all
dat[]
Histogram of the 5th column (temperatures)
- Code: Select all
hist(dat[,5])
Some statistics (e.g the maximum of the temperature)
- Code: Select all
db.stat(dat,fun="maxi",name=5)
Plot
- Code: Select all
plot(dat)
With the same scale for abscissa and ordinate.
- Code: Select all
plot(dat,asp=1)
I-3) Compute and plot the experimental variogram
10 lag of size 40, omnidirectional:
- Code: Select all
vario=vario.calc(dat,lag=10,nlag=40)
Various plots:
- Code: Select all
plot(vario)
plot(vario,npairdw=TRUE)
plot(vario,npairdw=TRUE,npairpt=TRUE)
Variogram cloud
- Code: Select all
v_cloud=cloud.calc(dat)
plot(v_cloud)
I-4) Fit a model
- Code: Select all
model=model.auto(vario)
The resulting 'model' is an object of type model
List of all the elementary models
- Code: Select all
melem.name()
To select a subset
- Code: Select all
melem.name(c(1,2,3))
Specify the elementary models for the automatic fitting
- Code: Select all
model.auto(vario,struct=melem.name(c(1,2,3)))
Change the order. The result can be different.
- Code: Select all
model.auto(vario,struct=melem.name(c(1,3,2)))
Indeed, it changes the initial guess.
- Code: Select all
model.auto(vario,struct=melem.name(c(1,3,2)),maxiter=0)
model.auto(vario,struct=melem.name(c(1,2,3)),maxiter=0)
model.auto(vario,struct=melem.name(c(1,2,3)),flag.noreduce=T)
With constraints
a) Display the result with verbosity to know the list of parameters and their rank
- Code: Select all
model.auto(vario,struct=melem.name(c(1,2,3)),flag.noreduce=T,verbose=T)
b) Then put your constraints with 'lower' and 'upper' arguments
You can put equality and inequality constraints.
- Code: Select all
model.auto(vario,struct=melem.name(c(1,2,3)),flag.noreduce=T,lower=c(0.02,NA,NA,NA,NA),upper=c(0.02,NA,NA,NA,NA))
Display the result:
- Code: Select all
model
I-5) 4 directions
We can compute directional variogram by specifying calculation directions.
In 2D, we provide the angles in degrees.
- Code: Select all
vario_4dir=vario.calc(dat,lag=15,nlag=15,dir=c(0,45,90,135))
vario_4dir
m4dir=model.auto(vario_4dir,struct=melem.name(c(1,3,2)))
We can display the results
- Code: Select all
plot(vario_4dir,npairdw=TRUE)
plot(m4dir,add=T,vario=vario_4dir)
I-7) Variogram map
- Code: Select all
vm=vmap.calc(dat,nx=10,dx=15)
model=vmap.auto(vm,struct=melem.name(c(1,3,2)),
flag.noreduce=T,lower=c(0.02,NA,NA,NA,NA),upper=c(0.02,NA,NA,NA,NA))
vario_4dir=vario.calc(dat,lag=15,nlag=15,dir=c(0,45,90,135))
plot(vario_4dir)
plot(model,vario=vario_4dir,add=T)
I-8) Variogram on grid
- Code: Select all
data(Exdemo_Scotland_Elevations)
grid=Exdemo_Scotland_Elevations
- Code: Select all
plot(grid,asp=1,col=topo.colors(100),pos.legend=1)
The use of vario.grid is much more efficient than vario.calc. But it is dedicated to data on grid (not for isolated points)
and the variogram is only computed along the grid axis (with no angular tolerance).
- Code: Select all
vario=vario.grid(grid,nlag=30)
plot(vario)
model.auto(vario,struct=melem.name(c(1,3,2)))
To force an isotropic model:
- Code: Select all
model.auto(vario,struct=melem.name(c(1,3,2)),auth.aniso=F)
II) Estimation
Return to the table of content
Will come soon.
III) Multivariate
Return to the table of content
III-1) Preamble
- Code: Select all
data(Exdemo_Scotland_Temperatures)
data(Exdemo_Scotland_Elevations)
dat=Exdemo_Scotland_Temperatures
grid=Exdemo_Scotland_Elevations
Creation of a variable indicating if temperature is available or if only the elevation is measured.
The argument 'loctype="NA"' indicates that the locator is NA.
Therefore, the variable is in the data base but currently inactive.
- Code: Select all
dat=db.add(dat,"Available_Temp"=(!is.na(dat[,5])),loctype="NA")
Neighborhood (unique)
- Code: Select all
unique.neigh=neigh.init(type=0,ndim=2)
III-2) Exploratory data analysis
We make the variable "Available_Temp" active as a selection.
- Code: Select all
dat=db.sel(dat,selold=6)
Therefore, the next commands will only concern data points with a temperature station
- Code: Select all
plot(dat)
db.stat(dat,fun="maxi",names=4)
Now we invert the selection. Therefore, the next commands will only concern data point without a temperature station.
- Code: Select all
dat[,6]=!dat[,6]
plot(dat,add=T,col=3)
db.stat(dat,fun="maxi",names=4)
Conclusion: the stations (giving the temperatures) are at a lower elevation that what we can find in Scotland.
Finally, we come back to the initial values of the variable "Available_Temp".
- Code: Select all
dat[,6]=!dat[,6]
Bivariate statistics
- Code: Select all
correlation(dat,icol1=4,icol2=5,col=c(2,3),pos.legend=1)
regression(dat,icol1=4,icol2=5,save.coeff= T,flag.draw=T)
III-3) Cross-variogram model
(same commands as for univariate case)
- Code: Select all
cross_vario=vario.calc(dat,lag=10,nlag=40)
m_model=model.auto(cross_vario,struct=melem.name(c(1,3,2)),flag.noreduce=T)
III-4) Co-kriging
a) Isotopic
We make the selection active. So only the pairs containing the temperature AND the elevations will be considered (isotopic case).
- Code: Select all
dat=db.locate(dat,6,"sel")
We perform the co-kriging (note that it provides results for all the variables).
- Code: Select all
resultm1=kriging(dat,grid,model=m_model,neigh=unique.neigh)
resultm1
plot(resultm1,name="Kriging.January_temp.estim",asp=1,pos.legend=5,col=topo.colors(100))
plot(dat,add=T)
plot(resultm1,name="Kriging.January_temp.stdev",asp=1,pos.legend=5,col=topo.colors(100))
plot(dat,add=T)
b) Heterotopic
We make the selection inactive.
- Code: Select all
dat=db.locate(dat,6)
- Code: Select all
resultm2=kriging(dat,grid,model=m_model,neigh=unique.neigh)
plot(resultm2,name="Kriging.January_temp.estim",asp=1,pos.legend=5,col=topo.colors(100))
plot(dat,add=T)
plot(resultm2,name="Kriging.January_temp.stdev",asp=1,pos.legend=5,col=topo.colors(100))
plot(dat,add=T)
c) Collocated co-kriging
(only works with a moving neighborhood. It will be extended to unique neighborhood soon)
- Code: Select all
dat=db.locate(dat,4:5,"z")
dat=db.locate(dat,6)
- Code: Select all
moving.neigh=neigh.input(ndim=2)
#Answers 2, 5, 20, n, n, 50
- Code: Select all
resultm3=kriging(dat,grid,model=m_model,neigh=moving.neigh,flag.colk=T)
plot(resultm3,name="Kriging.January_temp.estim", asp=1,pos.legend=5,col=topo.colors(100))
III-5) Kriging with an external drift
Put the variable Elevations with a locator "f" both in the data (dat) and in the grid of target locations (grid)
- Code: Select all
dat=db.locate(dat,4,"f")
grid=db.locate(grid,4,"f")
Perform the regression and put the result (residuals) in a new db: dat_r.
- Code: Select all
dat_r=regression(dat,flag.bivar=FALSE)
Compute and model the variogram of the residuals.
- Code: Select all
vario=vario.calc(dat_r,nlag=40,lag=10)
mres=model.auto(vario,struct=melem.name(c(5,3)),upper=c(NA,NA,NA,300))
Remark: compare the sill of the variable Temperature and the one of the residual.
Perform the kriging with external drift (see new parameter "uc")
- Code: Select all
resultm4=kriging(dat,grid,model=mres,uc=c("1","f1"),neigh=unique.neigh)
Display the results
- Code: Select all
plot(resultm4,name="Kriging.January_temp.estim",asp=1,pos.legend=5,col=topo.colors(100))
plot(dat,add=T,name=5)
Comparison (Collocated - External drift)
- Code: Select all
correlation(resultm3,7,6,resultm4,flag.diag=T)
IV) Simulations
Return to the table of content
Preamble
- Code: Select all
data(Exdemo_Scotland_Temperatures)
data(Exdemo_Scotland_Elevations)
dat=Exdemo_Scotland_Temperatures
grid=Exdemo_Scotland_Elevations
dat=db.locate(dat,4)
- Code: Select all
vario=vario.calc(dat,lag=10,nlag=40)
model=model.auto(vario, flag.noreduce=T,struct=melem.name(c(1,2,3)),lower=c(0.02,NA,NA,NA,NA),upper=c(0.02,NA,NA,NA,NA))
unique.neigh=neigh.init(type=0,ndim=2)
Compute a mean
- Code: Select all
m=mean(dat[,5],na.rm=T)
100 Conditional simulations
- Code: Select all
simu = simtub(dat,grid,model=model,mean=m,uc="",neigh=unique.neigh,nbtuba=1000,nbsimu=100)
Display some simulations
- Code: Select all
plot(simu,asp=1,col=topo.colors(100),name=6,zlim=c(-2,5),pos.legend=5)
plot(simu,asp=1,col=topo.colors(100),name=7,zlim=c(-2,5),pos.legend=5)
plot(simu,asp=1,col=topo.colors(100),name=8,zlim=c(-2,5),pos.legend=5)
Compute the pointwise means of all the simulations
- Code: Select all
result=db.compare(simu,fun="mean",name=6:105)
Comparison with simple kriging
- Code: Select all
result0=kriging(dat,grid,model=model,neigh=unique.neigh,uc="",mean=m)
correlation(result,106,6,result0,flag.diag=T)