vario.pgs {RGeostats}R Documentation

Calculate the variogram of the Underlying GRFs from the categories (e.g lithofacies)

Description

Calculate the variogram of the underlying GRFs from the categories (e.g lithofacies)

Usage

vario.pgs(db, dbprop = NA, rule = rule.input(), props = NA,
	  model.type = "residuals", flag.rho = FALSE, 
          dirvect = NA, lag = NA, toldis = 0.5, tolang = NA, 
          bench=0, cylrad=0, nlag = NA, breaks = NA, calcul = "cov",
          opt.code = 0, tolcode = 0, dates=NA, vario.add=NA) 

Arguments

db

The db-class containing the facies information (stored in the "z1" variable) used to calculate the experimental variogram.

dbprop

The db-class where the proportion variables must be located in the non-stationary case. If not defined, dbprop is set to the input Db.

rule

The rule-class structure which defines the Lithotype rule

props

Array giving the (constant) proportions of the different facies involved in the Lithotype rule

model.type

Character string specifying the model structure to use. One of "full", "symetrical" or "residuals" (default). See details.

flag.rho

When TRUE, the correlation coefficient linking the two underlying Gaussian Random Functions is inferred. Otherwise, the correlation coefficient as defined in the input rule-class is kept unchanged. Only flag.rho=FALSE is currently available.

dirvect

The directions in which the variogram must be calculated. For more information, refer to get.directions

lag

Array containing the distance lags for each calculation direction. If no corresponding lag is not defined, a default lag is calculated so that the maximum distance is equal to half of the field diagonal. For one distance, if the lag is not defined, it is set to the default lag.

toldis

Tolerance on the distance, expressed as a proportion of the lag. This parameter is only used in the case of regular lags

tolang

Array containing the angular tolerance for each calculation direction. If no angular tolerance is not defined, a default tolerance is calculated as one half of the angle. For one direction, if the angular tolerance is not defined, it is set to the default tolerance.

bench

Array containing the bench separation tolerance for each calculation direction. When bench is defined and positive, a pair is discarded if the separation distance between the two points along the third axis is larger than 'bench'.

cylrad

Array containing the cylrad tolerance for each calculation direction. When specified and positive, a pair is discarded if the separation distance in the plane orthogonal to the calculation distance is larger than 'cylrad'.

nlag

Array containing the number of lags for each calculation direction. If not defined, the defaulted number of lags is set to 10. For one direction, if the number of lags is not defined, it is set to 10.

breaks

Vector giving the boundaries of irregular distance lags. For multidrectional computations, breaks is either: a vector (The breaks are then the same for all directions) a list of vectors of the breaks for each direction. If provided, nlag is adapted to the breaks. Intervals are open left and closed right (e.g ]0;15]).

calcul

Character string giving the type of structure to compute:

  • vg : Variogram

  • cov : Covariance (Centered)

  • covnc : Non-centered Covariance

opt.code

Array containing the option concerning the use of the code for each direction. This option defines the strategy for the code usage when considering a pair of samples:

  • 0 : a pair of samples is selected whatever their codes (if active)

  • 1 : a pair is retained if the code is active in the data base and the codes of the two samples are closer than tolcode

  • 2 : a pair is retained only if the codes of the two samples are different

tolcode

Array giving the maximum distance between the codes of the two samples for each direction.

dates

When defined and if the 'time' locator is defined, this value corresponds to the bounds of the Date Intervals. It is organized as follows:

  • The lower bound of the first Date Interval

  • The upper bound of the first Date Interval

  • The lower bound of the second Date Interval

  • The upper bound of the second Date Interval

  • ...

vario.add

When defined, the current variogram computation is concatenated to the contents of the current variogram.

Details

The model.type argument indicates the structure of the multivariate covariance matrix to fit. The "residuals" option means that the second GRF is supposed to be of the form:

Y_2 = ρ Y_1 + √{1-ρ^2} R

where Y_1 is the first GRF and R is a second GRF independent from Y_1.

The "symmetrical" option indicates that the resulting covariance function is such as C12(h)=C21(h).

The "full" option indicates that there is no constraint.

The "residual" option has to be used if the two GRF's are supposed to be uncorrelated (it is then specified by the field "rho" of the input rule which is equal to 0 and by the input flag.rho which is set to FALSE).

Value

The output is the experimental Covariance of the GRFs.

The application also calculates the scores during this iterative procedure. The stack of these scores can be displayed through a call to the keypair mechanism after the function 'vario.pgs' has been run successfully. The syntax is:

stack <- get.keypair("vario.pgs_stack")

The resulting object 'stack' is organized as a matrix with the following information described per column:

Examples


