We will use the R packages listed below. They can be downloaded using the install.packages() command. In case you are running R on Windows, you may need to download Rtools to build the packages. We also install the esri2sf package from github to scrape the ArcGIS REST API.

# install.packages(c("icesSAG","sf","ggplot2","plotly","devtools"))
library(devtools)
install_github("yonghah/esri2sf")
# load dependencies
library(icesSAG)
library(sf)
library(esri2sf)
library(ggplot2)
library(plotly)

ICES offers web services to download tabular data as map services to access geospatial data.

1 Access ICES web services

In this excercise we will use the icesSAG package to access the web services of the ICES Stock Assessment Graphs database. In the script we get the summary data for all sandeel stocks published in 2018.

summary_data <- getSAG(stock="Sandeel", year=2018)
head(summary_data)
##   Year recruitment high_recruitment low_recruitment low_SSB SSB high_SSB
## 1 1970          NA               NA              NA      NA  NA       NA
## 2 1971          NA               NA              NA      NA  NA       NA
## 3 1972          NA               NA              NA      NA  NA       NA
## 4 1973          NA               NA              NA      NA  NA       NA
## 5 1974          NA               NA              NA      NA  NA       NA
## 6 1975          NA               NA              NA      NA  NA       NA
##   catches landings discards low_F  F high_F StockPublishNote Purpose Fage
## 1      NA        0       NA    NA NA     NA  Stock published  Advice <NA>
## 2      NA        0       NA    NA NA     NA  Stock published  Advice <NA>
## 3      NA        0       NA    NA NA     NA  Stock published  Advice <NA>
## 4      NA        0       NA    NA NA     NA  Stock published  Advice <NA>
## 5      NA        0       NA    NA NA     NA  Stock published  Advice <NA>
## 6      NA        0       NA    NA NA     NA  Stock published  Advice <NA>
##   fishstock recruitment_age AssessmentYear  units stockSizeDescription
## 1 san.27.6a              NA           2018 tonnes                 <NA>
## 2 san.27.6a              NA           2018 tonnes                 <NA>
## 3 san.27.6a              NA           2018 tonnes                 <NA>
## 4 san.27.6a              NA           2018 tonnes                 <NA>
## 5 san.27.6a              NA           2018 tonnes                 <NA>
## 6 san.27.6a              NA           2018 tonnes                 <NA>
##   stockSizeUnits fishingPressureDescription fishingPressureUnits
## 1           <NA>                       <NA>                 <NA>
## 2           <NA>                       <NA>                 <NA>
## 3           <NA>                       <NA>                 <NA>
## 4           <NA>                       <NA>                 <NA>
## 5           <NA>                       <NA>                 <NA>
## 6           <NA>                       <NA>                 <NA>

We now plot the recruitment evolution of the different sandeel stocks.

ggplot(summary_data[complete.cases(summary_data[c("Year", "recruitment")]),], 
       aes(x=Year, y=recruitment, group = fishstock, colour = fishstock)) +
    geom_line()

2 Access ICES map services

Along with the data services, ICES hosts services for spatial datasets. All of ICES reference layers such as the ICES Areas, ICES Statistical rectangles, ICES Ecoregions can be accessed through thier map server located at http://gis.ices.dk/gis/rest/services. This map server can be accessed either through the ArcGIS REST API or through OGC webservices.

In this tutorial we access vector data through the ArcGIS REST API using the esri2sf package which saves the data in an sf object. For the sake of this tutorial we will download the Datras Survey Areas available at http://gis.ices.dk/gis/rest/services/ICES_Datasets/DATRAS_GIS_full/MapServer/0, selecting only the AreaName and Description fields.

datras = esri2sf("http://gis.ices.dk/gis/rest/services/ICES_Datasets/DATRAS_GIS_full/MapServer/0",outFields = c("AreaName","Description"))
## [1] "Feature Layer"
## [1] "esriGeometryPolygon"
head(datras)
## Simple feature collection with 6 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -12.35178 ymin: 48.5 xmax: 10.30384 ymax: 59
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##    AreaName          Description                          geoms
## 1      VIIb  VIIb_Deep, 126-200m MULTIPOLYGON (((-10.78447 5...
## 2 BTS_total     BTS survey areas MULTIPOLYGON (((-4.186093 5...
## 3      VIIj VIIj_Slope, 201-600m MULTIPOLYGON (((-12.25435 5...
## 4      VIIb VIIb_Slope, 201-600m MULTIPOLYGON (((-12.35178 5...
## 5       VIa  VIa_Slope, 201-600m MULTIPOLYGON (((-9.134293 5...
## 6      VIIj VIIj_Medium, 81-125m MULTIPOLYGON (((-10.59832 5...

We now make an interactive map of the areas and description (note: on the website this is static)

# get european countries
datras_map = ggplot()+
            geom_sf(data = datras, aes(fill = AreaName, name=Description), size = 0.1)+
            borders(fill = "black", size = 0.1)+
            coord_sf(xlim=c(-20,30),ylim=c(35,65))+
            theme_bw()+
            theme(legend.position = "none",axis.title = element_blank())
ggplotly(datras_map,)