Applying Spectral Derivative Pigment (SDP) Phytoplankton Community Composition Algorithm to OCI data

Applying Spectral Derivative Pigment (SDP) Phytoplankton Community Composition Algorithm to OCI data#

Author(s): Anna Windle (NASA, SSAI), Ian Carroll (NASA, UMBC)
Adapted from code developed by: Max Danenhower (Bowdoin College), Sasha Kramer (Boston University)

Last updated: July 21, 2026

Summary#

This notebook applies the inversion algorithm described in Kramer et al., 2022 to estimate phytoplankton pigment concentrations from PACE OCI Rrs data. This algorithm, called Spectral Derivative Pigment (SDP), is currently being implemented in OBPG’s OCSSW software. This work was originally developed in MatLab by Sasha Kramer and subsequently translated to Python by Max Danenhower. This tutorial demonstrates how to apply the Python SDP algorithm to Level-2 (L2) PACE OCI data.

Learning Objectives#

At the end of this notebook you will know:

  • About searching for auxiliary data for the SDP algorithm

  • How to run the SDP algorithm on PACE OCI L2 data

1. Setup#

The SDP algorithm has been implemented as a Python package, so it can be easily installed, imported, and reused. While the package is not on PyPI or conda-forge, it can be installed directly from the source repository on GitHub. If you have followed the setup instructions, then sdp is available to import along with the other packages needed for this notebook.

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import earthaccess
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from sdp import sdp_from_pace

Set (and persist to your home directory on the host, if needed) your Earthdata Login credentials.

auth = earthaccess.login()

Assign “global” variables, which could be anything you want to define once and use consistently.

crs = ccrs.PlateCarree()

2. Access and open data#

We need the following data to run SDP:

  • Remote sensing reflectance (Rrs): PACE OCI L2 AOP data products

  • Sea surface salinity (SSS): JPL SMAP-SSS V5.0 CAP, 8-day running mean, level 3 mapped product from the NASA Soil Moisture Active Passive (SMAP) observatory

  • Sea surface temperature (SST): Group for High Resolution Sea Surface Temperature (GHRSST) Level 4 sea surface temperature

We can use earthaccess to find PACE OCI L2 data, checking in both the near-real-time and refined collections.

tspan = ("2026-05-05 17:35", "2026-05-05 17:35")
results = earthaccess.search_data(
    short_name=["PACE_OCI_L2_AOP_NRT", "PACE_OCI_L2_AOP"],
    temporal=tspan,
    count=1,
)
paths = earthaccess.open(results)

Let’s take a quick look at Rrs at 500nm:

datatree = xr.open_datatree(paths[-1])
rrs = datatree["geophysical_data"]["Rrs"]
for item in ("longitude", "latitude"):
    rrs[item] = datatree["navigation_data"][item]
rrs

<xarray.DataArray 'Rrs' (number_of_lines: 1710, pixels_per_line: 1272,
                         wavelength: 172)> Size: 1GB
[374120640 values with dtype=float32]
Coordinates:
    longitude   (number_of_lines, pixels_per_line) float32 9MB ...
    latitude    (number_of_lines, pixels_per_line) float32 9MB ...
  * wavelength  (wavelength) float32 688B 346.0 348.5 350.9 ... 716.8 719.3
Dimensions without coordinates: number_of_lines, pixels_per_line
Attributes:
    long_name:      Remote sensing reflectance
    units:          sr^-1
    standard_name:  surface_ratio_of_upwelling_radiance_emerging_from_sea_wat...
    valid_min:      -30000
    valid_max:      25000
fig, ax = plt.subplots(figsize=(8, 6), subplot_kw={"projection": crs})

da = rrs.sel({"wavelength": 500}, method="nearest")
da.plot(
    x="longitude",
    y="latitude",
    vmin=0,
    vmax=0.008,
    ax=ax,
    cbar_kwargs={"label": "Rrs ($sr^{-1}$)"},
)

ax.coastlines(resolution="10m")
ax.add_feature(cfeature.BORDERS, linestyle=":")
ax.gridlines(
    draw_labels=["left", "bottom"],
    linewidth=0.5,
    color="gray",
    alpha=0.5,
    linestyle="--",
)
plt.show()
../../_images/235a69635e833adb214d50d83c8b79fae31dd61dd7f27c1dde652ab76993edf1.png

Now, we can use sdp_from_pace to calculate phytoplankton pigment concentrations for this data. Let’s take a look at the docstring for this function:

?sdp_from_pace

You can see that this function accepts a bbox parameter for limiting calculations to a bounding box. The default is bbox=None, meaning the algorithm is applied to every single pixel in the L2 granule, which can take a significnat amount of time and may exceed available system memory. We will supply the bbox parameter with the following coordinates: 38 N, 35 S, -70 E, -67 W.

bbox = (-73.5, 37.5, -67, 40.5)

lon_min, lat_min, lon_max, lat_max = bbox
rect_lon = [lon_min, lon_max, lon_max, lon_min, lon_min]
rect_lat = [lat_min, lat_min, lat_max, lat_max, lat_min]
ax.plot(rect_lon, rect_lat, color="red", linewidth=2)

fig
../../_images/bcac7453c324e72fe7da3e165b8a2718b2807d09e09c36e715dc6c231045b6af.png

We need to find corresponding sea-surface salinity and temperature data for the SDP algorithm.

results = earthaccess.search_data(
    short_name="SMAP_JPL_L3_SSS_CAP_8DAY-RUNNINGMEAN_V5",
    temporal=tspan,
    count=1,
)
sss_paths = earthaccess.open(results)