## Not run: 
#This part contains several independent examples.
# Example 1: independent bigaussian model with stationary proportions
# Example 2: correlated bigaussian model ("residuals" type) 
  	     with estimation of the correlation coefficient
# Example 3: correlated bigaussian model ("symmetrical" type). 
  	     The interest is the particular form of the 
	     cross-correlation.
# Example 4: Shifted bigaussian ("full" type).
# Example 5: Non-stationary proportions

#########################################
#              Example 1		#
#########################################

#Simulation of 2 independent GRFs and one categorical variable:

set.seed(234)
data(Exdemo_simpgs_std.rule)
db=db.create(x1=runif(1000),x2=runif(1000),ndim=2)
model=model.create(vartype="Gaussian",range=.2,ndim=2)
model2=model.create(vartype="Exponential",range=.3,ndim=2)   
db=simtub(dbout=db,model=model,nbtuba=2000)
vario1=vario.calc(db,nlag=30)
db=simtub(dbout=db,model=model2,nbtuba=2000)
vario2=vario.calc(db,nlag=30)
db=db.locate(db,4:5,"z")
db=db.rule(db,rule=Exdemo_simpgs_std.rule)

#Call vario.pgs
vario=vario.pgs(db,rule=Exdemo_simpgs_std.rule,nlag=30)

#Graphical representation of the Rule
plot(Exdemo_simpgs_std.rule)
title("Composition rule")
dev.off()

#Graphical representation of the data and the variograms
split.screen(c(2,2))

#The data
screen(1)
plot(db[,2],db[,3],col=db[,6],pch=16,xlab="",ylab="",cex=.6)
title("Data: 3 categories")

#Legend for the variograms
screen(2)
plot(c(0,10),c(0,10),cex=0,xaxt="n",yaxt="n",xlab="",ylab="")
legend(0,5,legend=c("Models","Vario from GRFs","Vario from Categories"),
		col=c(6,4,1),lty=c(1,1,NA),pch=c(NA,NA,1),box.col=0)
title("Legend for the variograms")

#GRF1
screen(3)
v1=vario.reduce(vario=vario,varcols=1,flag.vario=TRUE)
plot(vario1,col=4,xlab="h",ylab=expression(gamma),reset=FALSE)
points(v1[1]$hh,v1[1]$gg)
plot(model,vario=v1,add=TRUE,col=6)
title("Variograms of the GRF #1")

#GRF2
screen(4)
v2=vario.reduce(vario=vario,varcols=2,flag.vario=TRUE)
plot(vario2,col=4,xlab="h",ylab=expression(gamma),reset=FALSE)
points(v2[1]$hh,v2[1]$gg)
plot(model2,vario=v2,pch=1,add=TRUE,col=6)
title("Variograms of the GRF #2")

rm(db,Exdemo_simpgs_std.rule,model,model2,v1,v2,vario,vario1,vario2)
par(mfrow=c(1,1))

#########################################
#              Example 2		#
#########################################


set.seed(23)
data(Exdemo_simpgs_std.rule)
db=db.create(flag.grid=TRUE,nx=5000,ndim=1)
r=-.7
model=model.create(vartype="Gaussian",range=20,sill=matrix(c(1,r,r,r^2),2),
ndim=1)
model=model.create(vartype="Exponential",range=50,sill=matrix(c(0,0,0,1-r^2),2),
model=model,ndim=1)   
Exdemo_simpgs_std.rule$rho=r
db=simtub(dbout=db,model=model,nbtuba=2000)
db=db.rule(db,rule=Exdemo_simpgs_std.rule)

#1st call of vario.pgs to determine rho

vario=vario.pgs(db,rule=Exdemo_simpgs_std.rule,nlag=5,lag=1,flag.rho=T)
Exdemo_simpgs_std.rule$rho=vario$vars[2]

#Second call with fixed and estimated rho

Exdemo_simpgs_std.rule$rho=vario$vars[2]
vario=vario.pgs(db,rule=Exdemo_simpgs_std.rule,nlag=50,lag=1,flag.rho=F)


plot(model,as.cov=TRUE,vario=vario)
plot(vario,add=TRUE,col=2)

plot(db[,2],db[,3],col=db[,6],pch=16,xlab="",ylab="",cex=.6)

#########################################
#              Example 3		#
#########################################



r1=-1
r2=0
s11=.4
s12=.4
range1=500

s21=.6
s22=.6
range2=10000

p1=.8
p2=p3=.09
p4=1-p1-p2-p3

rule=rule.create(c("S","T","F1","F2","T","F3","F4"))
model=model.create(melem.name(5),
sill=matrix(c(s11,sqrt(s11*s12)*r1,sqrt(s11*s12)*r1,s12),2),range=range1,ndim=2)
model=model.create(melem.name(5),
sill=matrix(c(s21,sqrt(s21*s22)*r2,sqrt(s21*s22)*r2,s22),2),range=range2,
model=model,ndim=2)
r=model.eval(model,0,as.cov=T,jvar=2)
rule$rho=r
curve(model.eval(model,x,as.cov=T,ivar=1,jvar=2,),0,6000,ylab="",
ylim=c(-1,1),n=1000)

