#######
# ICES CRR Handbook of geostatistics in R for fisheries and marine ecology
#
# GLOBAL AND LOCAL ESTIMATIONS USING TRANSITIVE GEOSTATISTICS 
#
# The data are supplied by Abdelmalek Faraj
# Institut National de Recherche Halieutique (http://www.inrh.ma/),
# Casablanca, Morocco. 
#
# Reference : Faraj, A., and Bez, N. 2007. Spatial considerations for the
# Dakhla stock of Octopus vulgaris: indicators, patterns and fisheries
# interactions. 
# ICES Journal of Marine Science, 64, 1820-1828.
#
# 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)

                                        # Load data
rg.load("Demo.octopus.morocco.db.data","db.data")
rg.load("Demo.octopus.morocco.poly.data","poly.data")
projec.toggle(0,verbose=FALSE)
db.data = db.locate(db.data,c("long","lat"),"x")  

strata = 11                             # Dimension of a regular strata in nm

                                        # Defining a projection
projec.define("mean",db=db.data)

                                        # Displaying data
plot(db.data,zmin=0.1,pch.low=3,cex.low=0.25,las=1,pch=21,col=1,inches=2.5,
     title="Juveniles of octopus 1999 (nb/nm2)",asp=0.8,flag.proj=FALSE)
plot(poly.data,las=1,add=T,flag.proj=FALSE)
map("worldHires", interior=F, add=T, fill=1, col=8)
demo.pause("Display Data")

                                        # Surface of influence for inner data
                                        # in square nautical miles 
db.data <- infl(db.data,nodes=400,extend=c(6,6),origin=c(-18,20.5),
                polygon=poly.data,plot=T,asp=0.8,
                title = "Influence Polygons")
demo.pause("Influence Polygons")

                                        # Covariogram computation
                                        # Lag = Regular strata 
lag     <- strata
nlag    <- 20
dirvect <- 50 + c(0,90)
covario.data <- vario.calc(db.data,lag=lag,nlag=nlag,calcul="covg",tolang=45,
                           dirvect=dirvect)

                                        # Calculate relative g(h) scaled by
                                        # Q (total abundance)
Q <- sum(db.data[,"JUV"]*db.data[,"Influence.Surface"])
covario.data = vario.transfo("v1/Q^2",covario.data)
variance = covario.data$vars

                                        # Remove the smallest lags 
covario.data@vardirs[[1]]@sw[c(nlag,nlag+2)] <- 0
covario.data@vardirs[[2]]@sw[c(nlag,nlag+2)] <- 0

                                        # Ajust the empirical covariogram with
                                        # variance constraints
model.covario <- model.auto(covario.data,struct=c(1,2),constraints=variance,
                            npairdw=1,las=1,inches=0.05,lwd=2,xlim=c(0,250),
                            title = "Relative Geometrical Covariogram")
demo.pause("Relative Geometric Covariogram")

                                        # Estimation variance 
projec.toggle(0,verbose=FALSE)  
CVtrue <- sqrt(strata^2*(variance -
                         model.cvv(v.mesh=strata,model=model.covario,
                                   seed=110366,ndisc=20)))

cat(paste("Abundance estimate = ",round(Q/10^6,0)," e+06",sep=""),"\n")
cat(paste("Estimation Coefficient of Variation = ", round(100*CVtrue,1),"%",
          sep=""),"\n")

                                        # Grid definition in geographical space
projec.toggle(0,verbose=FALSE) 
grid.kri <- db.grid.init(db.data,nodes=400,extend=c(6,6),origin=c(-18,20.5))
grid.kri <- db.polygon(grid.kri,poly.data)
projec.toggle(1,verbose=FALSE) 

                                        # Neighborhood definition
neigh.kri <- neigh.create(ndim=2,type=2,flag.aniso=TRUE,flag.rotation=TRUE,
                          nmini=10,nmaxi=30,radius=c(150,50),
                          rotmat=util.ang2mat(ndim=2,angles=50))

                                        # Transitive kriging
res <- kriging(db.data,grid.kri,model.covario,neigh.kri)

                                        # Threshold negative estimates
res[,"Kriging.JUV.estim"][res[,"Kriging.JUV.estim"]<0] <- 0

                                        # Map of the result
plot(res,col=rainbow(6,start=0.2,end=1),las=1,
     title="Octopus density - 1999",asp=0.8,flag.proj=FALSE)
plot(poly.data,las=1,add=T,flag.proj=FALSE)
map("worldHires", fill=T,col=grey(0.8),add=T)
plot(db.data,las=1,add=T,col=1,flag.proj=FALSE)
legend.image(range(db.extract(res,"Kriging.JUV.estim"),na.rm=T),
             position="bottomright",col=rainbow(6,start=0.2,end=1),
             ntdec=0,cex=0.75)
demo.pause("Displaying Results")
}