In [1]:
# load dependencies
## geospatial data handling
library(rgdal)
library(raster)
library(sp)
library(mapdata)
library(maptools)
library(ncdf4)
## general data handling
library(XML)
library(dplyr)
library(tidyr)
library(reshape2)
library(downloader)
## plotting
library(directlabels)
library(rasterVis)
library(ggplot2)
## for display reasons
library(knitr)
library(kableExtra)
library(IRdisplay)
library(repr)
options(repr.plot.wdith=6, repr.plot.height=8)
Loading required package: sp

rgdal: version: 1.4-4, (SVN revision 833)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
 Path to GDAL shared files: /usr/share/gdal/2.2
 GDAL binary built with GEOS: TRUE 
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: (autodetected)
 Linking to sp version: 1.3-1 

Loading required package: maps

Checking rgeos availability: FALSE
 	Note: when rgeos is not available, polygon geometry 	computations in maptools depend on gpclib,
 	which has a restricted licence. It is disabled by default;
 	to enable gpclib, type gpclibPermit()


Attaching package: ‘dplyr’


The following objects are masked from ‘package:raster’:

    intersect, select, union


The following objects are masked from ‘package:stats’:

    filter, lag


The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union



Attaching package: ‘tidyr’


The following object is masked from ‘package:raster’:

    extract



Attaching package: ‘reshape2’


The following object is masked from ‘package:tidyr’:

    smiths


Loading required package: lattice

Loading required package: latticeExtra

Loading required package: RColorBrewer


Attaching package: ‘ggplot2’


The following object is masked from ‘package:latticeExtra’:

    layer



Attaching package: ‘kableExtra’


The following object is masked from ‘package:dplyr’:

    group_rows



During the second part of the tutorial, we focus on the "Estuaires picards et mer d'Opale" Marine Protected Area (MPA) in France. The shapefile for this area can be downloaded from the World Database on Protected Areas (WDPA).

In [2]:
# download and unzip the shapefile
con = "http://d1gam3xoknrgr2.cloudfront.net/current/WDPA_WDOECM_protected_area_555559631_shp.zip"
download(con,"shapefile.zip", mode = "wb")
unzip("shapefile.zip",exdir="MPA2")
file.remove("shapefile.zip")
unzip("MPA2/WDPA_WDOECM_protected_area_555559631_shp0.zip", exdir = "MPA2")
# read it into R
mpa = readOGR("MPA2","WDPA_WDOECM_protected_area_555559631_shp-polygons")
# get the spatial extent of the MPA
xmin <- extent(mpa)@xmin
ymin <- extent(mpa)@ymin
xmax <- extent(mpa)@xmax
ymax <- extent(mpa)@ymax
TRUE
OGR data source with driver: ESRI Shapefile 
Source: "/home/jovyan/R-tutorial/MPA2", layer: "WDPA_WDOECM_protected_area_555559631_shp-polygons"
with 1 features
It has 30 fields

We make some maps to show the location of the MPA.

In [3]:
# Localisation at global scale
map("worldHires",col="light grey",fill=T)
points(coordinates(mpa),cex=2,col="blue",pch="+")
title(paste("Region of interest ",
            "( W-Lon",round(xmin,2),
            "S-Lat",round(ymin,2),
            "E-Lon",round(xmax,2),
            "N-Lat",round(ymax,2),")"),cex=.5)

# Localisation at local scale
plot(mpa,xlim=c(xmin-1,xmax+1),ylim=c(ymin-1,ymax+1),axes=T,col="red")
map("worldHires",add=T,col="light grey",fill=T)
plot(mpa,add=T,col="blue")
title(paste("MPA",mpa$ORIG_NAME),cex=.5)

We again extract the bathymetry using the function defined in part 1 of this tutorial.

