#######
# ICES CRR Handbook of geostatistics in R for fisheries and marine ecology
#
# R 3.0.2
# Rcpp 0.11.0
# RGeostats 11.0.2
# 
# INDICATOR CO-KRIGING USING ACOUSTIC DATA  
#
# The data are supplied by Cheikh-Baye Braham, 
# Institut Mauritanien de Recherche Oceanographique et des Peches (IMROP),
# Nouadhibou, Mauritania
# http://www.imrop.mr/
#
# Reference : Bez, N., and Braham, C.B. 2014. Indicator variables for a robust
# estimation of an acoustic index of abundance. 
# Canadian Journal of Fisheries and Aquatic Sciences, 71, 709-718. 
#
# Author: N. Bez, IRD 
#
#######
{
library(RGeostats)
liste.initiale = ls()
on.exit(demo.clean(liste.initiale))
par(ask=FALSE)
demo.start(FALSE)

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

                                        # Loading data
rg.load("Demo.acoustic.maur.db.data","db.data")
rg.load("Demo.acoustic.maur.poly.data","poly.data")

grid.kri <- db.grid.init(poly.data,margin=10,nodes=150)
grid.kri <- db.polygon(grid.kri,poly.data)
projec.define("mean",db=db.data)

                                        # Display the data
plot(db.sel(db.data,an==1),asp=1,flag.proj=FALSE,title="Data Set (Year #1)")
plot(poly.data,add=T,flag.proj=FALSE)
map("worldHires",add=T,fill=1,col=8)
demo.pause("Displaying Data")

                                        # Define cutoffs and build indicators
zcut <- as.numeric(quantile(db.data [,5][db.data[,5] >0]))
my.limits <- limits.create(zcut=zcut[-5],flag.zcut.int = F)
db.data <- db.indicator(db.data ,my.limits)

                                        # Variograhy
lag     <- c(5,10)
nlag    <- 15
dirvect <- c(0,90)
                                        # Annual variograms
for(i in 1:4){
  plot(vario.calc(db.sel(db.data,an==i),lag=lag,nlag=nlag,dirvect=dirvect),
       flag.norm=T,add=(i!= 1))
}
demo.pause("Experimental variograms")

                                        # Mean annual variogram
vario.data <- vario.calc(db.data,lag=lag,nlag=nlag,dirvect=dirvect,opt.code=1,
                         tolcode=0)

                                        # Model using structures that are linear
                                        # at origin (spherical or exponential)
model.vario <- model.auto(vario.data,struct = c(1,3,3,2,2),wmode=2,npairdw=1)
demo.pause("Modelling the mean annual variogram")

                                        # Co-Kriging

cat("Performing co_Kriging...\n")
neigh.kri <- neigh.create(nmini=10,nmaxi=50,radius=60)
kri.1 <- kriging(db.sel(db.data,an==1),grid.kri,model.vario,neigh.kri)

                                        # Truncating estimations within [0,1]
ranks = db.ident(kri.1,names="Kriging.Indicator*")
for(i in ranks){
  kri.1[,i][kri.1[,i] < 0] <- 0
  kri.1[,i][kri.1[,i] > 1] <- 1
}

                                        # Mapping the results 
plot(kri.1,asp=1,zlim=c(0,1),col=rainbow(4,start=0.2,end=1),flag.proj=FALSE,
     name="Kriging.Indicator.Total.1.estim",title="First Indicator")
legend.image(c(0,1),position="bottomleft",col=rainbow(4,start=0.2,end=1),
             digits=4,cex=0.75)
map("worldHires",add=T)
plot(db.sel(db.data,an==1),col="black",cex1=0.1,add=T,flag.proj=FALSE)
demo.pause("First indicator")

plot(kri.1,asp=1,zlim=c(0,1),col=rainbow(4,start=0.2,end=1),flag.proj=FALSE,
     name="Kriging.Indicator.Total.4.estim",title="Fourth Indicator")
legend.image(c(0,1),position="bottomleft",col=rainbow(4,start=0.2,end=1),
             digits=4,cex=0.75)
map("worldHires",add=T)
plot(db.sel(db.data,an==1),col="black",cex1=0.1,add=T,flag.proj=FALSE)
demo.pause("Fourth indicator")

                                        # Building indicators for intervals 
res <- kri.1
natt.first = res$natt + 1
res <- db.add(res, zero   = (1-res[,ranks[1]]))
res <- db.add(res, low    = (res[,ranks[1]]-res[,ranks[2]]))
res <- db.add(res, medium = (res[,ranks[2]]-res[,ranks[3]]))
res <- db.add(res, large  = (res[,ranks[3]]-res[,ranks[4]]))
res <- db.add(res, extrem = (res[,ranks[4]]))

                                        # Getting rank of most probable interval
res <- db.stat.simu(res,fun="maxi",names=natt.first:res$natt)
res <- db.add(res, class=rep(NA,res@nx[1]*res@nx[2]))
for(i in 1:5)
  res[,"class"][res[,natt.first+i-1]==res[,"maxi"]] <- i

plot(res,asp=1,col=rainbow(5,start=0.2,end=1),flag.proj=FALSE,
     title="Estimated Classes")
legend.image(c(1,5),position="bottomleft",col=rainbow(5,start=0.2,end=1),
             digits=4,cex=0.75)
map("worldHires",add=T)
plot(db.sel(db.data,an==1),col="black",cex1=0.1,add=T,flag.proj=FALSE)
demo.pause("Estimated Classes")
}