################################################################
############### Examples of Spatial local kriging #############
################################################################
require(GeoModels)
####
model="Gaussian"
# Define the spatial-coordinates of the points:
set.seed(759)
x = runif(1000, 0, 1)
y = runif(1000, 0, 1)
coords=cbind(x,y)
# Set the exponential cov parameters:
corrmodel = "GenWend"
mean=0; sill=1
nugget=0; scale=0.2
param=list(mean=mean,sill=sill,nugget=nugget,smooth=0,
scale=scale,power2=4)
# Simulation of the spatial Gaussian random field:
data = GeoSim(coordx=coords, corrmodel=corrmodel,
param=param)$data
# Maximum pairwise likelihood fitting of the space time random field:
start=list(scale=scale,sill=sill,mean=mean)
fixed=list(power2=4,smooth=0,nugget=0)
fit = GeoFit(data, coordx=coords, corrmodel=corrmodel,
start=start,fixed=fixed,
likelihood='Conditional', type='Pairwise',
neighb=3)
# locations to predict
loc_to_pred=matrix(runif(8),4,2)
################################################################
###
### Example 1. Comparing spatial kriging with local kriging for
### a Gaussian random field with GenWend correlation.
###
###############################################################
param=append(fit$param,fit$fixed)
pr=GeoKrig(fit,loc=loc_to_pred,mse=TRUE)
pr_loc=GeoKrigloc(fit,loc=loc_to_pred,neighb=100,mse=TRUE)
pr$pred;
pr_loc$pred
############################################################
#### Example: spatio temporal Gaussian local kriging ######
############################################################
require(GeoModels)
require(fields)
set.seed(78)
coords=cbind(runif(100),runif(100))
coordt=seq(0,5,0.25)
corrmodel="Matern_Matern"
param=list(nugget=0,mean=0,scale_s=0.2/3,scale_t=0.25/3,sill=2,
smooth_s=0.5,smooth_t=0.5)
data = GeoSim(coordx=coords, coordt=coordt,
corrmodel=corrmodel, param=param)$data
# Maximum pairwise likelihood fitting of the space time random field:
start = list(scale_s=0.2/3,scale_t=0.25,sill=2,mean=0)
fixed = list(smooth_s=0.5,smooth_t=0.5,nugget=0)
I=Inf
lower=list(scale_s=0,scale_t=0,sill=0,mean=-I)
upper=list(scale_s=I,scale_t=I,sill=I,mean=I)
fit = GeoFit(data, coordx=coords, coordt=coordt, model=model, corrmodel=corrmodel,
likelihood='Conditional', type='Pairwise',start=start,fixed=fixed,
optimizer="nlminb",lower=lower,upper=upper,
neighb=3,maxtime=1)
## four location to predict
loc_to_pred=matrix(runif(8),4,2)
## three temporal instants to predict
time=c(0.5,1.5,3.5)
pr=GeoKrig(fit,loc=loc_to_pred,time=time,mse=TRUE)
pr_loc=GeoKrigloc(fit,loc=loc_to_pred,time=time,
neigh=25,maxtime=1, mse=TRUE)
## full and local prediction
pr$pred
pr_loc$pred
############################################################
#### Example: spatio bivariate Gaussian local cokriging ######
############################################################
#set.seed(6)
#NN=1500 # number of spatial locations
#x = runif(NN, 0, 1);
#y = runif(NN, 0, 1)
#coords=cbind(x,y)
## setting parameters
#mean_1 = 2; mean_2= -1
#nugget_1 =0;nugget_2=0
#sill_1 =0.5; sill_2 =1;
### correlation parameters
#CorrParam("Bi_Matern")
#scale_1=0.2/3; scale_2=0.15/3; scale_12=0.5*(scale_2+scale_1)
#smooth_1=smooth_2=smooth_12=0.5
#pcol = -0.4
#param= list(nugget_1=nugget_1,nugget_2=nugget_2,
# sill_1=sill_1,sill_2=sill_2,
# mean_1=mean_1,mean_2=mean_2,
# smooth_1=smooth_1, smooth_2=smooth_2,smooth_12=smooth_12,
# scale_1=scale_1, scale_2=scale_2,scale_12=scale_12,
# pcol=pcol)
## simulation
#data = GeoSim(coordx=coords, corrmodel="Bi_Matern",model=model,param=param)$data
#fixed=list(mean_1=mean_1,mean_2=mean_2, nugget_1=nugget_1,nugget_2=nugget_2,
# smooth_1=smooth_1, smooth_2=smooth_2,smooth_12=smooth_12)
#start=list( sill_1=sill_1,sill_2=sill_2,
# scale_1=scale_1,scale_2=scale_2,scale_12=scale_12, pcol=pcol)
## estimation with maximum likelihood
#fit = GeoFit(data=data,coordx=coords, corrmodel="Bi_Matern",
# likelihood="Marginal",type="Pairwise",optimizer="BFGS",neighb=5,
#start=start,fixed=fixed)
###### co-kriging for the fist component ##############
#xx=seq(0,1,0.022)
#loc_to_pred=as.matrix(expand.grid(xx,xx))
#pr1 = GeoKrigloc(fit,which=1,mse=TRUE,loc=loc_to_pred,neighb=100)
#opar=par(no.readonly = TRUE)
#par(mfrow=c(1,2))
#zlim=c(-2.5,2.5)
#colour = rainbow(100)
#quilt.plot(coords,data[1,] ,col=colour,main = paste(" Fist component"))
#quilt.plot(loc_to_pred,pr1$pred,col=colour,
# main = paste(" Kriging first component"),ylab="")
#par(opar)
Run the code above in your browser using DataLab