GeoENV 2014: RGeostats Workshop

Events related to RGeostats are posted here (Training courses, Conferences, Publications or Congresses communication...)

GeoENV 2014: RGeostats Workshop

Postby Nicolas Desassis » Mon Jul 14, 2014 8:02 pm

You will find here all the documents related to the 2 days workshop on RGeostats held in Paris the 7 and 8th July 2014.

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)
Nicolas Desassis
 
Posts: 198
Joined: Thu Sep 20, 2012 4:22 pm

Return to Events

Who is online

Users browsing this forum: No registered users and 1 guest

cron