#######
# ICES CRR Handbook of geostatistics in R for fisheries and marine ecology
#
# This code performs Global estimation over a domain (polygon) 
# It uses the datasets Demo.herreggs.scot
#
# The data were supplied by Marine Scotland Science at the Marine Laboratory,
# Aberdeen, UK. 
#
# Author: P.Petitgas, Ifremer 
#######
{
library(RGeostats)
liste.initiale = ls()
on.exit(demo.clean(liste.initiale))
par(ask=FALSE)
demo.start(FALSE)

######
### Loading Data
######
                                        # limits of study area 
x1lim<-25.9; x2lim<-26.5; y1lim<-17.0; y2lim<-17.58 
                                        # read data file and polygon
rg.load("Demo.herreggs.scot.db.data","db.data")
rg.load("Demo.herreggs.scot.poly.data","poly.data")
projec.toggle(0,0)

                                        # select points inside polygon
db.data <- db.polygon(db.data,poly.data)
                                        # plot data and polygon
plot(db.data,name.prop="z1",xlab="",ylab="",
     xlim=c(x1lim,x2lim),ylim=c(y1lim,y2lim))
plot(poly.data,add=T,lty=1,density=0)
demo.pause("Data")
                                        # Data mean and variance inside polygon
zm <- mean(db.data[,4][db.data[,5]]) 
zv <- var(db.data[,4][db.data[,5]])*(sum(db.data[,5])-1)/sum(db.data[,5])
cat("Data\n")
cat("mean: ",zm,"    std: ",sqrt(zv),"   cv: ",sqrt(zv)/zm,"\n")

######
### Calculate experimental Variogram
######
cat("Defining the model...\n")

                                        # Calculate experimental variogram
Lag <- 0.05; Nlag <- 9  
vg <- vario.calc(db.data,dirvect=0,lag=Lag,toldis=0.5, nlag=Nlag, breaks = NA, 
                 calcul="vg",by.sample=FALSE,opt.code=0, tolcode=0, means=NA)

                                        # Plot experimental variogram
vario.plot(vg,npairdw=F,xlab="Distance (km)",ylab="Variogram",
           title="Experimental Variogram")
demo.pause("Experimental Variogram")

                                        # Fit a variogram model 
vg.init <- model.create("Nugget Effect",sill=50000,ndim=2)
vg.init <- model.create("Exponential",range=0.15,sill=370000,model=vg.init)
vg.fit  <- model.fit(vg, vg.init, niter=100, wmode=3, draw=T,
                     npairdw=F,xlab="Distance (km)",ylab="Variogram",
                     title="Variogram Model")
demo.pause("Model fitting")

#######
### Global estimation variance 2D
### Estimate = zone mean (over polygon)
#######
                                        # Define the discretization grid
gnx <- 100; gny <- 100;
gd.disc <- db.grid.init(obj=poly.data,nodes=c(gnx,gny))
gd.disc <- db.polygon(gd.disc,poly.data)
plot(gd.disc,pch=3,col=1);plot(db.data,add=T,pch=21,title="Discretization Grid")
plot(poly.data,add=T)
demo.pause("Discretization grid")

                                        # Global estimate = arithmetic mean
cat("Estimating the Global estimate by arithmetic mean...\n")
global.ma <- global(dbin=db.data, dbout=gd.disc, model = vg.fit, uc=c("1"),
                    polygon = poly.data, calcul = "arith", verbose=0)

                                        # Global estimate = kriged mean
cat("Estimating the Global Estimate by Kriging...\n")
global.mk <- global(dbin=db.data, dbout=gd.disc, model = vg.fit, uc=c("1"),
                    polygon = poly.data, calcul = "krige", flag.wgt=TRUE,
                    verbose=0)

                                        # Display kriging weights
db.data <- db.add(db.data,global.mk$wgt,loctype="w")
plot(db.data,name.prop="w",title="Kriging weights")
plot(poly.data,add=T)  
demo.pause("Kriging weights")

                                        # Theoretical process mean
cat("Estimating the Theoretical Mean...\n")
global.mt <- global(dbin=db.data, dbout=gd.disc, model = vg.fit, uc=c("1"),
                    polygon = poly.data, calcul = "mean", flag.wgt=TRUE,
                    verbose=0)

                                        # Summary of results
cat("\nGlobal estimation over the Polygon\n")
tab1 <- rbind(c(global.ma$zest,global.ma$cv),
              c(global.mk$zest,global.mk$cv),
              c(global.mt$zest,global.mt$cv) )
dimnames(tab1) <- list(c("arith mean","kriged zone mean","kriged process mean"),
                       c("Estimate","CV"))
print(round(tab1,3))

#######
###  Testing alternative sampling designs: regular and purely random
#######
                                # Regular grid design : create
x0 <- 25.9; y0 <- 17.0
dx <- 0.06; dy <- 0.06
nx <- 13;   ny <- 12
db.nw <- db.create(x0=c(x0,y0),dx=c(dx,dy),nx=c(nx,ny)) 
db.nw <- db.add(db.nw,loctype="z")
db.nw <- db.locate(db.nw,2:3,loctype="x")
db.nw <- db.add(db.nw,z1=rep(zm,db.nw$nech))
db.nw <- db.locate(db.nw,4,loctype="z");
db.nw <- db.polygon(db.nw,poly.data)

                                # Regular grid design : display
plot(db.nw,pch=3,xlab="km",ylab="km",flag.aspoint=TRUE,name.post=1,
     xlim=c(x1lim,x2lim),ylim=c(y1lim,y2lim),title="Regular Design")
plot(poly.data,add=T,lty=1,density=0)
demo.pause("Display of regular design")

                                # Regular grid design : global estimation
cat("Testing the Regular Design...\n")
global.syst <- global(dbin=db.nw, dbout=gd.disc, model = vg.fit, uc=c("1"),
                      polygon = poly.data, calcul = "arith", verbose=0)

                                # Summary of results
cat("\nTesting alternative sampling designs:\n")
tab2 <- rbind(c(global.ma$zest,global.ma$cv,sum(db.data[,5])),
              c(zm,global.syst$sse/zm,sum(db.nw[,5])),
              c(zm,sqrt(zv/sum(db.data[,5]))/zm,sum(db.data[,5])) )
dimnames(tab2) <- list(c("Data","Regular","Random"),c("Mean","CV","NB"))
print(round(tab2,3))
}