#### 
# ICES Cooperative Research Report
# Handbook of Geostatistics in R for fisheries and marine ecology
#
# Based on 
# ICES Training Course
# Application of Geostatistics to analyse spatially explicit Survey data
# in an Ecosystem Approach
# December 2013-2014, Fontainebleau, Centre de Geosciences, Mines ParisTech 
#
# Demo script:  Conditional simulations with the presence of zeroes         
#### 
{
library(RGeostats)
liste.initiale = ls()
on.exit(demo.clean(liste.initiale))
par(ask=FALSE)
demo.start(FALSE)

                                        # Load packages
library(mapdata)
library(maps)
projec.toggle(0,verbose=FALSE)

                                        # Loading data
rg.load("Demo.herring.len.scot.db.data","db.data")
rg.load("Demo.herring.len.scot.poly.data","poly.data")

                                        # Select points inside polygon
db.data <- db.polygon(db.data,poly.data,verbose=TRUE)

                                        # Print statistics and histogram
print(db.data,flag.stats=TRUE,names="m.length")
hist(db.extract(db.data,"m.length"),col=8,xlab="m.length",main="",nclass=20)
demo.pause("Histogram")

                                        # Define the projection 
projec.define(projection="mean", db=db.data)

                                        # Visualization 
plot(db.data,inches=5,asp=1/cos(mean(db.extract(db.data,"x2"))*pi/180),
     pos.legend=5,zmax=c(db.stat(db.data,"maxi")),include.bounds=FALSE,
     flag.proj=FALSE,title="Fish Length")
plot(poly.data,col=4,add=T,flag.proj=FALSE)
map("worldHires",add=T,fill=T,col=8)
demo.pause("Display Information")

                                        # Define the anamorphosis model
model.anam <- anam.fit(db.data,type="gaus",nbpoly=10,draw=T,
                       title = "Gaussian Anamorphosis")
demo.pause("Gaussian Anamorphosis")

                                        # Transform the data into Gaussian
db.data <- anam.z2y(db.data,anam=model.anam)

                                        # Modeling Gaussian variable Y
vario.data <- vario.calc(db.data)
model.vario <- model.auto(vario.data,struc=melem.name(c(1,3,12)),wmode=2,
                          npairdw=T,npairpt=F,inches=0.08,col="black",
                          title="Model of Gaussian Transformed")
demo.pause("Model of Gaussian Y")

                                        # Define simulation grid 
gnx <- 144 
gny <- 90
grid.simu <- db.grid.init(poly.data,nodes=c(gnx,gny))
grid.simu <- db.polygon(grid.simu,poly.data)
grid.simu

                                        # Display of the grid
plot(grid.simu,col=1,title="",pch="+",asp=1)
plot(db.data,inches=3,col=2,pch=19,add=T) 
plot(poly.data,col=4,add=T)
demo.pause("Display Grid")

                                        # Define the neighborhood
neigh.simu <- neigh.create(ndim=2,type=0) 
neigh.simu

                                        # Cond. simulation of gaussian variable
grid.simu <- simtub(dbin=db.data, dbout=grid.simu, model=model.vario,
                    neigh=neigh.simu, uc = "", mean = 0, seed = 29091978,
                    nbsimu = 1, nbtuba = 1000, radix = "Simu",
                    modify.target = TRUE)
print(grid.simu,flag.stats=TRUE,names="Simu.Gaussian.m.length.S1")

                                        # Transform gaussian cond simulation
                                        # into raw cond simulation
grid.simu <- anam.y2z(grid.simu,name="Simu.Gaussian.m.length.S1",
                      anam=model.anam)
print(grid.simu,flag.stats=TRUE,names="Raw.Simu.Gaussian.m.length.S1")

                                        # Display cond simulation of gaussian
plot(grid.simu,name="Simu.Gaussian.m.length.S1",
     asp=1/cos(mean(db.extract(db=db.data,names="x2"))*pi/180),
     title="Conditional Simulation of Gaussian Variable",
     pos.legend=5,flag.proj=F)
plot(poly.data,col=0,flag.proj=F,add=T)
map("worldHires",add=T,fill=T,col=8);box()
demo.pause("Display Cond. Simulation (Gaussian)")

                                        # Display conditional simulation of Z 
plot(grid.simu,name="Raw.Simu.Gaussian.m.length.S1",
     asp=1/cos(mean(db.extract(db=db.data,names="x2"))*pi/180),
     title="Conditional Simulation of Raw Variable",
     pos.legend=5,flag.proj=F)
plot(poly.data,col=0,flag.proj=F,add=T)
map("worldHires",add=T,fill=T,col=8);box()
demo.pause("Display Cond. Simulation (Raw)")
}