2. Comparison of chlorophyll data from different sensors
Different ocean color sensors have been launched since 1997 to provide continuous global ocean color data. Unfortunately, because of differences in sensor design and calibration, chlorophyll-a concentration values don’t match during their periods of overlap, making it challenging to study long-term trends.
As an example, we are going to plot time-series of mean chlorophyll a concentration from various sensors from 1997 to 2018 to look at the periods of overlap.
We are going to download data from Seawifs (1997-2010), MODIS (2002-present) and VIIRS (2012-present) and compare it to the ESA-CCI data which combines all 3 sensors into a homogeneous time-series.
First, let's load all the packages needed:
import urllib.request
import xarray as xr
import netCDF4 as nc
import pandas as pd
import numpy as np
import pyproj
from matplotlib import pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
np.warnings.filterwarnings('ignore')The OceanWatch website has a data catalog containing documentation and links to all the datasets available:
https://oceanwatch.pifsc.noaa.gov/doc.html
Navigate to the "Ocean Color" tab. From there you can access the different datasets using ERDDAP or THREDDS.
Go to ERDDAP to find the name of the dataset for monthly SeaWIFS data:
sw_chla_monthly_2018_0You should always examine the dataset in ERDDAP to check the date range, names of the variables and dataset ID, to make sure your
griddap calls are correct: https://oceanwatch.pifsc.noaa.gov/erddap/griddap/sw_chla_monthly_2018_0.htmlNotice also that for this dataset and others, the latitudes are ordered from North to South, which will affect the construction of the download URL. (ie. instead of selecting latitudes 0-40N, you need to request 40-0).
Test, test
- let's download data for a box around the Hawaiian Islands:
url='[https://oceanwatch.pifsc.noaa.gov/erddap/griddap/sw_chla_monthly_2018_0.nc?chlor_a[(1997-10-16T12:00:00Z):1:(2010-10-16T12:00:00Z)][(25):1:(15)][(198):1:(208)]'
urllib.request.urlretrieve(url, "sw.nc")`- let's use xarray to extract the data from the downloaded file:
sw_ds = xr.open_dataset('sw.nc',decode_cf=False)sw_ds.data_varsData variables:
chlor_a (time, latitude, longitude) float32 ...The downloaded data contains only one variable: chlor_a.
- let's compute the monthly mean over the region and extract the dates corresponding to each month of data:
swAVG=np.mean(sw_ds.chlor_a,axis=(1,2))
swdates=nc.num2date(sw_ds.time,sw_ds.time.units)sw_ds.close()url2="https://oceanwatch.pifsc.noaa.gov/erddap/griddap/aqua_chla_monthly_2018_0.nc?chlor_a[(2002-07-16T12:00:00Z):1:(2019-12-16T12:00:00Z)][(25):1:(15)][(198):1:(208)]"
urllib.request.urlretrieve(url2, "aq.nc")aq_ds = xr.open_dataset('aq.nc',decode_cf=False)
aqAVG=np.mean(aq_ds.chlor_a,axis=(1,2))
aqdates=nc.num2date(aq_ds.time,aq_ds.time.units)aq_ds.close()url3='https://oceanwatch.pifsc.noaa.gov/erddap/griddap/noaa_snpp_chla_monthly.nc?chlor_a[(2012-01-02T12:00:00Z):1:(2019-12-01T12:00:00Z)][(25):1:(15)][(198):1:(208)]'
urllib.request.urlretrieve(url3, "snpp.nc")snpp_ds = xr.open_dataset('snpp.nc',decode_cf=False)
snppAVG=np.mean(snpp_ds.chlor_a,axis=(1,2))
snppdates=nc.num2date(snpp_ds.time,snpp_ds.time.units)snpp_ds.close()url4="https://oceanwatch.pifsc.noaa.gov/erddap/griddap/esa-cci-chla-monthly-v4-2.nc?chlor_a[(1997-09-04):1:(2019-12-01T00:00:00Z)][(25):1:(15)][(198):1:(208)]"
urllib.request.urlretrieve(url4, "cci.nc")`cci_ds = xr.open_dataset('cci.nc',decode_cf=False)
cciAVG=np.mean(cci_ds.chlor_a,axis=(1,2))
ccidates=nc.num2date(cci_ds.time,cci_ds.time.units)
cci_ds.close()plt.figure(figsize=(12,5)) plt.plot(swdates,swAVG,label='sw',c='red',marker='.',linestyle='-') plt.plot(aqdates,aqAVG,label='aq',c='blue',marker='.',linestyle='-') plt.plot(snppdates,snppAVG,label='snpp',c='green',marker='.',linestyle='-')
plt.ylabel('Chl-a (mg/m^3)')
plt.legend()
plt.figure(figsize=(12,5))
plt.plot(ccidates,cciAVG, label='cci',c='black')
plt.scatter(swdates,swAVG,label='sw',c='red')
plt.scatter(aqdates,aqAVG,label='aq',c='blue')
plt.scatter(snppdates,snppAVG,label='snpp',c='green')
plt.ylabel('Chl-a (mg/m^3)')
plt.legend()