In [4]:
# Define a function to read in raster data from the EMODnet bathymetry WCS
getbathymetry<-function (name = "emodnet:mean", xmin = 15, xmax = 20.5, ymin = 30, ymax = 32.5){
  bbox <- paste(xmin, ymin, xmax, ymax, sep = ",")
                              
  con <- paste("https://ows.emodnet-bathymetry.eu/wcs?service=wcs&version=1.0.0&request=getcoverage&coverage=",
               name,"&crs=EPSG:4326&BBOX=", bbox,
               "&format=image/tiff&interpolation=nearest&resx=0.00208333&resy=0.00208333", sep = "")
  
  print(con)
  
  stop   
  nomfich <- paste(name, "img.tiff", sep = "_")
  nomfich <- tempfile(nomfich)
  download(con, nomfich, quiet = TRUE, mode = "wb")
  img <- raster(nomfich)
  img[img == 0] <- NA
  img[img < 0] <- 0
  names(img) <- paste(name)
  return(img)
}

# get the bathymetry data for the MPA
bathy_img <- getbathymetry(name = "emodnet:mean", xmin, xmax, ymin, ymax)
bathy <- as.data.frame(as(bathy_img, "SpatialPixelsDataFrame"))
[1] "https://ows.emodnet-bathymetry.eu/wcs?service=wcs&version=1.0.0&request=getcoverage&coverage=emodnet:mean&crs=EPSG:4326&BBOX=0.785219445000052,50.051728573,1.68869100100005,50.8176805940001&format=image/tiff&interpolation=nearest&resx=0.00208333&resy=0.00208333"

Working with data from EMODnet biology

Access polygon data with WFS

Using the EMODnet biology WFS we can get the zooplankton abundances for different grid cells within the MPA (bbox). To calculate statistics of the abundance data over a particular period and a particular season, we can download it as a geoJSON file.

In [5]:
# set the request parameters
reg.name = 36317
spec.name = 'Acartia_spp'
trend_start = 2004
trend_end = 2013
season.code = 1
bbox <- paste(xmin, ymin, xmax, ymax, sep = ",")
type.name='Emodnetbio:OOPS_products_vliz'

# build the WFS URL
url<-"http://geo.vliz.be/geoserver/Emodnetbio/ows?service=WFS&version=1.0.0&request=GetFeature&outputFormat=json"
full_url<-paste(url,"&viewParams=;startYearCollection:",trend_start,";mrgid:",reg.name,";scientificName:",spec.name,";season:",season.code,"&typeName=",type.name,"&BBOX=",bbox,sep="")

# get the data as a geoJSON file
library("rjson")
json_data <- fromJSON(file=full_url)
# to display the file
# str(json_data)
# to save
# geojson_write(json_data, file = 'MyJson.geojson') #bogue ?

test<-json_data[['features']]
mean=0 
l<- length(test)
for (i in 1:l){
tests<-test[[l]]
mean<-mean+tests$properties$measurementValue
}
mean_value_season1<-mean/l
mean_value_season1
0.450216

We can also visualise the extent of the zooplankton abundance grid cell by downloading the data as a zipped shapefile and plotting it on a map.

In [6]:
# build the WFS URL
url<-"http://geo.vliz.be/geoserver/Emodnetbio/ows?service=WFS&version=1.0.0&request=GetFeature&outputFormat=shape-zip"
full_url<-paste(url,"&viewParams=;startYearCollection:",trend_start,";mrgid:",reg.name,";scientificName:",spec.name,";season:",season.code,"&typeName=",type.name,"&BBOX=",bbox,sep="")

# get the data as a zipped shapefile
download(full_url, dest="shapefile.zip", mode="wb") 
unzip("shapefile.zip", overwrite = TRUE)
OOPS_products <- readOGR(dsn = ".", layer = "OOPS_products_vlizPolygon")

# plot on a map
map<-ggplot()+
  theme_bw()+
  theme(panel.grid.minor.y= element_blank(), panel.grid.minor.x = element_blank())+
  geom_raster(data=bathy,aes(x=x,y=y,fill=emodnet.mean),alpha=.75)+
  scale_fill_gradient(low = "white", high = "darkblue",name="Depth (m)")+
  geom_polygon(data=OOPS_products,aes(x=long,y=lat,group=group, fill="OOPS_products"),colour="green",fill="blue",alpha=.1)+
  coord_quickmap(xlim=range(xmin,xmax),ylim=range(ymin,ymax))+
ggtitle("OOPS_products polygons")+xlab("Longitude")+ylab("Latitude")
plot(map)

