DIY indicators

Calculating indicators tailor-made: Optional

If you still have time and are particularly interested in those indicators you can calculate them yourself. For a start, let's take a simple indicator to calculate ourselves. Let's define: the total surface area trawled by year. What we need to do for this exercise is:

1. clean the tacsat and eflalo data
2. merge them together
3. perform activity analyses and take fishing activity only
4. select bottom trawling gear only
5. define the width of the gear
6. define the grid

We have the choice to interpolate the data or not, but to keep things simple, we will not do that here. Obviously, the function to add a width to a gear, seen before, could be useful. However, we take less of a true spatial approach here (for most spatial analyses, your data structure is fine as vectors rather than as a grid!), whereby the 'addWidth' function would simply take too long to run.

data(tacsat);data(eflalo)

#Now clean the tacsat and eflalo dataset...
tacsatp <- mergeEflalo2Tacsat(eflalo,tacsat)
tacsatp$LE_GEAR <- eflalo$LE_GEAR[match(tacsatp$FT_REF,eflalo$FT_REF)]

#Again, the quick and dirty fishing activity analyses
tacsatp <- filterTacsat(tacsatp,st=c(2,6),hd=NULL,remDup=T)

#Select gears
table(tacsatp$LE_GEAR)
gears <- c("OTB","PTB","TBB")
tacsatp <- subset(tacsatp,LE_GEAR%in% gears)
gearWidth <- data.frame(LE_GEAR=gears,LE_WIDTH=c(87,2*87,24))

#Define the grid & subset the area
spatBound <- list(x=c(-5,10),y=c(48,62))
idxlon    <- which(tacsatp$SI_LONG >= spatBound$x[1] & tacsatp$SI_LONG <= spatBound$x[2])
idxlat    <- which(tacsatp$SI_LATI >= spatBound$y[1] & tacsatp$SI_LATI <= spatBound$y[2])
tacsatp   <- tacsatp[idxlon[which(idxlon %in% idxlat)],]
grd <- createGrid(xrange=c(unlist(spatBound$x)),
yrange=c(unlist(spatBound$y)),
resx=0.1,resy=0.05) #Let's take a tenth of an ICES rectangle for a start

So far so good, you've all seen these preparation steps a number of times now, but it's good to learn them by heart as you will need to go over these time and time again (look at Practicals2).

The next step is to start calculating the area trawled for each VMS ping. It might be handy here to store the results by year in some sort of data-frame again, so you can access it later, and also derive other results from it.

Area trawled = duration of trawling x width of gear x speed of a vessel. Let's work this equation out for each of the VMS pings in the dataset.

tacsatp <- sortTacsat(tacsatp)

#Calculate duration of trawling per ping
tacsatp <- intervalTacsat(tacsatp,level="trip",fill.na=TRUE)
tacsatp$SI_INTV <- tacsatp$INTV
gearWidth$LE_GEAR <- ac(gearWidth$LE_GEAR)
tacsatp <- merge(tacsatp,gearWidth,by="LE_GEAR")
tacsatp$TR_AREA <- (tacsatp$SI_INTV / 60) * (tacsatp$LE_WIDTH / 1000) * (tacsatp$SI_SP *1.852)

For the full tacsat dataset we have now calculated the area trawled per ping, however, not yet per year. We could easily split up the tacsat dataset into the two years (1800 and 1801) and aggregate over the trawled area column. But, as indicated before, if we want a bit more flexibility, it might be good to store the data first and perform calculations later-on.

#Turn the defined grid into a spatial dataframe
sg <- SpatialGrid(grid=grd)
spgDF <- as(sg,"SpatialGridDataFrame")
slotNames(spgDF)

#We can simply add a dataframe to this object (but make sure it has the same size as all coordinates considered!
spgDF@data <- data.frame("Y1800"=rep(0,nrow(coordinates(spgDF))),
"Y1801"=rep(0,nrow(coordinates(spgDF))),
surface=rep(0,nrow(coordinates(spgDF))))
head(spgDF@data)

#Now align the VMS pings to the grid
#We have to turn the coordinates of the tacsat dataset into spatial points first
for(iYr in c(1800:1801)){
subTacsat <- subset(tacsatp,year(SI_DATIM)==iYr)
sp <- SpatialPoints(coordinates(subTacsat[,c("SI_LONG","SI_LATI")]))

#The spatial grid dataframe is a long list (nrow(spgDF@data)) and the idx tells me which row number (associated with a grid cell) the VMS ping is located in.
idx <- over(sp,sg)
spgDF@data[ac(aggregate(subTacsat$TR_AREA,          by=list(idx),FUN=sum,na.rm=T)[,1]), paste("Y",iYr,sep="")]    <- aggregate(subTacsat$TR_AREA,  by=list(idx),FUN=sum,na.rm=T)[,2]
}

#Have a look at some data (where Y1800 is not zero)
head(spgDF@data[which(spgDF@data$Y1800>0),])

#And calculate the indicator value
colSums(spgDF@data[,c("Y1800","Y1801")])

However, we just set this system up to be a bit more flexible and use the data for other questions too. Can we now easily move from total area trawled to percentage of the area untrawled? Or could we even see which area has been trawled in 1800 but not in 1801?

#Add surface of each gridcell to the dataframe
spgDF@data$surface <- surface(spgDF)$cellArea

#Percentage trawled
idxun <- which(spgDF@data[,"Y1800"] == 0)
sum(spgDF@data$surface[idxun]) / sum(spgDF@data$surface)

#Area trawled in 1800 but not in 1801
idx <- which(spgDF@data[,"Y1800"] > 0 & spgDF@data[,"Y1801"] == 0)
sum(spgDF@data$surface[idx]) / sum(spgDF@data$surface)

All in all, a lot to play with in a rather flexible way.