EMODnetWFS Case Study: Accessing and mapping EMODnet data
Source:vignettes/articles/emodnetwfs.Rmd
emodnetwfs.RmdIntroduction
The package was designed to make EMODnet vector data layers easily
accessible in R. The package allows users to query information on and
download data from all available EMODnet
Web Feature Service (WFS) endpoints directly into their R working
environment. Data are managed as sf objects which
are currently the state-of-the-art in handling of vector spatial data in
R. The package also allows user to specify the coordinate reference
system of imported data.
Data Product
Installation
You can install the development version of EMODnetWFS
from GitHub with:
remotes::install_github("EMODnet/EMODnetWFS")Explore the EMODnet WFS services with R
For this tutorial we will make use of the sf,
dplyr and mapview packages. The simple
features sf package is a well known standard for dealing
with geospatial vector data. The package dplyr is a strong
library for data manipulation. This package also loads
magrittr’s pipe operator %>% (you could
also use the base
pipe), which allows to write pipelines in R. To visualize
geometries, mapview will create quick interactive maps.
Run this line to install these packages:
install.packages(c("sf", "dplyr", "mapview"))With the EMODnetWFS package, we can explore and combine
the data served by the EMODnet lots through OGC Web Feature
Services or WFS.
Imagine we are interested in seabed substrates. The first step is to
choose what EMODnet lot can provide with these data. For that, we can
check the services available with the emodnet_wfs()
function.
library(EMODnetWFS)
library(mapview)
library(dplyr)
library(sf)
emodnet_wfs()
#> service_name
#> 1 bathymetry
#> 2 biology
#> 3 biology_occurrence_data
#> 4 chemistry_cdi_data_discovery_and_access_service
#> 5 chemistry_cdi_distribution_observations_per_category_and_region
#> 6 chemistry_contaminants
#> 7 chemistry_marine_litter
#> 8 geology_coastal_behavior
#> 9 geology_events_and_probabilities
#> 10 geology_marine_minerals
#> 11 geology_sea_floor_bedrock
#> 12 geology_seabed_substrate_maps
#> 13 geology_submerged_landscapes
#> 14 human_activities
#> 15 physics
#> 16 seabed_habitats_general_datasets_and_products
#> 17 seabed_habitats_individual_habitat_map_and_model_datasets
#> service_url
#> 1 https://ows.emodnet-bathymetry.eu/wfs
#> 2 https://geo.vliz.be/geoserver/Emodnetbio/wfs
#> 3 https://geo.vliz.be/geoserver/Dataportal/wfs
#> 4 https://geo-service.maris.nl/emodnet_chemistry/wfs
#> 5 https://geo-service.maris.nl/emodnet_chemistry_p36/wfs
#> 6 https://nodc.ogs.trieste.it/geoserver/Contaminants/wfs
#> 7 https://www.ifremer.fr/services/wfs/emodnet_chemistry2
#> 8 https://drive.emodnet-geology.eu/geoserver/tno/wfs
#> 9 https://drive.emodnet-geology.eu/geoserver/ispra/wfs
#> 10 https://drive.emodnet-geology.eu/geoserver/gsi/wfs
#> 11 https://drive.emodnet-geology.eu/geoserver/bgr/wfs
#> 12 https://drive.emodnet-geology.eu/geoserver/gtk/wfs
#> 13 https://drive.emodnet-geology.eu/geoserver/bgs/wfs
#> 14 https://ows.emodnet-humanactivities.eu/wfs
#> 15 https://prod-geoserver.emodnet-physics.eu/geoserver/ows
#> 16 https://ows.emodnet-seabedhabitats.eu/geoserver/emodnet_open/wfs
#> 17 https://ows.emodnet-seabedhabitats.eu/geoserver/emodnet_open_maplibrary/wfsThe column service_name shows services available, while
service_url has the corresponding base url to perform a WFS
request. The Seabed portal should have the data we are looking for. A
WFS client can be created by passing the corresponding
service_name to the function
emodnet_init_wfs_client(). The layers available to this WFS
client are consulted with emodnet_get_wfs_info().
seabed_wfs_client <- emodnet_init_wfs_client(service = "seabed_habitats_general_datasets_and_products")
#> Loading ISO 19139 XML schemas...
#> Loading ISO 19115 codelists...
#> ✔ WFS client created successfully
#> ℹ Service: "https://ows.emodnet-seabedhabitats.eu/geoserver/emodnet_open/wfs"
#> ℹ Version: "2.0.0"
emodnet_get_wfs_info(wfs = seabed_wfs_client)
#> # A tibble: 36 × 9
#> # Rowwise:
#> data_source service_name service_url layer_name title abstract class format
#> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
#> 1 emodnet_wfs seabed_habita… https://ow… art17_hab… 2013… "Gridde… WFSF… sf
#> 2 emodnet_wfs seabed_habita… https://ow… art17_hab… 2013… "Gridde… WFSF… sf
#> 3 emodnet_wfs seabed_habita… https://ow… art17_hab… 2013… "Gridde… WFSF… sf
#> 4 emodnet_wfs seabed_habita… https://ow… art17_hab… 2013… "Gridde… WFSF… sf
#> 5 emodnet_wfs seabed_habita… https://ow… art17_hab… 2013… "Gridde… WFSF… sf
#> 6 emodnet_wfs seabed_habita… https://ow… art17_hab… 2013… "Gridde… WFSF… sf
#> 7 emodnet_wfs seabed_habita… https://ow… art17_hab… 2013… "Gridde… WFSF… sf
#> 8 emodnet_wfs seabed_habita… https://ow… art17_hab… 2013… "Gridde… WFSF… sf
#> 9 emodnet_wfs seabed_habita… https://ow… art17_hab… 2018… "Gridde… WFSF… sf
#> 10 emodnet_wfs seabed_habita… https://ow… biogenic_… Biog… "This l… WFSF… sf
#> # ℹ 26 more rows
#> # ℹ 1 more variable: layer_namespace <chr>Each layer is explained in the abstract column. We can
see several layers with the information provided by the EU member states
for the Habitats
Directive 92/43/EEC reporting. We will select the layers about
coastal lagoons, mudflats and sandbanks with their respective
layer_name.
habitats_directive_layer_names <- c("art17_hab_1110", "art17_hab_1140", "art17_hab_1150")
emodnet_get_layer_info(
wfs = seabed_wfs_client,
layers = habitats_directive_layer_names
)
#> # A tibble: 3 × 9
#> # Rowwise:
#> data_source service_name service_url layer_name title abstract class format
#> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
#> 1 emodnet_wfs https://ows.em… seabed_hab… art17_hab… 2013… "Gridde… WFSF… sf
#> 2 emodnet_wfs https://ows.em… seabed_hab… art17_hab… 2013… "Gridde… WFSF… sf
#> 3 emodnet_wfs https://ows.em… seabed_hab… art17_hab… 2013… "Gridde… WFSF… sf
#> # ℹ 1 more variable: layer_namespace <chr>We are now ready to read the layers into R with
emodnet_get_layers(). EMODnetWFS reads the geometries as
simple features (See sf package) transformed to 4326 by default. Specifying another map
projection is possible by passing a EPGS code or projection string with
emodnet_get_layers(crs = "your projection"). The argument
reduce_layers = TRUE stack all the layers in one single
tibble. Default is FALSE and returns a list of sf objects, one per
layer.
habitats_directive_layers <- emodnet_get_layers(
wfs = seabed_wfs_client,
layers = habitats_directive_layer_names,
reduce_layers = TRUE
)
class(habitats_directive_layers)
#> [1] "sf" "data.frame"
glimpse(habitats_directive_layers)
#> Rows: 221
#> Columns: 9
#> $ gml_id <chr> "art17_hab_1110.13", "art17_hab_1110.22", "art17_h…
#> $ habitat_code <chr> "1110", "1110", "1110", "1110", "1110", "1110", "1…
#> $ ms <chr> "DK", "ES", "ES", "PT", "PT", "PL", "DK", "FR", "U…
#> $ region <chr> "ATL", "MAC", "MMAC", "MMAC", "MATL", "MBAL", "MBA…
#> $ cs_ms <chr> "U2+", "U1+", "U1+", "XX", "U1-", "U1-", "U1-", "U…
#> $ country_code <chr> "Denmark", "Spain", "Spain", "Portugal", "Portugal…
#> $ habitat_code_uri <chr> "http://dd.eionet.europa.eu/vocabulary/art17_2018/…
#> $ habitat_description <chr> "Sandbanks which are slightly covered by sea water…
#> $ geom <MULTISURFACE [m]> MULTISURFACE (POLYGON ((420..., MULTI…Run the following code to have a quick look at the layers geometries
# Transform to Polygon geometry type from Multisurface
if (unique(st_geometry_type(habitats_directive_layers)) == "MULTISURFACE") {
habitats_directive_layers <- habitats_directive_layers %>%
st_cast(to = "GEOMETRYCOLLECTION") %>%
st_collection_extract(type = "POLYGON")
}
# Visualize
map <- mapview(habitats_directive_layers, zcol = "habitat_description", burst = TRUE)
mapFurthermore, we can get data from other EMODnet lots and combine them. The Human Activities portal provides the maritime boundaries of the European Union state members. This time we will not initiate a WFS client, but we will use the service name. The WFS client will be generated on the fly.
Same as before, we have a look at the layers available first.
emodnet_get_wfs_info(service = "human_activities")
#> ✔ WFS client created successfully
#> ℹ Service: "https://ows.emodnet-humanactivities.eu/wfs"
#> ℹ Version: "2.0.0"
#> # A tibble: 99 × 9
#> # Rowwise:
#> data_source service_name service_url layer_name title abstract class format
#> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
#> 1 emodnet_wfs human_activit… https://ow… activelic… Acti… The dat… WFSF… sf
#> 2 emodnet_wfs human_activit… https://ow… advisoryc… Advi… This sh… WFSF… sf
#> 3 emodnet_wfs human_activit… https://ow… aquacultu… Advi… This sh… WFSF… sf
#> 4 emodnet_wfs human_activit… https://ow… baltic Advi… This sh… WFSF… sf
#> 5 emodnet_wfs human_activit… https://ow… blacksea Advi… This sh… WFSF… sf
#> 6 emodnet_wfs human_activit… https://ow… longdista… Advi… This sh… WFSF… sf
#> 7 emodnet_wfs human_activit… https://ow… market Advi… This sh… WFSF… sf
#> 8 emodnet_wfs human_activit… https://ow… mediterra… Advi… This sh… WFSF… sf
#> 9 emodnet_wfs human_activit… https://ow… northsea Advi… This sh… WFSF… sf
#> 10 emodnet_wfs human_activit… https://ow… northwest… Advi… This sh… WFSF… sf
#> # ℹ 89 more rows
#> # ℹ 1 more variable: layer_namespace <chr>The layer_name for the maritime
boundaries seems to be maritimebnds. This dataset was
developed based on the official data
provided by the European Environmental Agency and the Maritime Boundaries
Database compiled by MarineRegions.org (Flanders Marine Institute,
2019).
The sitename variable specifies the type of boundary
each feature represents. For illustration purposes, we will filter our
request to return only Territorial
Seas.
maritime_boundaries <- emodnet_get_layers(
service = "human_activities",
layers = "maritimebnds",
reduce_layers = TRUE,
cql_filter = "sitename='Territory sea (12 nm)'"
)
#> ✔ WFS client created successfully
#> ℹ Service: "https://ows.emodnet-humanactivities.eu/wfs"
#> ℹ Version: "2.0.0"
glimpse(maritime_boundaries)
#> Rows: 64
#> Columns: 13
#> $ gml_id <chr> "maritimebnds.54", "maritimebnds.55", "maritimebnds.56", "m…
#> $ objectid <dbl> 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 39,…
#> $ mblszotpid <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
#> $ localid <dbl> 49036, 49042, 49063, 49064, 49065, 49066, 49087, 49099, 491…
#> $ sitename <chr> "Territory sea (12 nm)", "Territory sea (12 nm)", "Territor…
#> $ legalfound <date> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
#> $ legalfou_1 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
#> $ country <chr> "United Kingdom", "France", "France", "France", "France", "…
#> $ nationalle <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
#> $ nutscode <chr> "UK", "FR", "FR", "FR", "FR", "FR", "ES", "CY", "UK", "UK",…
#> $ mblsds_mbl <chr> "In www.marineregions.org", "In www.marineregions.org", "In…
#> $ shape_leng <dbl> 141.545470, 70.017156, 5.116398, 1.949563, 2.072975, 1.8378…
#> $ the_geom <MULTICURVE [°]> MULTICURVE (LINESTRING (-13..., MULTICURVE (LINE…Add the maritime boundaries to the previous map with this line. To do
so, cats to MULTILINESTRING geometry type, transform to the
same crs as the habitats directive sf (only required
because of the next step), and crop using the bounding box of the
habitats directive sf.