In this script we aim to explore different available datasets and tools
for producing a EU marine areas maps. It does not intend to be
exhaustive, but rather to make a sample script that I often (re)use and
might be useful to others. From a short quick map to put in data to
higher resolution maps, that may require downloading the data from the
internet.
From these bases, we can than add our data or costumize the maps
accordingly our taste and desires.
This first chunk provides data to use in most chunks and other
scripts. Part of this was obtained with chatGPT aid.
# Establish our three scales of geographic limits (bbox)
eu.lim = data.frame(xlim=c(-25,35), ylim=c(30,70))
eu.bbox <- st_bbox(c(xmin=eu.lim$xlim[1],
xmax=eu.lim$xlim[2],
ymin=eu.lim$ylim[1],
ymax=eu.lim$ylim[2]), crs=4326)
pt.lim = data.frame(ylim=c(36.6, 42), xlim=c(-10, -7.0))
pt.bbox <- st_bbox(c(xmin=pt.lim$xlim[1],
xmax=pt.lim$xlim[2],
ymin=pt.lim$ylim[1],
ymax=pt.lim$ylim[2]))
loc.lim = data.frame(ylim=c(40.6, 40.9), xlim=c(-8.8, -8.4))
loc.bbox <- st_bbox(c(
xmin=loc.lim$xlim[1],
xmax=loc.lim$xlim[2],
ymin=loc.lim$ylim[1],
ymax=loc.lim$ylim[2]))
# Filter for European EEZs (example using EU ISO codes)
eu.iso3 <- c("ALB", "BEL", "BGR", "HRV", "CYP", "CZE", "DNK", "EST", "FIN", "FRA","DEU", "GRC", "HUN", "IRL", "ITA", "LVA", "LTU", "LUX", "MLT", "NLD","POL", "PRT", "ROU", "SVK", "SVN", "ESP", "SWE", "GBR", "MNE", "NOR", "ISL")
# there is a gap between Lithuania and Poland - it belongs to Russia I think OMG
eu.eez <- c("Albania", "Belgium", "Bosnia and Herzegovina", "Bulgaria", "Croatia",
"Cyprus", "Denmark", "Estonia", "Finland", "France", "Germany", "Greece",
"Iceland", "Ireland", "Italy", "Latvia", "Lithuania", "Malta", "Monaco",
"Montenegro", "Netherlands", "Norway", "Poland", "Portugal", "Romania",
"Slovenia", "Spain", "Sweden","Ukraine", "United Kingdom","Russia")
# European Community
ec.eez <- c("Belgium","Bulgaria", "Croatia","Cyprus", "Denmark", "Estonia", "Finland", "France", "Germany", "Greece","Ireland", "Italy", "Latvia", "Lithuania", "Malta", "Monaco", "Netherlands", "Poland", "Portugal", "Romania","Slovenia", "Spain", "Sweden")
# EUn.1# EU country borders
eu.nocoast <- c("Austria", "Belarus", "Czech Republic", "Hungary", "Kosovo","Liechtenstein", "Luxembourg", "North Macedonia", "Moldova","Serbia", "Slovakia", "Switzerland", "Vatican", "Jersey","Guernsey")
# Define vector of EC countries (ISO or names)
ec.countries <- c("Belgium", "France", "Germany", "Portugal", "Spain", "Italy", "Ireland", "Denmark", "Netherlands", "Sweden", "Finland","Estonia", "Latvia", "Lithuania", "Poland", "Croatia","Slovenia", "Greece", "Malta", "Cyprus", "Bulgaria","Romania", "Czech Republic", "Slovakia", "Austria", "Hungary", "Luxembourg")
# select ICES coastal areas
ices.coastal.areas <- sort(unique(c(
"27.1", "27.2.a.2", "27.3.a.20", "27.3.a.21", "27.3.a.22", "27.3.b",
"27.3.b.23", "27.3.c.22", "27.3.d.24", "27.3.d.25", "27.3.d.26",
"27.3.d.27", "27.3.d.28.1", "27.3.d.28.2", "27.3.d.29", "27.3.d.30",
"27.3.d.32", "27.4.a", "27.4.b", "27.4.c", "27.4.d", "27.5.a",
"27.5.b", "27.6.a", "27.6.b", "27.7.a", "27.7.b", "27.7.c", "27.7.d",
"27.7.e", "27.7.f", "27.7.g", "27.7.h", "27.7.j.2", "27.8.a",
"27.8.b", "27.8.c", "27.9.a"
)))
# "27.2.a.1",
# Exclude: 27.9.b.2
# Define the list of GSA areas that have EU borders (based on knowledge of regions)
gsa.eu <- c("GSA 1", "GSA 5", "GSA 6", "GSA 7", "GSA 8", "GSA 9","GSA 10", "GSA 11.1", "GSA 11.2", "GSA 15", "GSA 16", "GSA 17","GSA 18", "GSA 19", "GSA 20", "GSA 22", "GSA 23", "GSA 25", "GSA 28", "GSA 29")
# Exclude: "GSA 4", "GSA 2","GSA 12", "GSA 14", "GSA 21.1","GSA 24", "GSA 26", "GSA 3", "GSA 13", "GSA 21.2", "GSA 21.3""GSA 27", "GSA 30"
Just a general map. Note this as some issues latter on - I need to
get a better resolution for this.
require(rnaturalearth)
# extract data
eu.map <- ne_countries(scale = 'large',
type = 'map_units', returnclass = 'sf') %>%
filter(
continent == 'Europe' | admin == "Turkey" | admin == "Cyprus"
#!admin %in% landlocked_countries
) %>%
dplyr::select(admin, adm0_a3, iso_a3, geounit, continent, region_un, region_wb, subregion, pop_est, abbrev, sovereignt, type) %>% st_transform(crs=4326) %>%
rename(iso3 = adm0_a3)
# #eu.map
# eu.map %>%
# st_crop(eu.bbox) %>%
# mapview(zcol="iso3")
# Static map
eu.map %>%
ggplot() +
geom_sf(aes(fill=iso3), col="grey", na.rm =TRUE)+
geom_sf_text(aes(label = iso3), size=2) +
coord_sf(xlim=c(-24,35), ylim=c(35,70))+
#scale_fill_brewer(palette = "Set3")+
#scale_fill_distiller(palette = "Set3", direction = 1, na.value = "transparent")+
ggtitle("Iso3") +
xlab("Longitude") + ylab("Latitude") +
theme(legend.position = "none")
Emodnet data downloaded from link;
“The present layer compiles the marine regions and subregions listed
in Article 4 of the Marine Strategy Framework Directive (MSFD), together
with other surrounding seas of Europe. The MSFD marine regions and
subregions map was developed to support DG Environment and EU Member
States in their implementation of the MSFD. It represents the current
state of understanding of the marine regions and subregions and is
subject to amendment in light of any new information which may be
produced”
File:Emodnet_HA_OtherManagementAreas_MSFDReportingUnits_20190422/Emodnet_HA_OtherManagementAreas_MSFDReportingUnits_20190422.shp\
# marine boundaries from emodnet
msfd.regions = st_read(paste0(MapsDir, "Emodnet_HA_OtherManagementAreas_MSFDReportingUnits_20190422/Emodnet_HA_OtherManagementAreas_MSFDReportingUnits_20190422.shp"), quiet = TRUE)
msfd.regions <- msfd.regions %>% st_transform(crs=4326) %>% select(ID, NAME)
# eu.regions <- eu.regions %>% st_make_valid()
# names(msfd.regions)
# msfd.regions %>% head()
#eu.regions$NAME
# msfd.regions %>%
# st_crop(xmin = -25, xmax = 35, ymin = 32, ymax = 70) %>%
# mapview(zcol="NAME", legend=FALSE)
# Static map
msfd.regions %>%
ggplot() +
geom_sf(aes(fill=NAME), col="grey", na.rm =TRUE) +
geom_sf_text(aes(label = NAME), size=2) +
coord_sf(xlim=c(-24,35), ylim=c(35,70))+
xlab("Longitude") + ylab("Latitude") +
theme(legend.position = "none")
Marineregions.org data downloaded from link;
This is a global map with 572 features; I have cutted to include only
our interest areas (~Europe).
This file would be: Maritime Boundaries Geodatabase: Maritime
Boundaries and Exclusive Economic Zones (200NM), version 12. Note there
are two shapes inthere.
# # marineregions EZZ
dim(eez <- st_read(paste0(MapsDir,"/marineregions/World_EEZ_v12_20231025/eez_v12.shp"), quiet = TRUE))
## [1] 285 32
#eez_boundaries_v12.shp
# filter for EUR
dim(eez <- eez %>%
filter(SOVEREIGN1 %in% eu.eez) %>%
rename(iso3=ISO_TER1) %>%
select(TERRITORY1, SOVEREIGN1, iso3, GEONAME))
## [1] 94 5
# 94
# make valid
eez <- eez %>% st_make_valid()
# cut by interest area
eez <- eez %>%
st_transform(crs=4326)
#st_crop(eu.bbox)
# dynamic map
# eez %>% mapview(zcol="SOVEREIGN1")
# Static map
eez %>%
filter(iso3!= "RUS") %>%
ggplot() +
geom_sf(aes(fill=iso3), col="grey", na.rm =TRUE) +
geom_sf_text(aes(label = iso3), size=2) +
coord_sf(xlim=c(-24,35), ylim=c(35,70))+
xlab("Longitude") + ylab("Latitude") +
theme(legend.position = "none")
Marineregions.org data downloaded from link;
The first file would be the intercection between IHO and EEZ.
File:Intersect_EEZ_IHO_v5_20241010/Intersect_EEZ_IHO_v5_20241010.shp\
# # marineregions EZZ intercepts with iho
dim(eez.iho <- st_read(paste0(MapsDir,"/marineregions/Intersect_EEZ_IHO_v5_20241010/Intersect_EEZ_IHO_v5_20241010.shp"), quiet = TRUE))
## [1] 572 35
#572 35
# names(eez.iho)
# head(eez.iho)
# filter for EUROPE only and sleect relevant cols
dim(eez.iho <- eez.iho %>%
filter(SOVEREIGN1 %in% eu.eez) %>%
rename(iso3=ISO_TER1) %>%
select(TERRITORY1, SOVEREIGN1, iso3, MarRegion, IHO_Sea, EEZ))
## [1] 181 7
# 181
# make valid
eez.iho <- eez.iho %>% st_make_valid()
# cut by interest area
eez.iho <- eez.iho %>%
st_transform(crs=4326) %>%
st_crop(eu.bbox)
# interactive map
# eez.iho %>% mapview(zcol="SOVEREIGN1",legend=FALSE)
# example of a couple of countries
# eez.iho %>%
# #filter(SOVEREIGN1 %in% c("Malta","Italy")) %>%
# #filter(SOVEREIGN1 %in% c("France","Italy")) %>% filter(SOVEREIGN1 %in% c("Portugal")) %>%
# mapview()
# Static map
eez.iho %>%
ggplot() +
geom_sf(aes(fill=iso3), col="grey", na.rm =TRUE) +
geom_sf_text(aes(label = iso3), size=2) +
coord_sf(xlim=c(-24,35), ylim=c(35,70))+
xlab("Longitude") + ylab("Latitude") +
theme(legend.position = "none")
European Environmental Assessment areas were downloaded from EEA
webpage directly.
“The EEA marine assessment areas are harmonised units used for
statistics. The EEA established the assessment areas for marine
assessments to avoid jurisdictional issues on disputed maritime borders.
The assessment areas dataset allows for a consistent geospatial approach
for providing and updating spatial statistics of relevance for the
implementation of EU legislation and progress monitoring of EU
policies.”
File:
EEA_Marine_Assessment_Areas/EEA_Marine_Assessment_Areas_LAEA_20221003.shp
# EEA Assessment areas
dim(eea.areas <- st_read(paste0(MapsDir,"EEA_Marine_Assessment_Areas/EEA_Marine_Assessment_Areas_LAEA_20221003.shp"), quiet = TRUE))
## [1] 10 12
# note it has this third weird geometry which is actually empty- weird behaviour when using mac, not occur in windows: Remove the Z dimension (3rd coordinate)
# eea.areas <- st_zm(eea.areas)
# note also: Projected CRS: ETRS89-extended / LAEA Europe
# names(eea.areas)
# eea.areas %>%
# st_drop_geometry() %>%
# kable() %>%
# kable_classic()
# make valid
eea.areas <- eea.areas %>% st_make_valid()
# transform into common crs
eea.areas <- eea.areas %>%
st_transform(crs=4326)
# no need to crop, already EUR
# %>% st_crop(eu.bbox)
# interactive map
# eea.areas %>% mapview(zcol="subregionN",legend=FALSE)
# Static map
eea.areas %>%
ggplot() +
geom_sf(aes(fill=subregionN), col="grey", na.rm =TRUE) +
geom_sf_text(aes(label = subregionN), size=2) +
coord_sf(xlim=c(-24,35), ylim=c(35,70))+
xlab("Longitude") + ylab("Latitude") +
theme(legend.position = "none")
Ices sub-areas were kindly provided to me by Josefine (DTU-Aqua). But
we can also download ICES statistical areas from marineregions.org
and problably from ICES.
File:
ICES_Areas_20160601_cut_dense_3857/ICES_Areas_20160601_cut_dense_3857.shp
dim(ices.areas <- st_read(paste0(MapsDir,"/ICES_Areas_20160601_cut_dense_3857/ICES_Areas_20160601_cut_dense_3857.shp"), quiet = TRUE))
## [1] 66 12
# 66 12
head(ices.areas %>% st_drop_geometry())
## OBJECTID_1 OBJECTID Major_FA SubArea Division SubDivisio Unit Area_Full
## 1 1 1 27 3 d 27 <NA> 27.3.d.27
## 2 2 2 27 3 d 25 <NA> 27.3.d.25
## 3 3 3 27 3 d 24 <NA> 27.3.d.24
## 4 4 4 27 3 b 23 <NA> 27.3.b.23
## 5 5 5 27 3 c 22 <NA> 27.3.c.22
## 6 6 6 27 5 a 2 <NA> 27.5.a.2
## Area_27 Area_km2 area_fdi
## 1 3.d.27 30398.142 27.3.d.27
## 2 3.d.25 43707.615 27.3.d.25
## 3 3.d.24 24167.589 27.3.d.24
## 4 3.b.23 2280.994 27.3.b.23
## 5 3.c.22 17862.100 27.3.c.22
## 6 5.a.2 362099.240 27.5.a.2
# transform crs
ices.areas <- ices.areas %>% st_transform(crs=4326)
# make it valid: takes time.
ices.areas <- ices.areas %>% st_make_valid()
## plot it
ices.areas %>%
filter(Area_Full %in% ices.coastal.areas) %>%
mapview(zcol="Area_Full")