db=db.create(nx=c(200,200),dx=c(100,100))
db=simtub(dbout=db,model=model,nbtuba=1000)
db=db.rule(db,rule=rule,props=c(p1,p2,p3,p4))
plot(db[,2],db[,3],pch=16,col=db[,6])

vario=vario.pgs(db,rule=rule,lag=1,nlag=50,flag.rho=F,model.type="symmetrical",props=c(p1,p2,p3,p4))
plot(vario)
plot(model,vario=vario,add=T,col=2)


#########################################
#              Example 4		#
#########################################


set.seed(2233)
data(Exdemo_simpgs_std.rule)
model=model.create(melem.name(2),range=50,ndim=1)
data(Exdemo_simpgs_std.rule)
Exdemo_simpgs_std.rule$rho=model.eval(model,10,as.cov=T)
db=db.create(nx=5010)
db=simtub(dbout=db,model=model,nbtuba=2000)
db=db.create(nx=5000,z1=db[1:5000,3],z2=db[11:5010,3])
db=db.rule(db,rule=Exdemo_simpgs_std.rule)
Exdemo_simpgs_std.rule$rho=0
vario=vario.pgs(db,rule=Exdemo_simpgs_std.rule,lag=1,nlag=5,flag.rho=T,model.type="full")
vario=vario.pgs(db,rule=Exdemo_simpgs_std.rule,lag=1,nlag=50,flag.rho=F,model.type="full")
plot(vario)
x11()
plot(vario[1]$hh[103:201],vario[1]$gg[103:201],xlab="h",ylab=expression(gamma),type="l")
abline(v=0)
curve(model.eval(model,x+10,as.cov=T),add=T,col=2)


#########################################
#              Example 5		#
#########################################



set.seed(23)
data(Exdemo_simpgs_std.rule)
db=db.create(nx=5000)
r=-.7
model=model.create(vartype="Gaussian",range=20,
sill=matrix(c(1,r,r,r^2),2))
model=model.create(vartype="Exponential",range=50,
sill=matrix(c(0,0,0,1-r^2),2),model=model)   
Exdemo_simpgs_std.rule$rho=r
db=simtub(dbout=db,model=model,nbtuba=2000)

#Simulation of the proportions
model.prop=model.create(vartype="Gaussian",range=1000)
db.prop=db.create(nx=5000)
db.prop=simtub(dbout=db.prop,model=model.prop,seed=12213,nbtuba=2000,nbsimu=3)
proportions=exp(db.prop[,3:4])/(apply(exp(db.prop[,3:4]),1,sum)+1)
proportions=cbind(proportions,1-apply(proportions[,1:2],1,sum))

db=db.add(db,proportions,locnames=c("p1","p2","p3"))
db=db.rule(db,rule=Exdemo_simpgs_std.rule)


par(mfrow=c(1,5))
plot(c(0,1),c(1,5000),cex=0,xlab="Proportions",ylab="",xaxp=c(0,1,2),yaxp=c(0,5000,5),cex.lab=1.2,cex.axis=1.4)
polygon(c(0,db[,5],0),c(1,1:5000,5000),col=1)
polygon(c(db[5000:1,5],1-db[,7]),c(5000:1,1:5000),col=2)
polygon(c(1-db[,7],1,1),c(1:5000,5000,1),col=3)

db.prop=db.prop.thresh(db.prop,db,rule=Exdemo_simpgs_std.rule,model=model)
plot(db[,3],1:5000,type="l",xlab="1st GRF",ylab="",xaxp=c(3,3,2),yaxt="n",cex.lab=1.2,cex.axis=1.4)
lines(db.prop[,6],1:5000,col=2,lwd=1.5)

plot(c(0,1),c(1,5000),cex=0,xlab="Facies",ylab="",xaxt="n",yaxt="n",cex.axis=1.4)
from=c(1,.5+(1:4999)[db[-1,8]!=db[-5000,8]])
to=c(from[-1],5000)
for(i in 1:length(from))
polygon(c(0,1,1,0),c(from[i],from[i],to[i],to[i]),col=db[1,8],db[-1][db[-1,8]!=db[-5000,8],8][i])
vario=vario.pgs(db,rule=Exdemo_simpgs_std.rule,lag=1,nlag=50,flag.rho=F)
plot(vario)
plot(model,vario=vario,add=T,col=2)

## End(Not run)

[Package RGeostats version 11.1.1 Index]