results = earthaccess.search_data(
    short_name="MUR-JPL-L4-GLOB-v4.1",
    temporal=tspan,
    count=1,
)
sst_paths = earthaccess.open(results)

3. Run SDP on L2 Data#

Now we’re ready to run SDP. Let’s name what we want the outfile to be:

output_file = paths[-1].full_name.split("/")[-1]
output_file = output_file.replace(".nc", "_SDP_pigments.nc")
output_file

'PACE_OCI.20260505T173406.L2.OC_AOP.V3_2.NRT_SDP_pigments.nc'

And let’s run it! There will be some RuntimeWarnings that need to be considered.

sdp_from_pace(
    paths[-1],
    output_file,
    sss_file=sss_paths[-1],
    sst_file=sst_paths[-1],
    bbox=bbox,
)

starting interpolation
north 40.5 south 37.5 east -67 west -73.5
interpolation complete: 3.9881784060271457
rrs_interp shape: (110331, 374)
(70338,) (70338,) (70338, 301)
0
calculating residuals
ray availble resources {'memory': 1828832052.0, 'node:__internal_head__': 1.0, 'object_store_memory': 783785164.0, 'CPU': 3.0, 'node:192.168.18.110': 1.0}
residuals calculated 569.2999944640324
calculating 2nd derivative
2nd derivative calculated 0.1553171779960394
(70338, 299)
running coefs
coefs run complete 6.139696105034091
sdp calculation complete

/srv/conda/envs/notebook/lib/python3.13/site-packages/sdp/Kramer_hyperRrs.py:334: RuntimeWarning: overflow encountered in power
  bbp = (443 / wavelengths.reshape(-1, 1)) ** bbp_s
(run_batch pid=17673) /srv/conda/envs/notebook/lib/python3.13/site-packages/sdp/Kramer_hyperRrs.py:200: RuntimeWarning: invalid value encountered in power
(run_batch pid=17673)   aph = A * IOPs[0]**B
(run_batch pid=17674) /srv/conda/envs/notebook/lib/python3.13/site-packages/sdp/Kramer_hyperRrs.py:203: RuntimeWarning: invalid value encountered in divide
(run_batch pid=17674)   x = bb / (a + bb)
(run_batch pid=17675) /srv/conda/envs/notebook/lib/python3.13/site-packages/sdp/Kramer_hyperRrs.py:200: RuntimeWarning: invalid value encountered in power
(run_batch pid=17675)   aph = A * IOPs[0]**B
(run_batch pid=17674) /srv/conda/envs/notebook/lib/python3.13/site-packages/sdp/Kramer_hyperRrs.py:200: RuntimeWarning: invalid value encountered in power
(run_batch pid=17674)   aph = A * IOPs[0]**B
/srv/conda/envs/notebook/lib/python3.13/site-packages/sdp/Kramer_hyperRrs.py:405: RuntimeWarning: invalid value encountered in divide
  rrsP = bb / (a + bb)

4. Open and Plot SDP output#

Let’s open the new file:

ds = xr.open_dataset(output_file)
ds

<xarray.Dataset> Size: 12MB
Dimensions:          (number_of_lines: 369, pixels_per_line: 299)
Coordinates:
  * number_of_lines  (number_of_lines) int64 3kB 0 1 2 3 4 ... 365 366 367 368
  * pixels_per_line  (pixels_per_line) int64 2kB 0 1 2 3 4 ... 295 296 297 298
    longitude        (number_of_lines, pixels_per_line) float32 441kB ...
    latitude         (number_of_lines, pixels_per_line) float32 441kB ...
Data variables: (12/13)
    chla             (number_of_lines, pixels_per_line) float64 883kB ...
    chlb             (number_of_lines, pixels_per_line) float64 883kB ...
    chlc             (number_of_lines, pixels_per_line) float64 883kB ...
    zea              (number_of_lines, pixels_per_line) float64 883kB ...
    dvchla           (number_of_lines, pixels_per_line) float64 883kB ...
    butfuco          (number_of_lines, pixels_per_line) float64 883kB ...
    ...               ...
    allo             (number_of_lines, pixels_per_line) float64 883kB ...
    neo              (number_of_lines, pixels_per_line) float64 883kB ...
    viola            (number_of_lines, pixels_per_line) float64 883kB ...
    fuco             (number_of_lines, pixels_per_line) float64 883kB ...
    chlc3            (number_of_lines, pixels_per_line) float64 883kB ...
    perid            (number_of_lines, pixels_per_line) float64 883kB ...

You can see all the different pigment concentrations derived for our L2 granule. Let’s plot them:

fig, axs = plt.subplots(
    nrows=7,
    ncols=2,
    figsize=(10, 12),
    sharex=True,
    sharey=True,
    constrained_layout=False,
)

variables = tuple(ds)
axs = axs.ravel().tolist()
ax = axs.pop(1)
ax.set_visible(False)

for i, ax in enumerate(axs):
    da = ds[variables[i]]
    da.plot(
        x="longitude",
        y="latitude",
        cmap="viridis",
        shading="auto",
        vmax=da.quantile(0.99),
        ax=ax,
        cbar_kwargs={"label": r"$\mathrm{%s\ [mg\ m^{-3}]}$" % variables[i]},
    )
    ax.set_xlim(bbox[0], bbox[2])
    ax.set_ylim(bbox[1], bbox[3])
    ax.set_xlabel("")
    ax.set_ylabel("")

plt.tight_layout()
plt.show()
../../_images/a566276a7abb51b42fc08f578bdef0745789c7134e866f0ce520edeb15d7c4f9.png