Recreating the ICEP from Sarelli et al (2018) with CMEMS
Introduction
Overview
This notebook aims to reproduce the results in Sarelli et al (2018): creating an Index of Coastal Eutrophication (ICEP) using the Copernicus Marine Environment Monitoring Service (CMEMS). The authors developped an automated methodology to estimate ICEP from the environmental variables:
- nitrate concentrations
- phosphate concentrations
- silica concentrations
- chlorophyll concentrations
- the euphotic zone depth
Motivations
One of the United Nations (UN) Sustainable Development Goal (SDG) is focuses on the oceans: "SDG 14 - Conserve and sustainably use the oceans, seas and marine resources for sustainable development". The very first of its 10 targets is "Indicator 14.1.1: Index of coastal eutrophication and floating plastic debris density". The ICEP is however not precisely defined and Sarelli et al (2018), among others, have proposed methodologies to estimate ICEP from Earth Observation (EO) systems like Copernicus' Marine Service.
The Index for Coastal Eutrophication Potential (ICEP) is an indicator for the potential of riverine nutrient export to sustain new production of non-diatoms phytoplankton biomass, including Harmful Algal Blooms. It uses nitrogen (N), phosphorus (P) and silica (Si) values of riverine loadings.
Another reference is “3.10 Coastal Eutrophication Potential” in Hofste et al, 2019 (World Resources Institute).
Note that the Copernicus data is only available at a resolution of 1/12$^\circ$: it is too coarse to quantity most conservation efforts (unless they are at the oceanic mesoscale, which is rare). It is however a great approach to obtain a baseline ICEP for a region.
Workflow: 01 Retrieving the CMEMS data
First we retrieved the CMEMS data: this can be manually downloaded via CMEMS My Ocean Viewer using the “Atlantic-Iberian Biscay Irish- Ocean Biogeochemical Analysis and Forecast” product.
This is the product page: NON ASSIMILATIVE Hindcast - IBI_MULTIYEAR_BGC_005_003.
- [x] To do: retrieve the data using the Copernicus API through the Motu client.
Here, we use the motuclient to retrieve the data. NOTE: you need a registered username!
References
Sarelli, A., Sykas, D., Miltiadou, M., Bliziotis, D., Spastra, Y. and Ieronymaki, M., 2018, August. A novel automated methodology that estimates the United Nations (UN) Sustainable Development Goal (SDG) 14.1. 1.: index of coastal eutrophication using the Copernicus Marine Environment Monitoring Service (CMEMS). In Sixth International Conference on Remote Sensing and Geoinformation of the Environment (RSCy2018) (Vol. 10773, p. 1077302). International Society for Optics and Photonics.doi:10.1117/12.2326160
Hofste, R.W., Kuzma, S., Walker, S., Sutanudjaja, E.H., Bierkens, M.F., Kuijper, M.J., Sanchez, M.F., Van Beek, R., Wada, Y., Rodríguez, S.G. and Reig, P., 2019. Aqueduct 3.0: Updated decision-relevant global water risk indicators. World Resources Institute: Washington, DC, USA. Available online at: https://www.wri.org/publication/aqueduct-30.
import os
Source directory where CMEMS file(s) are saved
source_dir = "data_CMEMS"
if not os.path.exists(source_dir):
os.mkdir(source_dir)
print("The downloaded files will be saved in the following directory:")
print(source_dir)
# utilities packages
import datetime
from netCDF4 import Dataset
import pandas as pd
# plotting packages
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import cmocean
import seaborn as sns
import cartopy.crs as ccrs
import cartopy.feature as cfeature
ICEP calculations definitions
from eutrophication_utils import *
Year as variable to save results to their respective folders
year = 2017
Out directory to save results
out_dir = "Sarelli_et_al2018_Aug" + str(year) + "_example_out"
if not os.path.exists(out_dir):
os.mkdir(out_dir)
save_plot = False # saves plots to out_dir
show_all = True # show all the plots, at all steps
Workflow
01 Retrieving the CMEMS data
We use the motuclient:
--motu https://my.cmems-du.eu/motu-web/Motu \
--service-id IBI_MULTIYEAR_BGC_005_003-TDS \
--product-id cmems_mod_ibi_bgc_my_0.083deg-3D_P1D-m \This refers to the specific CMEMS product mentioned above. Note that for dates between 2020-05-02 12:00:00 and 2022-05-14 12:00:00 we would use IBI_ANALYSISFORECAST_BGC_005_004-TDS instead.
To illustrate how to recreate the results from Sarelli et al (2018), we first focus on the "Week 1" subpanel in Figure 6: the first week of August 2017. The authors did not provide exact spatio-temporal bounds, so I guessed the following:
--date-min '2017-08-03 00:00:00' \
--date-max '2017-08-10 23:59:59' \
--longitude-min -7 --longitude-max 3.5 --latitude-min 49.5 --latitude-max 55.0 \We also download surface data only:
--depth-min 0.5057600140571594 --depth-max 0.5057600140571594We download all the following variables at once:
--variable no3 \
--variable po4 \
--variable si \
--variable chl \
--variable zeu \Finally we save the downloaded file to the specified out-dir and out-name.
!python3 -m motuclient \
--motu https://my.cmems-du.eu/motu-web/Motu \
--service-id IBI_MULTIYEAR_BGC_005_003-TDS \
--product-id cmems_mod_ibi_bgc_my_0.083deg-3D_P1D-m \
--date-min '2017-08-03 00:00:00' \
--date-max '2017-08-10 23:59:59' \
--depth-min 0.5057600140571594 \
--depth-max 0.5057600140571594 \
--variable no3 \
--variable po4 \
--variable si \
--variable chl \
--variable zeu \
--out-dir 'data_CMEMS' \
--out-name 'chl_no3_po4_si_zeu_201708.nc' \
--longitude-min -7 --longitude-max 3.5 --latitude-min 49.5 --latitude-max 55.0 \
--user <YOUR_USERNAME> --pwd <YOUR_PASSWORD>
Import local data from source_dir.
fnames = os.listdir(source_dir)
We want to import all the variables from all the files automatically. We know that latitude, longitude, depth and time will be recurrent independent variables. We want to import the data as "the variables that are NOT these independent variables".
indep_var = ["depth", "latitude", "time", "longitude"]
We loop through the files, then loop through the netCDF4 variables, find the variable of interest and make a dictionary out of it. Note: in a netCDF4 file with Python, we can get all the variables names with fh.variables.keys(). Same for dimensions.
for fname in fnames:
# load the netCDF4 file
fh = Dataset(os.path.join(source_dir,fname), 'r')
# now we loop through individual variables
for variable in fh.variables.keys():
# if this is the variable (eg chlorophyll, nitrate)
if variable not in indep_var:
# transform the dict_keys intro attributes with a 'variable" dictionary
exec(variable+' = xr.open_dataset(os.path.join(source_dir,fname))')
print("The variable " + variable + " was imported.")
print("The resolution is: ")
print("* " + "{:.4f}".format((chl.latitude[1]-chl.latitude[0]).data) + " degrees of latitude;")
print("* " + "{:.4f}".format(((chl.latitude[1]-chl.latitude[0]).data)*111) + " km;")
print("* " + "{:.4f}".format(((chl.latitude[1]-chl.latitude[0]).data)*111*0.6213712) + " mi.")
# This is helpful in e.g. case po4.latitude.shape != chl.latitude.shape
# dims = fh.dimensions
# # now we get individual dimensions
# for dim in dims.keys():
# temp_dim = fh.dimensions[dim].name
# temp_num = fh.dimensions[dim].size
# exec('N'+temp_dim + '= ' + str(temp_num))
# print(dim)
02 Water Area Classification Workflow
"The first step of the methodology contains the normalization of all the products to the parameters used in the algorithm. Concentration of chlorophyll-a and nutrients (silica, phosphate and nitrates) are provided directly from the CMEMS product. For the calculation of Secchi Depth (SD) as a measure of water transparency, a formula that combines euphotic zone depth with the SD of a water region is used according to the equations"
For the euphotic zone depth ($Z_{eu}$) and the Secchi Depth ($Z_{SD}$) in meters as a measure of water transparency:
linear_fn = True
2.1 Get index for nitrate, phosphorus, silicate concentrations
The concentrations of nitrogen N ($C_N$), of silica Si ($C_{Si}$) and of phosphorus ($C_P$) are in mmol/m$^3$.
print(po4.po4.long_name + ": " + po4.po4.units)
print(no3.no3.long_name + ": " + no3.no3.units)
print(si.si.long_name + ": " + si.si.units)
We now verify that the coordinates are the same. We can inspect visually with plt.plot or just use boolean.
if (po4.longitude == no3.longitude).all() & (po4.longitude == si.longitude).all():
print("The longitude coordinates are compatible!")
else:
print("Issues with longitude compatibility. Check out https://docs.xarray.dev/en/stable/generated/xarray.combine_by_coords.html to auto-magically combine the datasets!")
if (po4.latitude == no3.latitude).all() & (po4.latitude == si.latitude).all():
print("The latitude coordinates are compatible!")
else:
print("Issues with latitude compatibility. Check out https://docs.xarray.dev/en/stable/generated/xarray.combine_by_coords.html to auto-magically combine the datasets!")
if (po4.depth == no3.depth).all() & (po4.depth == si.depth).all():
print("The depth coordinates are compatible!")
else:
print("Issues with depth coordinate compatibility. Check out https://docs.xarray.dev/en/stable/generated/xarray.combine_by_coords.html to auto-magically combine the datasets!")
ind_N = get_indN(no3.no3.isel(depth = 0, time=0),
si.si.isel(depth = 0, time=0),
po4.po4.isel(depth = 0, time=0))
2.2 Get index for chlorophyll
Mass concentration of chlorophyll ($chl$) in mg/m$^3$.
print(chl.chl.long_name + ": " + chl.chl.units)
ind_chl = get_ind_chl(chl.chl.isel(depth = 0, time=0))
2.3 Get index for Secchi Depth
First we get the Secchi Depth from the Euphotic Zone Depth.
Z_SD = get_Z_SD(zeu);
It is likely that this is the main source of discrepancy between our results and the plots from Sarelli et al (2018).
ind_Z_SD = get_ind_ZSD(Z_SD.zeu.isel(time=0))
03 Get eutrophication index
index_eutroph = get_index(ind_N, ind_chl, ind_Z_SD)
Save output to out_dir.
index_eutroph.to_netcdf(out_dir + "/" + "index_eutroph_" + str(year) + ".nc")
04 Plots
Create longitude, latitude, datetime arrays etc. for ease of plotting.
Lon, Lat = np.meshgrid(zeu.longitude, zeu.latitude)
ts = pd.to_datetime(str(zeu.time[0].time.data))
d = ts.strftime('%Y-%m-%d')
Plots of input variables
Discover the beauty of xarray.
if show_all:
po4.po4.isel(depth = 0).plot(cmap = cmocean.cm.amp, col="time", col_wrap=4);
if show_all:
no3.no3.isel(depth = 0).plot(cmap = cmocean.cm.matter, col="time", col_wrap=4, vmax = 50);
if show_all:
si.si.isel(depth = 0).plot(cmap = cmocean.cm.turbid, col="time", col_wrap=4);
if show_all:
chl.chl.isel(depth = 0).plot(cmap = cmocean.cm.algae, col="time", col_wrap=4);
if show_all:
zeu.zeu.plot(cmap = sns.cubehelix_palette(start=.5, rot=-.75, as_cmap=True), col="time", col_wrap=4);
if show_all:
Z_SD.zeu.plot(cmap = sns.cubehelix_palette(start=.5, rot=-.75, as_cmap=True), col="time", col_wrap=4);
Plots of the different indices
plt.pcolor(ind_N, cmap = sns.cubehelix_palette(as_cmap=True))
plt.gca().set_aspect('equal', adjustable='box')
plt.colorbar()
plt.tight_layout()
plt.title("2.1 Index for nutrients", fontsize = 20);
plt.pcolor(ind_chl, cmap = sns.cubehelix_palette(start=2, rot=0, dark=0, as_cmap=True))
plt.gca().set_aspect('equal', adjustable='box')
plt.colorbar()
plt.tight_layout()
plt.title("2.2 Index for chlorophyll", fontsize = 20);
plt.pcolor(ind_chl, cmap = sns.cubehelix_palette(start=.5, rot=-.5, as_cmap=True))
plt.gca().set_aspect('equal', adjustable='box')
plt.colorbar()
plt.tight_layout()
plt.title("2.3 Index for euphotic depth", fontsize = 20);
ICEP plots
pcm = plt.pcolor(index_eutroph, cmap = 'Spectral_r');
plt.colorbar(pcm);
plt.title("Out-of-the-box results for ICEP");
fig, axs = plt.subplots(2, 3, figsize=(20,10))
# Ind N
pcm = axs[0,0].pcolor(ind_N, cmap = sns.cubehelix_palette(as_cmap=True))
axs[0,0].set_title("Index N")
axs[0,0].set_aspect('equal', adjustable='box')
fig.colorbar(pcm, ax=axs[0,0], orientation="horizontal")
# Ind chlorophyll
pcm = axs[0,1].pcolor(ind_chl, cmap = sns.cubehelix_palette(start=.5, rot=-.75, as_cmap=True))
axs[0,1].set_title("Index Chlorophyll")
axs[0,1].set_aspect('equal', adjustable='box')
fig.colorbar(pcm, ax=axs[0,1], orientation="horizontal")
# Ind depth
pcm = axs[0,2].pcolor(ind_Z_SD, cmap = sns.cubehelix_palette(start=.5, rot=-.5, as_cmap=True))
axs[0,2].set_title("Index Secchi depth")
axs[0,2].set_aspect('equal', adjustable='box')
fig.colorbar(pcm, ax=axs[0,2], orientation="horizontal")
# make single eutrophication index plot
axs[1,0].remove()
axs[1,2].remove()
pcm = axs[1,1].pcolor(index_eutroph, cmap = 'Spectral_r')
axs[1,1].set_title("Eutrophication level")
axs[1,1].set_aspect('equal', adjustable='box')
fig.colorbar(pcm, ax=axs[1,1], orientation="horizontal")
plt.suptitle("Water quality indices and ICEP", fontsize = 20);
fig.tight_layout()
if save_plot:
plt.savefig(out_dir + "/" + "Water quality indices and ICEP")
# This plots the matrix with ALL indices, including ICEP==zero (e.g. offshore)
# following the colormap from Sarelli et al (2018).
# This does not highlight actionable data for eutrophication as ICEP is coastal
# and offshore values will tend to zero. It however helps check the subsequent
# calculations when we aim to create a time series.
if show_all:
plt.figure(figsize=(10,7))
# Creates the map
mymap = plt.axes(projection=ccrs.PlateCarree())
mymap.add_feature(cfeature.LAND, alpha=0.5)
mymap.add_feature(cfeature.OCEAN, alpha=0.5)
mymap.add_feature(cfeature.BORDERS, linestyle='-', alpha=0.5)
mymap.add_feature(cfeature.RIVERS)
mymap.xaxis.set_visible(True)
mymap.yaxis.set_visible(True)
# Plots the data onto map
plt.scatter(Lon[index_eutroph<=2], Lat[index_eutroph<=2],
marker = 's', s= 10,
transform=ccrs.PlateCarree(),
c = 'limegreen')
plt.scatter(Lon[index_eutroph==3], Lat[index_eutroph==3],
marker = 's', s= 10,
transform=ccrs.PlateCarree(),
c = 'yellow')
plt.scatter(Lon[(index_eutroph>3)&(index_eutroph<=5)], Lat[(index_eutroph>3)&(index_eutroph<=5)],
marker = 's', s= 10,
transform=ccrs.PlateCarree(),
c = 'orange')
plt.scatter(Lon[index_eutroph>5], Lat[index_eutroph>5],
marker = 's', s= 10,
transform=ccrs.PlateCarree(),
c = 'red')
# Plot labels
plt.ylabel("Latitude", fontsize=14)
plt.xlabel("Longitude", fontsize=14)
plt.legend()
# Another custom plot for fun.
plt.figure(figsize=(20,14))
# Creates the map
mymap = plt.axes(projection=ccrs.PlateCarree())
mymap.add_feature(cfeature.LAND, alpha=0.5)
mymap.add_feature(cfeature.OCEAN, alpha=0.5)
mymap.add_feature(cfeature.BORDERS, linestyle='-', alpha=0.5)
mymap.add_feature(cfeature.RIVERS)
mymap.xaxis.set_visible(True)
mymap.yaxis.set_visible(True)
# Plots the data onto map
pcm = plt.scatter(Lon[index_eutroph>0], Lat[index_eutroph>0], c = index_eutroph.data[index_eutroph>0],
cmap = 'RdYlGn_r', transform=ccrs.PlateCarree())
plt.title("Eutrophication level on " + d, fontsize = 20);
plt.legend(*pcm.legend_elements(), title = "Index");
# Plot labels
plt.ylabel("Latitude", fontsize=14);
plt.xlabel("Longitude", fontsize=14);
if save_plot:
plt.savefig(out_dir + "/" + "ICEP_map_custom-legend")
# This plots the matrix with selected indices, as we try to recreate
# the colormap in Figure 6 from Sarelli et al (2018).
plt.figure(figsize=(20,14))
# Creates the map
mymap = plt.axes(projection=ccrs.PlateCarree())
mymap.add_feature(cfeature.LAND, alpha=0.5)
mymap.add_feature(cfeature.OCEAN, alpha=0.5)
mymap.add_feature(cfeature.BORDERS, linestyle='-', alpha=0.5)
mymap.add_feature(cfeature.RIVERS)
mymap.xaxis.set_visible(True)
mymap.yaxis.set_visible(True)
# Plots the data onto map
plt.scatter(Lon[(index_eutroph>0)&(index_eutroph<=2)], Lat[(index_eutroph>0)&(index_eutroph<=2)],
marker = 's', s= 10,
transform=ccrs.PlateCarree(),
c = 'limegreen')
plt.scatter(Lon[index_eutroph==3], Lat[index_eutroph==3],
marker = 's', s= 10,
transform=ccrs.PlateCarree(),
c = 'yellow')
plt.scatter(Lon[(index_eutroph>3)&(index_eutroph<=5)], Lat[(index_eutroph>3)&(index_eutroph<=5)],
marker = 's', s= 10,
transform=ccrs.PlateCarree(),
c = 'orange')
plt.scatter(Lon[index_eutroph>5], Lat[index_eutroph>5],
marker = 's', s= 10,
transform=ccrs.PlateCarree(),
c = 'red')
#Custom legend
custom_lines = [Line2D([0], [0], color='limegreen', lw=4),
Line2D([0], [0], color='yellow', lw=4),
Line2D([0], [0], color='orange', lw=4),
Line2D([0], [0], color='red', lw=4)]
plt.legend(custom_lines, ['Index <= 2: non-problem areas',
'Index = 3: tendency for eutrophication events',
'4 <= Index <=5: possibility of eutrophication events',
'Index = 6: problem areas'],
title = "Legend")
# Plot labels
plt.ylabel("Latitude", fontsize=14);
plt.xlabel("Longitude", fontsize=14);
if save_plot:
plt.savefig(out_dir + "/" + "ICEP_map_SarelliEtAl2018-legend")
I think we are very close!
import session_info
session_info.show()