###########
# 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 
#
# Author:       Mathieu Woillez, Ifremer
###########
{
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)

                                        # Load data
rg.load("Demo.herring.sa.scot.db.data","db.data")
rg.load("Demo.herring.sa.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="sa")
hist(db.extract(db.data,"sa"),col=8,xlab="sa",nclass=100,main="Histogram of Initial Data")
demo.pause("Histogram")

                                        # Visualization 
plot(db.data,inches=5,asp=1/cos(mean(db.extract(db.data,"x2"))*pi/180),
     pos.legend=5,zmin=0,zmax=c(db.stat(db.data,"maxi")),include.bounds=FALSE)
plot(poly.data,add=T)
map("worldHires",add=T,fill=T,col=8)
demo.pause("Data location")

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

                                        # Compute areas of influence (weights)
cat("Calculating the Areas of Influence...\n")
db.data <- infl(db.data,nodes=c(300,300),origin=c(-4,57.5),
                extend=c(6.5,4.5),radius=12.5,
                polygon=poly.data,plot=T,asp=1,title="Areas of Influence")
demo.pause("Areas of Influence")

                                        # Define the anamorphosis model
model.anam <- anam.fit(db.data,type="emp",ndisc=db.data$nactive,sigma2e=800,
                       draw=T,title="Gaussian Anamorphosis")
demo.pause("Gaussian Anamorphosis")

                                        # Transform the data into Gaussian
db.data <- anam.z2y(db.data,anam=model.anam)
db.data <- db.rename(db.data,name="Gaussian.sa",newname="Yp")

print(db.data,flag.stats=TRUE,names="Yp")
#ycut <- qnorm(sum(db.extract(db.data,"sa") == 0) / db.data$nactive)
Y <- db.extract(db.data,"Yp")
# To avoid numerical uncertainty, the Gaussian value is considered as the
# minimum of Y, rather than the theoretical value derived from proportion of 0
ycut = range(Y)[1] 
Y[Y <= ycut] <- ycut
db.data <- db.replace(db.data,"Yp",Y)
print(db.data,flag.stats=TRUE,names="Yp")
hist(db.extract(db.data,"Yp"),col=8,xlab="Yp",nclass=100,main="Histogram of Gaussian transformed data")
demo.pause("Histogram of Yp")

# We calculate the experimental of the gaussian transformed data 
# (including the 0 trasnformed values)
# We derive the experimental of the underlying Gaussian variable
# where the zero-data have been expanded to negative possible values

n.H <- 50
vario.Yp <- vario.calc(db.data,lag=2.5,nlag=50)
vario.Y  <- vario.trans.cut(vario.Yp,ycut,n.H)
model.vario.Y <- model.auto(vario.Y,struc=melem.name(c(1,2,3)),draw=F)
plot(vario.Y ,npairdw=T,npairpt=F,inches=0.08,col="red",
     title="Gaussian Model")
plot(vario.Yp,npairdw=T,npairpt=F,inches=0.08,col="black",add=T)
plot(model.vario.Y,add=T,col="red")
demo.pause("Gaussian Model")


                                        # Define the moving neighborhood
neigh.simu <- neigh.create(ndim=2,type=2,nmini=5,nmaxi=100,radius=60)
neigh.simu

                                        # Define interval limits for the gibbs
Ymax <- db.extract(db.data,name="Yp",flag.compress=F)
Ymin <- db.extract(db.data,name="Yp",flag.compress=F)
Ymin[Ymin <= ycut] <- -10
db.data<-db.add(db.data,Ymax)
db.data<-db.locate(db.data,db.data$natt,"upper")
db.data<-db.add(db.data,Ymin)
db.data<-db.locate(db.data,db.data$natt,"lower")

## Simulating gaussian values below ycut at datapoints 
## while honouring the gaussian variable model and conditional on the other
## gaussian data values:

## A Gibbs sampler simulates a gaussian value at each point between its
## boundaries ymin and ymax:
## within [-10; ycut] where raw data is 0;
## and ymin = ymax = y where y is known

cat("Performing the Gibbs sampler (100 iterations)...\n")
db.data <- gibbs(db = db.data, model = model.vario.Y, seed = 232132, 
                 nboot = 10, niter = 100, flag.norm=FALSE, percent=0,
								 radix = "Gibbs", modify.target = TRUE)
db.data<-db.rename(db.data,"Gibbs.G1","Y")
print(db.data,flag.stats=TRUE,names="Y")

                                        # Visualization Gibbs sampling step
histYp <- hist(db.extract(db.data,name="Yp"),plot=F,breaks=seq(-4,4,.1))
hist(db.extract(db.data,name="Yp"),proba=T,breaks=seq(-4,4,.1),
     xlab="Y+",col=8,main="",
     xlim=c(-4,4),ylim=c(0,ceiling(max(histYp$density))))
demo.pause("Histogram")

hist(db.extract(db.data,name="Y"),proba=T,breaks=seq(-4,4,.1),
     xlab="Y",main="",col=8,
     xlim=c(-4,4),ylim=c(0,ceiling(max(histYp$density))))
lines(seq(-4,4,0.1),dnorm(seq(-4,4,0.1),0,1),col=2)
demo.pause("Completed Histogram")

# Overlay the experimental variogram calculated on expanded data
# with the model of underlying expanded information
plot(vario.calc(db.data,lag=2.5,nlag=50),npairdw=T,npairpt=F,
     ylab=expression(gamma),main="",xlab="Distance (nmi)",inches=0.08)
plot(model.vario.Y,add=T,col=2)
demo.pause("Gaussian Model")

## For each simulation, gaussian values are now defined at all datapoints
## They will be used for classical gaussian simulation

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

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

cat("Performing Conditional Simulations...\n")
grid.simu <- simtub(dbin=db.data, dbout=grid.simu, model=model.vario.Y,
                    neigh=neigh.simu, uc = "", mean = 0, seed = 232132,
                    nbsimu = 1, nbtuba = 200,
                    radix = "Simu", modify.target = TRUE)
grid.simu <- db.rename(grid.simu,"Simu.Y.S1","Simu.Y")
print(grid.simu,flag.stats=TRUE,names="Simu.Y")
hist(db.extract(grid.simu,"Simu.Y"),col=8,xlab="Y",nclass=100,main="Histogram of Simulated Gaussian expanded values")

                                        # Transform gaussian cond simulation
                                        # into raw cond simulation
grid.simu <- anam.y2z(grid.simu,name="Simu.Y",anam=model.anam)
print(grid.simu,flag.stats=TRUE,names="Raw.Simu.Y")
hist(db.extract(grid.simu,"Raw.Simu.Y"),col=8,xlab="Raw",nclass=100,main="Histogram of Simulated Raw values")


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

                                        # Display conditional simulation of Z 
pal2 <- colorRampPalette(c("cyan", "yellow", "red", "black"), bias=4)
plot(poly.data,col=4,asp=1/cos(mean(db.extract(db.data,"x2"))*pi/180),
     flag.proj=F,title="Conditional Simulation (Raw)")
plot(grid.simu,name="Raw.Simu.Y",col=pal2(100),pos.legend=5,flag.proj=F,add=T)
map("worldHires",add=T,fill=T,col=8);box()
demo.pause("Simulation (Raw)")

                                        # Display conditional simulation of Z>0
Raw.Simu.Y.sup0 <- grid.simu@items$Raw.Simu.Y
Raw.Simu.Y.sup0[Raw.Simu.Y.sup0 <= 0.00] <- NA
grid.simu <- db.add(grid.simu,Raw.Simu.Y.sup0)
plot(poly.data,col=4,asp=1/cos(mean(db.extract(db.data,"x2"))*pi/180),
     flag.proj=F,title="Conditional Simulation (positive)")
plot(grid.simu,name="Raw.Simu.Y.sup0",col=pal2(100),pos.legend=5,
     flag.proj=F,add=T)
map("worldHires",add=T,fill=T,col=8);box()
demo.pause("Cond Simulation Z>0")
}