####### ## ICES CRR Handbook of Geostatistics for Fisheries and Marine Ecology ## ## This code performs global estimation in 1D using acoustic transect sums ## The approach applies to acoustic surveys with regularly spaced parallel ## transects ## The code uses the data files where fish biomass per transect ## was obtained by summing the densities along the transects ## ## Author: P.Petitgas, Ifremer ####### { library(RGeostats) liste.initiale = ls() on.exit(demo.clean(liste.initiale)) par(ask=FALSE) demo.start(FALSE) ####### ## 1D data (transect sums) ####### # Load data file rg.load("Demo.anchovy.bob.1d.db.data","db.data") nrad <-db.data$nsamples # Nb of transects aa <- 1 # Inter-transect (arbitrary) distance # Transform db.data into regular grid db.datagrid <- db.grid.init(db.data,nodes=nrad,flag.regular=T) db.datagrid <- migrate(db.data,db.datagrid,flag.fill=2,name="Tr.biomass") # Display information plot(db.data,pch=20,type="b",title="Biomass", xlab="S <---- Transects ----> N", ylab="Biomass per transect") plot(db.datagrid,add=TRUE,col="red") demo.pause("Display information") ####### ## Estimation globale 1d : intrinsic method + geometric error ####### # Experimental variogram 1D vg <- vario.calc(db.data,calcul="vg") # Semi-automatic fit (sills only) vg.init <- model.create(vartype="Nugget Effect",sill=1.5E+05,ndim=1) vg.init <- model.create(vartype="Spherical",range=4.5, model=vg.init) vg.fit <- model.fit(vg, vg.init, niter=100, wmode=3,draw=FALSE) # Automatic fit vg.auto <- model.auto(vg,struc=c("Nugget Effect","Spherical"), xlab="Distance", ylab="variogram",draw=FALSE) # Variogram and model representations plot(vg,npairpt=0,npairdw=T,xlab="Distance", ylab="variogram",inches=.025, title="1-D Variogram") plot(vg.fit ,add=TRUE,col="blue") plot(vg.auto,add=TRUE,col="red") demo.pause("1-D variogram") # Choose model vgmod <- vg.auto perc <- vgmod[1]$sill/(vgmod[1]$sill + vgmod[2]$sill) cat("Percent of nugget in total sill =",perc,"\n") # Estimation variance 1D gloa <- global(db.datagrid,calcul="arith",model=vgmod,ndisc=100,verbose=0) s2 <- var(db.data[,"Tr.biomass"],na.rm=T)*(nrad-1)/nrad d2geom <- s2*(aa^2/6)/(aa*nrad)^2 cv.geom <- sqrt(d2geom)/gloa$zest # Geometric error (limits of 1D) cv.tot <- sqrt(gloa$sse^2+d2geom)/gloa$zest # Total error CV cat("Mean=",round(gloa$zest,3), "CVest=",round(gloa$sse/gloa$zest,3),"\n") cat("CV.est=",round(gloa$cv,3)," CV.geom=",round(cv.geom,3),"\n") cat("Qtot=",gloa$zest*nrad*aa," CV=",round(cv.tot,3),"\n") ####### ## Alternative sampling effort: ## Estimation variance for other ('nk') inter-transect distances ## (geometric error neglected) ####### nk <- 9 # Loop on alternative inter-Tr distances sse <- numeric(nk) for (k in 1:nk) { ak <- k*0.25*aa # new inter-transect distance nrk <- round(nrad*aa/ak,0) # new nb of transects d2geom <- s2*(ak^2/6)/(ak*nrk)^2 # variance geometric error # variance estimation error db.datagrid <- db.grid.init(db.data,nodes=nrk,flag.regular=T) db.datagrid <- migrate(db.data,db.datagrid,flag.fill=2,name="Tr.biomass") d2estim <- global(db.datagrid,calcul="arith",model=vgmod,ndisc=100, verbose=0)$sse^2 sse[k] <- sqrt(d2estim+d2geom) } plot(0.25*aa*(1:nk),sse/gloa$zest,type="b", xlab="Multiplier of Inter-transect Distance",ylab="Estimation CV") demo.pause("Alternative Inter-transect distances") }