1.6.1. Impact of using successive satellites on Ocean Colour time series#
Production Date: 16-03-2024
Produced by: Chiara Volta (ENEA, Italy) and Salvatore Marullo (CNR, Italy)
🌍 Use case: Assess potential decrease in marine fish stock due to oceanic desertification#
❓ Quality assessment question#
Does the consistency of chlorophyll-a and remote sensing reflectance get affected by the succession of satellite sensors over time?
The Ocean Colour dataset version 6.0, as produced for the Copernicus Climate Change Service (C3S), includes the variables, mass concentration of chlorophyll-a and remote sensing reflectance (Rrs) for six wavelengths from October 1997 to present [1].
Chlorophyll-a concentration is typically used as proxy for algal biomass and fisheries production [2] [3] and is derived through specific algorithms which rely on remote sensing reflectance (Rrs) [4].
The dataset combines measurements carried out by six satellite sensors [1]: SeaWiFS, MERIS, MODIS-Aqua, VIIRS, OLCI-3A and OLCI-3B. No sensor was operational over the whole temporal coverage of the datasets and only OLCI-3A and OLCI-3B are still operational.
Here, the goal is to assess if the stability of Ocean Colour time series is affected by the transition of different satellite sensors over time.
📢 Quality assessment statement#
These are the key outcomes of this assessment
The Ocean Colour dataset version 6.0:
provides a time consistent record of remote sensing reflectance at 443 and 560 nm, and derived chlorophyll-a concentrations, which is not affected by the succession of different satellite sensors. However, data prior to March 2002 should be carefully evaluated.
would overestimate chlorophyll-a concentrations before March 2002, especially when high-latitude, more productive oceanic regions are included in the analysis.
indicates decreasing chlorophyll-a trends in the more productive regions analysed, while the most oligotrophic area exhibits an increasing trend.
📋 Methodology#
This notebook provides an assessment of the ability of Ocean Colour dataset version 6.0 to describe the temporal variability of chlorophyll-a concentrations and Rrs at 443 and 560 nm over different oceanic regions.
The analysis and results are detailed in the sections below:
1. Choose the data to use and set up the code
Import required packages
Define parameters to be analysed (time period, variables, regions) and data request
2. Area average function and data retrieval
Area average function
Retrieve daily averages for the selected variables, time period and regions
3. Plot and statistics functions
Compute statistics and define plot function
Define satellite missions
4. Display and discuss results
Display results
Discussion
📈 Analysis and results#
1. Choose the data to use and set up the code#
Import required packages#
Besides the standard libraries used to manage and analyse multidimensional arrays, the C3S EQC custom functions c3s_eqc_automatic_quality_controlis imported to download data and calculate statistics.
Show code cell source
import matplotlib.patheffects as pe
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymannkendall as mk
import xarray as xr
from c3s_eqc_automatic_quality_control import diagnostics, download, utils
plt.style.use("seaborn-v0_8-notebook")
Define parameters to be analysed (time period, variables, regions) and data request#
The analysis performed in this notebook focuses on the time series of chlorophyll-a concentration and Rrs at 443 and 560 nm, the most commonly used bands to derive chlorophyll-a data, over a 25-year period (January 1998 - December 2022), in four regions: the global ocean, the oceanic region between 50°S and 50°N, the region where the North Atlantic SubTropical Gyre (NASTG) is located [5] and a 10°x10° area, within the South Pacific Gyre (SPG), over which the inter-sensor bias-correction process was validated [6].
Show code cell source
# Time period
start = "1998-01"
stop = "2022-12"
# Variables
variables = [
"chlor_a",
"Rrs_443",
"Rrs_560",
]
assert set(variables) <= {"chlor_a"} | {
f"Rrs_{wl}" for wl in (443, 560)
}
# Regions
regions = {
"Global Ocean": {
"lon_slice": slice(-180, 180),
"lat_slice": slice(90, -90),
},
"Latitudes 50$^o$S-50$^o$N": {
"lon_slice": slice(-180, 180),
"lat_slice": slice(50, -50),
},
"NASTG": {
"lon_slice": slice(-80, 0),
"lat_slice": slice(50, 0),
},
"SPG": {
"lon_slice": slice(-120, -110),
"lat_slice": slice(-20, -30),
},
}
for region, slices in regions.items():
# Enforce sorting as original data
for k, v in slices.items():
assert v.start >= v.stop if k == "lat_slice" else v.start <= v.stop, (region, k)
# Define data request
collection_id = "satellite-ocean-colour"
request = {
"projection": "regular_latitude_longitude_grid",
"version": "6_0",
"format": "zip",
}
2. Area average function and data retrieval#
Area average function#
The regionalised_spatial_weighted_mean function, which accounts for the varying surface area at different latitudes, is defined to calculate the area averages over the selected regions. Area averages for chlorophyll-a concentration are computed using log-transformed daily data, and then unlogged [7].
Show code cell source
def regionalised_spatial_weighted_mean(
ds, variable, lon_slice, lat_slice, log_transformed, vmin, vmax
):
da = ds[variable]
da = utils.regionalise(da, lon_slice=lon_slice, lat_slice=lat_slice)
if vmin is not None:
da = da.where(da > vmin)
if vmax is not None:
da = da.where(da < vmax)
if log_transformed:
with xr.set_options(keep_attrs=True):
da = 10 ** diagnostics.spatial_weighted_mean(np.log10(da))
da.attrs["long_name"] = da.attrs["long_name"].replace(
" (not log-transformed)", ""
)
else:
da = diagnostics.spatial_weighted_mean(da)
return da.to_dataset(name=variable)
Data retrieval#
The regionalised_spatial_weighted_mean is applied to the selected variables and regions, and the time series are downloaded as a single array. Chlorophyll-a concentrations outside the range 0.01-100 mg m-3 are excluded from the download and the analysis [8].
Show code cell source
datasets = []
for variable in variables:
for region, slices in regions.items():
requests = download.update_request_date(
request
| {
"variable": "remote_sensing_reflectance"
if variable.startswith("Rrs")
else "mass_concentration_of_chlorophyll_a"
},
start=start,
stop=stop,
stringify_dates=True,
)
ds = download.download_and_transform(
collection_id,
requests,
transform_func=regionalised_spatial_weighted_mean,
transform_func_kwargs=slices
| {
"variable": variable,
"vmin": 1.0e-2 if variable == "chlor_a" else None,
"vmax": 1.0e2 if variable == "chlor_a" else None,
"log_transformed": variable == "chlor_a",
},
chunks={"year": 1, "month": 1, "variable": 1},
quiet=True,
)
datasets.append(ds.expand_dims(region=[region]))
ds = xr.merge(datasets).compute()
3. Plot and statistics functions#
Compute statistics and define plot function#
Functions to calculate 48-month rolling mean and linear trend for the selected variables and regions are embedded in the plot function. Linear trend equations and their statistical significance (p) are specified in the legend of each figure. Time-averaged values are also computed for the selected variables and regions.
Show code cell source
def plot_timeseries(da):
fig, ax = plt.subplots(1, 1)
# Compute running mean
da_daily = da.resample(time="D").mean()
da_yearly = da.resample(time="12ME").mean()
da_running = (
da.resample(time="MS")
.mean()
.rolling(time=48, center=True, min_periods=1)
.mean()
)
# Compute linear trend
trend, h, p, z, tau, s, var_s, slope, intercept = mk.original_test(
da_yearly.squeeze()
)
# Plots lines
da_daily.plot(
label="daily", ms=4, color="tab:grey", ls=" ", marker=".", ax=ax, zorder=1
)
da_running.plot(
add_legend=False,
label="48-month running",
color="r",
ls="-",
marker=" ",
ax=ax,
path_effects=[pe.Stroke(linewidth=4, foreground="w"), pe.Normal()],
zorder=3,
)
ax.plot(
da_yearly["time"],
np.arange(da_yearly.sizes["time"]) * slope + intercept,
label=f'Linear regression:\ny={slope:.1e}t+{intercept:.1e}\n{p=:.1e}',
color="k",
ls="--",
marker=" ",
path_effects=[pe.Stroke(linewidth=4, foreground="w"), pe.Normal()],
zorder=3,
)
# Final settings
plt.legend(loc="upper left", bbox_to_anchor=(1, 1))
plt.grid()
return fig, ax
averages = ds.mean(dim='time')
Define satellite missions#
The name and data acquiring period are specified for each satellite sensor used to produce the dataset [1]. This information is displayed in the time series figures (see below).
Show code cell source
missions = {
'mission1': {
'mission': 'SeaWiFS',
'start': np.datetime64('1998-01-01'),
'stop': np.datetime64('2010-12-31'),
'linestyle': '-',
},
'mission2': {
'mission': 'MERIS',
'start': np.datetime64('2002-03-01'),
'stop': np.datetime64('2012-04-30'),
'linestyle': '-' ,
},
'mission3': {
'mission': 'MODIS-Aqua',
'start': np.datetime64('2002-05-01'),
'stop': np.datetime64('2019-12-31'),
'linestyle': '-' ,
},
'mission4': {
'mission': 'VIIRS',
'start': np.datetime64('2011-11-01'),
'stop': np.datetime64('2019-12-31'),
'linestyle': '-' ,
},
'mission5': {
'mission': 'OLCI-3A',
'start': np.datetime64('2016-05-01'),
'stop': np.datetime64('2022-12-31'),
'linestyle': '-' ,
},
'mission6': {
'mission': 'OLCI-3B',
'start': np.datetime64('2018-05-01'),
'stop': np.datetime64('2022-12-31'),
'linestyle': '-' ,
},
}
4. Display and discuss results#
Display results#
Time series of chlorophyll-a concentrations and Rrs at 443 and 560 nm are displayed for each region analysed, together with the duration of every satellite mission used to retrieve data. Linear trend equations and their statistical significance are indicated in the legends. Data are plotted on the same zoomed-in scale to focus on long-term chlorophyll-a and Rrs trends in the four regions analysed and facilitate their comparison.
Show code cell source
for variable in variables:
for region in regions:
fig, ax = plot_timeseries(ds[variable].sel(region=[region]))
if variable == "chlor_a":
ax.set_title(f"Chlorophyll-a concentration in seawater\n{region=:s}")
ax.set_ylabel("mg m-3")
if variable == "Rrs_443":
ax.set_title(f"Sea surface reflectance at 443 nm\n{region=:s}")
ax.set_ylabel(f"{ds[variable].attrs['units']}")
if variable == "Rrs_560":
ax.set_title(f"Sea surface reflectance at 560 nm\n{region=:s}")
ax.set_ylabel(f"{ds[variable].attrs['units']}")
ax.set_ylim(ds[variable].quantile([0.001, 0.999]).values.tolist())
for mission_name, mission in missions.items():
mission_start = mission['start']
mission_stop = mission['stop']
linestyle = mission['linestyle']
plt.axvline(mission_start, color='black', linestyle=linestyle, lw=1.0)
plt.axvline(mission_stop, color='black', linestyle=linestyle, lw=1.0)
plt.annotate('', xy=(np.datetime64('1998-01-01'), ds[variable].quantile(0.980).values), xytext=(np.datetime64('2010-12-31'), ds[variable].quantile(0.980).values), arrowprops=dict(arrowstyle='<->', color='black', lw=1.2))
plt.annotate('', xy=(np.datetime64('2002-03-01'), ds[variable].quantile(0.984).values), xytext=(np.datetime64('2012-04-30'), ds[variable].quantile(0.984).values), arrowprops=dict(arrowstyle='<->', color='black', lw=1.2))
plt.annotate('', xy=(np.datetime64('2002-05-01'), ds[variable].quantile(0.991).values), xytext=(np.datetime64('2019-12-31'), ds[variable].quantile(0.991).values), arrowprops=dict(arrowstyle='<->', color='black', lw=1.2))
plt.annotate('', xy=(np.datetime64('2011-11-01'), ds[variable].quantile(0.993).values), xytext=(np.datetime64('2019-12-31'), ds[variable].quantile(0.993).values), arrowprops=dict(arrowstyle='<->', color='black', lw=1.2))
plt.annotate('', xy=(np.datetime64('2016-05-01'), ds[variable].quantile(0.995).values), xytext=(np.datetime64('2022-12-31'), ds[variable].quantile(0.995).values), arrowprops=dict(arrowstyle='<->', color='black', lw=1.2))
plt.annotate('', xy=(np.datetime64('2018-05-01'), ds[variable].quantile(0.9975).values), xytext=(np.datetime64('2022-12-31'), ds[variable].quantile(0.9975).values), arrowprops=dict(arrowstyle='<->', color='black', lw=1.2))
plt.text(np.datetime64('1998-01-01'), ds[variable].quantile(0.987).values.tolist(), 'SeaWiFS', ha='left', va='center')
plt.text(np.datetime64('2002-05-01'), ds[variable].quantile(0.987).values.tolist(), 'MERIS', ha='left', va='center')
plt.text(np.datetime64('2002-05-01'), ds[variable].quantile(0.993).values.tolist(), 'MODIS-Aqua', ha='left', va='center')
plt.text(np.datetime64('2011-11-01'), ds[variable].quantile(0.995).values.tolist(), 'VIIRS', ha='left', va='center')
plt.text(np.datetime64('2016-05-01'), ds[variable].quantile(0.996).values.tolist(), 'OLCI-3A', ha='left', va='center')
plt.text(np.datetime64('2018-05-01'), ds[variable].quantile(0.998).values.tolist(), 'OLCI-3B', ha='left', va='center')
plt.show()
Discussion#
Results show that the Ocean Colour dataset version 6.0 captures both the seasonal and interannual variability of the variables analysed in the four selected regions.
Chlorophyll-a:
Chlorophyll-a timeseries in the global ocean and the region between 50°S and 50°N show similar seasonal and interannual variabilities, and trends. Long-term trends obtained through a 48-month rolling mean show a maxima around 2002 in both the global ocean and the oceanic region between latitudes 50°S and 50°N, followed by a sharp decrease until 2004. Afterwards, a progressive decrease until around 2015 and a slight increase until 2022 are observed.
The largest chlorophyll-a seasonal variability is observed in the NASTG’s region, due to the strong seasonal expansion/contraction cycle of the gyre [5], together with a rather constant long-term decreasing trend.
The lowest chlorophyll-a concentrations and seasonal variability are observed in the SPG’s region and indicate a high stability in the area [6]. Linear trends are statistically significant (p<0.05) in all regions analyzed here, indicating that chlorophyll-a is decreasing everywhere except in the SPG. These results align with the recent paper [9], where decreasing (increasing) chlorophyll-a is associated with increasing (decreasing) sea surface temperature trends.
Rrs:
Compared to chlorophyll-a, average Rrs at 443 nm shows reverse patterns and statistically significant (p<0.05) trends in all selected regions, which are consistent with results in [9]. It is worth noting that results for the SPG’s region are consistent with the post-bias correction validation [6] and show no suspicious artefacts along with the satellite succession.
Rrs at 560 nm trends are rather linear in all the regions analysed here, except the NASTG, where a progressive increase between 1998 and 2017, followed by a rapid decrease until 2022, is observed. Linear trends calculated here for Rrs at 560 nm are not statistically significant (p>0.05), although an increasing tendency is observed everywhere, except in the SPG where Rrs at 560 nm decreases over time.
Relatively high chlorophyll-a concentrations and low Rrs at 443 nm are observed when SeaWiFS was the only active satellite (i.e., from January 1998 to March 2002) in both the global ocean and the oceanic region between 50°S and 50°N, but not in the NASTG’s and SPG’s regions. Overall, these results suggest that, prior to March 2002, Rrs at 443 would be underestimated, especially when the analysis is extended to include high-latitude, more productive oceanic zones, and trigger chlorophyll-a overestimates. Despite the extensive time series analyzed here and the relatively short time series length required to detect climate-driven variations in low-latitude and/or oligotrophic regions, such as the NASTG and the SPG [9] [10], it is not yet possible to unequivocally associate the trends calculated here with climate change. An additional analysis, performed using the same notebook, but excluding SeaWiFS data prior to March 2002, indicates that they do not significantly modify the long-term trends computed for the variables analysed here nor their statistics in all selected regions.
ℹ️ If you want to know more#
Key resources#
CDS catalogue entry used in this notebook is the Ocean colour daily data from 1997 to present derived from satellite observations
Data download is from CDS
Product Documentation is available on the CDS website
Code libraries used:
C3S EQC custom functions,
c3s_eqc_automatic_quality_control, prepared by B-Open
Further readings:
Ocean Colour Monitoring Portals: CMEMS, ESA-GlobColour Project, ESA-OC-CCI, ESA-Sentinel 3, IMOS, NASA, NOAA CoastWatch
Groom, S., et al. (2019). Satellite Ocean Colour: Current Status and Future Perspective. Frontiers in Marine Science, 6.
References#
[1] Jackson, T., et al. (2023). C3S Ocean Colour Version 6.0: Product User Guide and Specification. Issue 1.1. E.U. Copernicus Climate Change Service. Document ref. WP2-FDDP-2022-04_C3S2-Lot3_PUGS-of-v6.0-OceanColour-product.
[2] Behrenfeld, M.J., & Falkowski, P.G. (1997). A consumer’s guide to phytoplankton primary productivity models. Limnology and Oceanography, 42, 1479-1491. 9.
[3] Friedland, K.D., et al. (2012). Pathways between Primary Production and Fisheries Yields of Large Marine Ecosystems. PLoS ONE, 7(1), e28945.
[4] Jackson, T., et al. (2022). C3S Ocean Colour Version 6.0: Algorithm Theoretical Basis Document. Issue 1.1. E.U. Copernicus Climate Change Service. Document ref. WP2-FDDP-2022-04_C3S2-Lot3_ATBD-of-v6.0-OceanColour-product.
[5] Leonelli, F.E., et al. (2022). Ultra-oligotrophic waters expansion in the North Atlantic Subtropical Gyre revealed by 21 years of satellite observations. Geophysical Research Letters, 49, e2021GL096965.
[6] Jackson, T., et al. (2023) C3S Ocean Colour Version 6.0: Product Quality Assessment Report. Issue 1.1. E.U. Copernicus Climate Change Service. Document ref. WP2-FDDP-2022-04_C3S2-Lot3_ PQAR-of-v6.0-OceanColour-product.
[7] Campbell, J.W. (1995). The lognormal distribution as a model for bio-optical variability in the sea. Journal of Geophysical Research, 100, 13237-13254.
[8] Sathyendranath, S., et al. (2019). An Ocean-Colour time series for use in climate studies: the experience of the Ocean-Colour Climate Change Initiative (OC-CCI). Sensors, 19, 4285.
[9] Cael, B.B., et al. (2023). Global climate-change trends detected in indicators of ocean ecology. Nature, 619, 551-554.
[10] Henson, S.A., et al. (2010). Detection of anthropogenic climate change in satellite records of ocean chlorophyll and productivity. Biogeosciences, 7, 621-640.