Practicals9 - nielshintzen/vmstools GitHub Wiki
In many cases, VMS and Logbook data are used to create maps of fishing effort in specific areas or over time. Although R facilitates making plots, making one that looks good may be hard work. Here we present one potential way of doing it, and you can tailor the specifications to your own needs. See it as a template which you can use over and over, that is at least how it's currently being used.
Setting your picture settings in advance, at the top of your script, might seem a lot of extra work. However, once you've spefied all your settings at once, you'll make sure that in each figure you design afterwards you will use the exact same settings.
library(vmstools)
library(RColorBrewer)
library(sf)
library(ggplot2)
cl <- 1.1 #cex.lab
ca <- 1 #cex.axis
ct <- 1
fonts <- 2 #font
xl <- list(label="Longitude",font=fonts,cex=cl) #x-label
yl <- list(label="Latitude",font=fonts,cex=cl) #y-label
zl <- list(font=fonts,cex=cl) #z-label (if any)
colintens <- brewer.pal(6,"YlOrRd") #colour for fishing intensity
colland <- brewer.pal(9,"PiYG")[8] #colour of land
colgrey <- brewer.pal(9,"Greys") #colour of grey shades
figtype <- "jpeg" #figure extension type
parmar <- rep(2,4) #panel settings
paroma <- (c(6,6,2,2)+0.1) #panel settings
reso <- 1 #dpi setting (1=100dpi)
The next step is setting the spatial extent of the figure, creating a grid to project the VMS data onto, and potentially loading a few additional area sets (such as N2000 areas). The code below shows you how to create such a grid, how to convert your VMS data in a way that it can easily be projected on a grid, and also defines some other areas of interest.
#- Define grid cell area
resx <- 0.25
resy <- 0.125
#-Load the data and convert your data into a 'SpatialPointsDataFrame'
data(tacsat); data(eflalo)
coords <- st_as_sf(tacsat,coords=c("SI_LONG","SI_LATI"))
st_crs(coords) <- 4326
A 'SpatialPointsDataFrame' is the same as a 'SpatialPoints' object, only this one allows you to store per coordinate a lot of data inside a data.frame.
data(ICESareas)
data(europa)
#-Define the area of interest and some reference area
areaInt <- subset(ICESareas,Area_27 %in% c("4.a","4.b","4.c"))
ggplot() + geom_sf(data=st_geometry(areaInt)) + theme_bw()
areaRef <- subset(ICESareas,Area_27 == "4.b")
#-Obtain outer region of my areas and VMS positions. This helps to create maps later on with the same dimensions.
bbox <- cbind(matrix(st_bbox(areaInt),ncol=2),
matrix(st_bbox(areaRef),ncol=2))
rownames(bbox) <- c("x","y")
spatBound <- list(xrange = c(floor(range(bbox["x",])[1]),ceiling(range(bbox["x",])[2])),
yrange = c(floor(range(bbox["y",])[1]),ceiling(range(bbox["y",])[2])))
grd <- createGrid(spatBound$x,spatBound$y,resx,resy,type="GridDF",exactBorder=T)
st_crs(grd) <- 4326
#- If you rather work with a hexagon style grid compared to a rectangular grid, you can use the code below.
bbox <- st_sfc(st_polygon(list(rbind(c(spatBound$x[1],spatBound$y[1]),
c(spatBound$x[2],spatBound$y[1]),
c(spatBound$x[2],spatBound$y[2]),
c(spatBound$x[1],spatBound$y[2]),
c(spatBound$x[1],spatBound$y[1])))))
grd <- st_make_grid(bbox,cellsize=c(resx,resy),square=F)
grd <- st_as_sf(data.frame(grID=1:length(grd)),grd)
st_crs(grd) <- 4326
We have loaded a file called 'europa' before, this is a sf object which allows easy plotting
data(europa)
All areas are defined now. However, once we start plotting the results, quite often we find out that we need some colour scale to present the results at. If for example the data spans from 0 to 1000, some steps like 0-10,10-25,25-50,50-100,100-200,200-500 might look good. Below we define these ranges first, which makes our lives easy when plotting and adding legends to these plots. Linear scales work too, but very often VMS and logbook data are not quite linear in values. Log scales are often more difficult to interpret so that is the reason why we define our own scale here which looks log-scale, but has nice rounded numbers
cutbreaksval <- list(ALL = c(-1,0,0.05,0.1,0.2,0.5,1,500))
legval <- list(ALL = c("0","0 <= 0.05","0.05 <= 0.1", "0.1 <= 0.2","0.2 <= 0.5","0.5 <= 1","> 1"))
#- Potentially, divide your data by a certain constant and add this constant to the legend title
valdiv <- 1000 * 60 #combination of converting from minutes to hours and getting a legend value per 1000 hours
unitval <- c('x 1000 hours per year')
Since we've defined all figure settings and areas, we can start to make some simple plots. Let's start with a figure of the study areas
baseMap <- ggplot() +
geom_sf(data=europa,colour=1,fill=colgrey[5]) +
xlim(spatBound$xrange) +
ylim(spatBound$yrange) +
labs(x=xl$label,y=yl$label) +
theme_bw() +
theme(plot.margin = unit(c(10,10,10,10),"mm"),
axis.title = element_text(face = xl$font,size = rel(xl$cex)))
baseMap +
geom_sf(data=areaInt,colour=NA,fill="yellow") +
geom_sf(data=areaRef,colour=NA,fill="lightblue") +
geom_sf(data=europa,colour=1,fill=colgrey[5])
ggsave(filename=paste("./output/","Map_interestAreas.tiff",sep=""),plot=last_plot(),width = (17/2.54), height = (25/2.54), units = "in", dpi=reso*100)
Note in the picture above that ICES area 4.b was first plotted in yellow and thereafter, on top of that, in blue. Nicely enough, when we use ggsave the symbol for 'degree' is printed correctly in the file while on some Windows machine it shows up as an 'unknown' square on screen.
Since we've now got the study area available, we can add some simple figures on fishing intensity. In principle, the code to plot swept area, value of an area, trawling frequency or intensity is all the same. The only difference will pop up in the column you'll try to aggregate / summarize over.
#- Create a dummy column 'data' in our grid
grd$data <- 0
#-Create column to aggregate over (most often, this column already exists and is output from your previous analyses)
tacsat <- intervalTacsat(tacsat,level="vessel",fill.na=T)
idx <- st_over(coords,grd)
tacsat$gridID <- idx
#- Here we aggregate data to the grid cell. You can aggregate any column you like, we use INTV as an example here.
grd$data[an(names(table(idx)))] <- aggregate(tacsat$INTV,by=list(tacsat$gridID),FUN=sum,na.rm=T)$x
We have specified the binns we use for the legend and colours above, but it might be handy if you see where I got the range in values from. Below, you can see the 'grd@data[an(names(table(idx))),1]' piece of code. This selects only those coordinates from the data.frame that are associated with a value > 0. All in all, this is a simple vector and allows me to see what values it spans.
#Look at value ranges:
range(grd$data[an(names(table(idx)))])
#-Here we cut the range of values up into colour binns
cols <- c(NA,colintens)[cut(grd$data[an(names(table(idx)))]/valdiv,breaks=cutbreaksval$ALL)]
baseMap +
geom_sf(data=grd[an(names(table(idx))),],colour=NA,aes(fill=cols)) +
geom_sf(data=europa,colour=1,fill=colgrey[5]) +
#-Add a legend
guides(fill=guide_legend(title=unitval)) +
scale_fill_manual(values = rev(colintens),labels=rev(legval$ALL))
#- Save the plot
ggsave(filename=paste("./output/","Map_griddedAreas.tiff",sep=""),plot=last_plot(),width = (17/2.54), height = (25/2.54), units = "in", dpi=reso*100)
If you prefer not to plot it straight into R, but rather export it to e.g. ArcGIS, then below you'll find a routine to export the 'grdPols' we created above to a shapefile.
grd$colour <- "white"
grd$colour[an(names(table(idx)))] <- cols
st_write(grd,paste0("./output/grid.shp"))
- Try to open the shapefile in your GIS package and investigate if it looks the same.
Creating tables out of the VMS and Logbook data is often based on 'point in polygon' types of calculations where GPS positions inside a reference area or area of interest are aggregated. Here we show how you could also merge logbook records with VMS records in the same type of analyses.
Let's start with calculating the time spent in areas IV and IVb in specific by VMS and Logbook data.
data(tacsat); data(eflalo)
#Calculate effort associated with each VMS ping
tacsat <- intervalTacsat(tacsat,level="vessel",fill.na=T)
tacsat <- tacsat[which(is.na(tacsat$SI_LONG)==F & is.na(tacsat$SI_LATI)==F),]
#Calculate the effort in an eflalo dataset
eflalo <- effort(eflalo,unit="mins",by='byRow')
#Because some trip numbers are replicated, we take only unique trip numbers
# (a single trip can be associated with multiple eflalo records, that means that also time out of harbour and in harbour are replicated)
eflalo <- eflalo[!duplicated(eflalo$FT_REF),]
#Assign GPS positions to eflalo
eflalo$SI_LONG <- ICESrectangle2LonLat(eflalo$LE_RECT)$SI_LONG
eflalo$SI_LATI <- ICESrectangle2LonLat(eflalo$LE_RECT)$SI_LATI
eflalo <- eflalo[which(is.na(eflalo$SI_LONG)==F & is.na(eflalo$SI_LATI) == F),]
coordsTacsat <- st_as_sf(tacsat,coords=c("SI_LONG","SI_LATI"))
coordsEflalo <- st_as_sf(eflalo,coords=c("SI_LONG","SI_LATI"))
st_crs(coordsTacsat) <- 4326
st_crs(coordsEflalo) <- 4326
#Determine tacsat GPS position inside area of interest (IVa, IVb and IVc)
# and determine which midpoint of eflalo rectangles are within the area of interest
inTacsat <- st_over(coordsTacsat,areaInt)
inEflalo <- st_over(coordsEflalo,areaInt)
#Calculate here the total minutes spend by vessels in these areas
sum(tacsat$INTV[ which(is.na(inTacsat)==F)],na.rm=T) +
sum(eflalo$EFFORT[which(is.na(inEflalo)==F)],na.rm=T)
Easier would it be if we could combine tacsat and eflalo straight away. This especially pays of if you start to divide ICES rectangles up into e.g. 1/4th 1/16th or even small areas. The midpoints of each of these 1/4th or 1/16th ICES rectangles can be treated similarly.
eflaloTrim <- eflalo[,c("SI_LONG","SI_LATI","EFFORT")]
colnames(eflaloTrim)[3] <- "INTV"
comb <- rbind(tacsat[,c("SI_LONG","SI_LATI","INTV")],
eflaloTrim)
coordsComb <- st_as_sf(comb,coords=c("SI_LONG","SI_LATI"))
st_crs(coordsComb) <- 4326
inComb <- st_over(coordsComb,areaInt)
#Calculate the total minutes spend by vessels from tacsat and eflalo combined
sum(comb$INTV[which(is.na(inComb)==F)],na.rm=T)
#Now do it for IVb only
inComb <- st_over(coordsComb,areaRef)
sum(comb$INTV[which(is.na(inComb)==F)],na.rm=T)
#Or we define an arbitrary polygon ourselves
arbPol <- st_as_sfc(list(lonLat2SFPolygons(SI_LONG=c(2,3,3.5,2),SI_LATI=c(53,54,56.4,53))))
st_crs(arbPol) <- 4326
baseMap + geom_sf(data=arbPol,colour=NA,fill="1")
inArbPol <- st_over(coordsComb,arbPol)
sum(comb$INTV[which(is.na(inArbPol)==F)],na.rm=T)
R in its basic form does not have an interactive graphic capability, allowing you to zoom in and out, like most bona fide GIS systems. It can be done, however, using the library leaflet. A student (ask for credential with the course instructors) from a previous course was kind enough to share her take on a leaflet map, please find it below.
data(tacsat);data(eflalo)
tacsatp <- mergeEflalo2Tacsat(eflalo,tacsat)
tacsatp$IDX <- 1:nrow(tacsatp)
#A quick and dirty assessment of fishing activity: speeds between 1 and 6 are fishing
tacsatFilter <- filterTacsat(tacsatp,st=c(1,6),hd=NULL,remDup=T)
tacsatp$SI_STATE <- 0
tacsatp$SI_STATE[tacsatFilter$IDX] <- 1
#I needed the activity distinction because I only want to dispatch catches or values over fishing pings!
tacsatEflalo <- splitAmongPings(tacsat=tacsatp,eflalo=eflalo,
variable="all",level="day",conserve=T)
#Add a grid notation (csquaraes)
tacsatEflalo$LE_SPAT <- CSquare(tacsatEflalo$SI_LONG,tacsatEflalo$SI_LATI,degrees=0.1)
aggregate(tacsatEflalo$LE_KG_PLE,by=list(tacsatEflalo$LE_SPAT),FUN=sum,na.rm=T)
#In a quick picture, it looks like this
x11()
# check this overlay problem!
plotTools(tacsatEflalo,level="gridcell",xlim=c(-23,23),ylim=c(40,64),zlim=NULL,log=F,gridcell=c(0.1,0.1),color=NULL,control.tacsat=list(clm="LE_KG_PLE"))
# In leaflet it looks like this
library(leaflet)
library(tidyverse)
tacsatEflalo %>%
select(LE_SPAT, LE_KG_PLE) %>%
group_by(LE_SPAT) %>%
summarise(PLE=sum(LE_KG_PLE, na.rm=TRUE)) %>%
filter(!is.na(PLE) & PLE>0) -> tacsatToMap
# Add C-squares
tacsatToMap$MID_LATI=(CSquare2LonLat(tacsatToMap$"LE_SPAT", 0.1))$SI_LATI
tacsatToMap$MID_LONG=(CSquare2LonLat(tacsatToMap$"LE_SPAT", 0.1))$SI_LONG
tacsatToMap$dw_LATI=tacsatToMap$MID_LATI-0.05
tacsatToMap$dw_LONG=tacsatToMap$MID_LONG-0.05
tacsatToMap$gw_LATI=tacsatToMap$MID_LATI+0.05
tacsatToMap$gw_LONG=tacsatToMap$MID_LONG+0.05
colorNumeric(palette = heat.colors(10), domain = tacsatToMap$PLE, reverse=TRUE) -> pal #
colorNumeric(palette = heat.colors(10), domain = tacsatToMap$PLE, reverse=TRUE) -> pal
addTiles = function (map, urlTemplate = "http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png",
attribution = NULL, layerId = NULL, group = NULL, options = tileOptions())
{
options$attribution = attribution
if (missing(urlTemplate) && is.null(options$attribution))
options$attribution = paste("© OpenStreetMap",
"contributors, CC-BY-SA")
invokeMethod(map, getMapData(map), "addTiles", urlTemplate,
layerId, group, options)
}
leaflet(tacsatToMap) %>%
addTiles() %>%
fitBounds(~min(dw_LONG), ~min(dw_LATI), ~max(dw_LONG), ~max(dw_LATI)) %>%
setView(lng = 2, lat = 53, zoom = 7) %>%
addRectangles(lng1=~dw_LONG, lat1=~dw_LATI, lng2=~gw_LONG, lat2=~gw_LATI,fill=TRUE,color = ~pal(PLE),stroke=FALSE, fillOpacity=0.5 )