# export to a shapefile which can then be visualised in e.g. QGIS
writeOGR(OOPS_products,"OOPS_products","OOPS_products_bbox",driver="ESRI Shapefile",overwrite_layer = TRUE)
OGR data source with driver: ESRI Shapefile 
Source: "/home/jovyan/R-tutorial", layer: "OOPS_products_vlizPolygon"
with 62 features
It has 10 fields
Regions defined for each Polygons

Alternatively, we can also download the zooplankton abundance grid cell as geoJSON.

In [7]:
# function to read GeoJSON
getgeojson<-function(name="Emodnetbio:OOPS_products_vliz",xmin=-1,xmax=1,ymin=49,ymax=50){
  reg.name = 36317 #label="Greater North Sea"
  trend_start = 1958 #startDate
  trend_end = 2013 #endDate
  spec.name = 'Acartia_spp'
  trend_start = 2004
  trend_end = 2013
  season.code = 1
  bbox<-paste(xmin,ymin,xmax,ymax,sep=",")
  
  url<-"http://geo.vliz.be/geoserver/Emodnetbio/ows?service=WFS&version=1.0.0&request=GetFeature&outputFormat=json"
  con<-paste(url,"&viewParams=;startYearCollection:",trend_start,";mrgid:",reg.name,";scientificName:",spec.name,";season:",season.code,"&typeName=",name,"&BBOX=",bbox,sep="")
  print(con)
  temp <- tempfile()
  download(con,temp)
  ogrInfo(dsn=temp)
  layer<-readOGR(dsn=temp)
  return(layer)
}
OOPS_products<-getgeojson(name="Emodnetbio:OOPS_products_vliz",xmin,xmax,ymin,ymax)

map <- ggplot()+
  theme_bw()+
  theme(panel.grid.minor.y= element_blank(), panel.grid.minor.x = element_blank()) +
  geom_raster(data=bathy,aes(x=x,y=y,fill=emodnet.mean),alpha=.75) +
  scale_fill_gradient(low = "white", high = "darkblue",name="Depth (m)") +
  geom_polygon(data=OOPS_products,aes(x=long,y=lat,group=group,fill="OOPS_products"),colour="green",fill="blue",alpha=.1) +
  coord_quickmap(xlim=range(xmin,xmax),ylim=range(ymin,ymax)) +
  ggtitle("OOPS_products polygons")+xlab("Longitude")+ylab("Latitude")

plot(map)
[1] "http://geo.vliz.be/geoserver/Emodnetbio/ows?service=WFS&version=1.0.0&request=GetFeature&outputFormat=json&viewParams=;startYearCollection:2004;mrgid:36317;scientificName:Acartia_spp;season:1&typeName=Emodnetbio:OOPS_products_vliz&BBOX=0.785219445000052,50.051728573,1.68869100100005,50.8176805940001"
OGR data source with driver: GeoJSON 
Source: "/tmp/Rtmp4Mike4/file15821153d44", layer: "file15821153d44"
with 62 features
It has 12 fields, of which 1 list fields
Regions defined for each Polygons

Access time series data with WFS

Using WFS, we can also access the zooplankton abundance time series as a csv file.

In [8]:
# determine parameters
reg.name = 36317 #label="Greater North Sea"
spec.name = 'Acartia_spp' #species
trend_start = 1958 #startDate
trend_end = 2013 #endDate

# construct WFS url
mainurl<-"http://geo.vliz.be/geoserver/Emodnetbio/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=Emodnetbio:OOPS_summaries&outputFormat=csv&SQL_FILTER=startYearCollection-endYearCollection=0"
fullurl<-paste(mainurl,"ANDscientificName='",spec.name,"'ANDmrgid=",reg.name,sep="")

# get the data
datans <- read.csv(file= fullurl, header=TRUE)

We now analyse the time series seasonality and trends using the stl() loess seasonal decomposition function.

In [9]:
levels(datans$scientificName)[levels(datans$scientificName)=="RatioLS copepods"] <- "Large/Small copepod ratio"
seldata<-datans[with(datans, order(datans$startYearCollection,datans$season)),]
# convert to a time series object
tscf<-ts(data=seldata$avg_ab, start=trend_start, end=trend_end, frequency=4)
# decompose the time series
fitstl<-stl(tscf,s.window="periodic",t.window=15)
plot(fitstl)