Background

This code accompanies the publication “Making the most of data-poor fisheries: low cost mapping of small island fisheries to inform policy” (in review) Marine Policy. It can be used to generate maps of relative fishing intensity such as those in figures 2-4 in the manuscript.

Load required packages and data

The following r packages are required to create the maps:

#Required packages
library(raster)
library(rgdal)
library(tmap)
library(colorRamps)
library(rasterVis)
library(grid)

Import the fishing area shapefile, along with other shapefiles to include in the plots (e.g. no-fishing areas, major landmarks, land outline). The attribute data for each fishing area polygon includes:

# Here **dsn** is the folder location of the shapefiles
FGArea <- readOGR(dsn,"FishingGrounds 1Apr17")      # fishing area polygons with attribute data
# Other shapefiles 
Nofishing <- readOGR(dsn,"NoFishingAreasAll_erase") # no fishing areas
towns <- readOGR(dsn,"MajorTowns")                  # major towns and communities
land <- readOGR(dsn,"LandPolyine")                  # polygon with outline of the land area
plot(FGArea)

Preparing the spatial data

Ensure that all the data have the same projected coordinate system. For Barbados, we use WGS84 UTM 21N, and for convenience, we set our units to km.

proj1 <- CRS('+proj=utm +zone=21 +datum=WGS84 +units=km +no_defs +ellps=WGS84 +towgs84=0,0,0')
FGArea <- spTransform(FGArea,CRS=proj1)
Nofishing <- spTransform(Nofishing,CRS=proj1)
towns <- spTransform(towns,CRS=proj1)
land <- spTransform(land,CRS=proj1)

Plot and inspect the shapefiles to ensure that they line up accurately

plot(FGArea);plot(Nofishing, add=T, col='black');plot(land,add=T,col='light gray');plot(towns,add=T,pch=3,col='red')

Create an empty raster grid with the properties of the desired raster map.

r <- raster(resolution=c(0.01,0.01),        # cell size: here 10m x 10m as units in km
            ext=extent(FGArea),             # extent/size of raster: set to fishing polygon area
            crs=crs(FGArea))                # coordinate system

Create a raster that transfers the values from the fishing area shapefile to the raster grid

TotYd <- rasterize(FGArea,                  # fishing area shapefile
                   r,                       # gridded raster
                   FGArea$Yield/1000,       # values to assign to the grid cells (Yield (kg) to metric tons)
                   fun=sum)                 # sum all values from overlapping polygon features
Creating maps

Example plot using the base plotting functions

plot(TotYd,col=matlab.like2(50),ext=extent(FGArea),axes=F)
scalebar(10, divs = 4,type = "bar", adj=c(0.7, -0.5))

Example plot using the tmap** package (Figure 2 in main paper)**
mypalette <- matlab.like2(100)              # colour scheme for map

png(file="Figure1.png", width = 8, height = 6, units = 'in',res = 600) # create hi-res png file

# Creates map object by progressively adding layers
  # Fishing ranges polygon
tm_shape(TotYd, unit.size = 1)+              
  tm_raster(palette = mypalette,title="Fishing intensity /n(mt/km/yr)",
            style='cont',auto.palette.mapping=F) +
  # Add land outline
  tm_shape(land)+
  tm_polygons(border.col='black',col='white')+
  # Add no fishing areas
  tm_shape(Nofishing)+
  tm_polygons('Id',labels = 'No fishing areas',
              palette='grey',border.col='black',title='')+
  tm_legend(position=c('right','middle'),title.size = 1,scale = 1.2)+
  # Add communities (with customized label locations)
  tm_shape(towns)+
  tm_dots(size=0.09)+
  tm_text('Name', size=1,ymod=0.3,
          xmod=c(3.6,5.25,2.75,4.5,-1,3.25,3.75),just='right')+
  # Add legend, compass, and scale bar
  tm_legend(position=c('right','top'))+
  tm_compass(position=c(0.9,0.9))+
  tm_scale_bar(position=c(0,0),width = 0.3)

 dev.off()     # end graphic file creation
Figure 2. Variability in annual commercial fishing intensity (mt km-2yr-1) for reef-associated species on the nearshore shallow shelf in Barbados.

Figure 2. Variability in annual commercial fishing intensity (mt km-2yr-1) for reef-associated species on the nearshore shallow shelf in Barbados.

Example plot using the RasterVis package (Figure 3 in main paper)

Create a raster brick to combine yield estimates for each season

# Create gridded rasters using seasonal yield data
NPlgcSeasonYd <- rasterize(FGArea,r,FGArea$NPlgcSeasonYd/1000,fun=sum)  
PlgcSeasonYd <- rasterize(FGArea, r,FGArea$PlgcSeasonYd/1000,fun=sum)
# Set zero values as blank (optional)
NPlgcSeasonYd[(NPlgcSeasonYd)==0] <- NA                                 
PlgcSeasonYd[(PlgcSeasonYd)==0] <- NA
# Create raster brick from the 2 rasters
Seasonbrick <- brick(NPlgcSeasonYd,PlgcSeasonYd,CRS('+proj=longlat +ellps=WGS84'))

Plot season graph using RasterVis.

#(*Optional*:  Define colour legend (warning: operates independently of map colours). Here the range of values and break points from the data are used)
minSeason <- min(minValue(Seasonbrick))
maxSeason <- max(maxValue(Seasonbrick))
miat <- seq(minSeason,maxSeason,length.out=100)
myTheme <- rasterTheme(region = matlab.like2(100),layout.heights=list(bottom.padding=4))
myColorkey <- list(at=miat, ## where the colors change
                   labels=list(at=seq(0,5,1),cex=1.3),space="bottom" ## where to print labels
)

# Create map object using the levelplot function
pSeason <- levelplot(Seasonbrick,                                         # rasterbrick object
                     par.settings=myTheme,                                # plot options (optional)
                     names.attr=c('Pelagic off-season','Pelagic season'), # panel names
                     par.strip.text = list(cex = 1.3),                    # size of text in panels
                     scales=list(draw=FALSE),                             # surpress marginal graphics
                     colorkey=myColorkey)                                 # settings for the colour legend

# add additional layers
pSeason + 
  spplot(land, lwd=2, col='black',fill='white') +
  spplot(Nofishing, lwd=1, col='black',fill='grey')


#----- Using Trellis for other labels (optional)
trellis.focus("legend", side="bottom", clip.off=F, highlight=FALSE)
grid.text(expression(paste("Fishing intensity (mt km"^"2",")")), 0.5, 0, hjust=0.5, vjust=1.5,rot = 0,gp=gpar(fontsize=14))
trellis.unfocus()

trellis.focus("panel", column=1, row=1,clip.off=F, highlight=T)
grid.text(expression('A'), 0.9, 0.95, hjust=0.5, vjust=1,rot = 0,gp=gpar(fontsize=14))
trellis.unfocus()

trellis.focus("panel", column=2, row=1, clip.off=F, highlight=T)
grid.text(expression('B'), 0.9, 0.95, hjust=0.5, vjust=1,rot = 0,gp=gpar(fontsize=14))
trellis.unfocus()
Figure 3. Spatial variation in reef fishing intensity (mt km-2) on Barbados’ nearshore shallow shelf during the (A) pelagic fishing off-season (Jun-Nov) and the (B) pelagic fishing season (Dec-May).

Figure 3. Spatial variation in reef fishing intensity (mt km-2) on Barbados’ nearshore shallow shelf during the (A) pelagic fishing ‘off-season’ (Jun-Nov) and the (B) pelagic fishing season (Dec-May).