Subsetting PACE OCI data

Subsetting PACE OCI data#

Authors: Anna Windle (NASA, SSAI), Carina Poulin (NASA, SSAI), Ian Carroll (NASA, UMBC)

Last Updated: July 1, 2026

The following notebooks are prerequisites for this tutorial.

Learn with OCI: Data Access

An Earthdata Login account is required to access data from the NASA Earthdata system, including NASA ocean color data.

Summary#

This tutorial demonstrates how to subset PACE OCI L2 and L3M data. L2 data is subsetted here using harmony-py, while L3M data is subsetted using xarray.

Harmony is a web service that allows you to customize many NASA datasets, including the ability to subset, reproject, and reformat files. Data can be subsetted for a geographic region, a temporal range and by variable. In some cases, data can be “reprojected” from its native coordinate reference system (CRS) to the coordinate reference system relevant to your analysis. Data can also be reformatted from its native file format to a format that is more relevant for your application. These services are collectively called transformation services. However, not all services and transformations are available for all datasets. You will learn how to discover which services are available for a given dataset.

Harmony services can be used in multiple ways:

through a graphical user interface (GUI) while downloading applicable granules from Earthdata Search,
by direct requests to Harmony’s RESTful API, or, as in this tutorial,
using the harmony-py Python package.

The Python package handles NASA Earthdata Login (EDL) authentication and optionally integrates with the CMR Python Wrapper by accepting collection results as a request parameter. It’s convenient for scientists who wish to use Harmony from Jupyter notebooks. After this tutorial, you can dive deeper into harmony-py on ReadTheDocs.

This tutorial demonstrates how to subset and reformat PACE OCI data from the NASA Earthdata Cloud using harmony-py. It is modelled off of tutorials developed by NSIDC and NASA-Openscapes.

Learning Objectives#

At the end of this notebook you will know how to:

Use the harmony-py Python library to spatially and temporally subset PACE OCI Level-2 data
Download the subsetted data to your local directory
Stream the subsetted data
Open and plot subsetted L2 PACE OCI data
Open and subset PACE OCI Level-3 mapped data using xarray

1. Setup#

Begin by importing all of the packages used in this notebook.

from datetime import datetime
from pathlib import Path

import earthaccess
import matplotlib.pyplot as plt
import numpy as np
import rasterio
import xarray as xr
from harmony import (
    BBox,
    CapabilitiesRequest,
    Client,
    Collection,
    JobsRequest,
    LinkType,
    Request,
)
from IPython.display import JSON
from rasterio.enums import Resampling

2. Earthdata authentication and Harmony client initalization#

The earthaccess package lets us easily authenticate with NASA’s Earthdata Login (EDL) and pass an EDL token to Harmony.

auth = earthaccess.login()
harmony_client = Client(token=earthaccess.get_edl_token())

3. Discover subsetting capabilities for Level-2 PACE OCI data#

Define which data set you’d like to access using either the dataset short_name or its collection concept_id. For this example, we’ll use PACE OCI Level-2 Regional Ocean Biogeochemical Properties Data (PACE_OCI_L2_BGC).

capabilities_request = CapabilitiesRequest(short_name="PACE_OCI_L2_BGC")
capabilities = harmony_client.submit(capabilities_request)
capabilities

{'conceptId': 'C4124887072-OB_CLOUD',
 'shortName': 'PACE_OCI_L2_BGC',
 'variableSubset': True,
 'bboxSubset': True,
 'shapeSubset': True,
 'temporalSubset': True,
 'concatenate': False,
 'reproject': False,
 'outputFormats': ['application/netcdf'],
 'services': [{'name': 'podaac/l2-subsetter',
   'href': 'https://cmr.earthdata.nasa.gov/search/concepts/S1962070864-POCLOUD',
   'capabilities': {'subsetting': {'temporal': True,
     'bbox': True,
     'variable': True,
     'shape': True,
     'dimension': True},
    'output_formats': ['application/netcdf']}}],
 'variables': [{'name': 'geophysical_data/carbon_phyto',
   'href': 'https://cmr.earthdata.nasa.gov/search/concepts/V4198030748-OB_CLOUD'},
  {'name': 'geophysical_data/carbon_phyto_unc',
   'href': 'https://cmr.earthdata.nasa.gov/search/concepts/V4198030761-OB_CLOUD'},
  {'name': 'geophysical_data/chlor_a',
   'href': 'https://cmr.earthdata.nasa.gov/search/concepts/V4198030774-OB_CLOUD'},
  {'name': 'geophysical_data/chlor_a_unc',
   'href': 'https://cmr.earthdata.nasa.gov/search/concepts/V4198030781-OB_CLOUD'},
  {'name': 'geophysical_data/pic',
   'href': 'https://cmr.earthdata.nasa.gov/search/concepts/V4198030803-OB_CLOUD'},
  {'name': 'geophysical_data/poc',
   'href': 'https://cmr.earthdata.nasa.gov/search/concepts/V4198030812-OB_CLOUD'}],
 'capabilitiesVersion': '2'}

You can see here under ['services'] that this dataset can be subsetted using the 'podacc/l2-subsetter'

4. Build and submit a request#

The Harmony service is based around user-submitted jobs. Using the service is a two-step process of (1) creating a specific type of request, and (2) passing the request to the harmony_client in exchange for some response from Harmony.

One type of request, the JobsRequest, let’s you check for any previously submitted jobs. You can filter these jobs using labels, which we will learn how to apply below.

First step, create a request:

request = JobsRequest(labels=["help-hub-tutorial"])

Second, get a response from your request submission:

response = harmony_client.submit(request)
JSON(response)

<IPython.core.display.JSON object>

For your first time through this tutorial, you shouldn’t see any existing jobs. If you’ve already submitted a job with this label (e.g. because you are re-running this tutorial), then the response includes information about that existing job.

Let’s continue to build a request for the subsetting job we want the service to run. Using the “labels” keyword tags jobs so that we can easily find and re-use results from this job later.

request = Request(
    collection=Collection(id="PACE_OCI_L2_BGC"),
    spatial=BBox(-76.75, 36.97, -75.74, 39.01),
    temporal={"start": datetime(2025, 7, 1), "stop": datetime(2025, 8, 1)},
    variables=["geophysical_data/chlor_a"],
    labels=request.labels,
)

Depending on whether we got a non-zero count from the JobsRequest, we can either submit the subsetting request or dig into the job info we retrieved. Either way, we come up with a job_id that uniquely identifies the subsetting job.

if response["count"]:
    job_id = response["jobs"][0]["jobID"]
else:
    job_id = harmony_client.submit(request)
job_id

'a954b3a4-5689-4d63-85df-a5797ca1bd0b'

The job_id can be used to check on the details of your job progress. If you are logged into Earthdata, you can see your jobs here: https://harmony.earthdata.nasa.gov/workflow-ui or look at your process by running this line:

harmony_client.wait_for_processing(job_id, show_progress=True)

 [ Processing:   0% ] |                                                   | [/]
 [ Processing: 100% ] |###################################################| [|]

Note: Requesting a job can take a variable amount of time. In our experience, running this particular job has ranged from 2 to 24 minutes.

job_summary = harmony_client.result_json(job_id)

print("Number of granules:", job_summary["numInputGranules"])
print("Original Data Size:", job_summary["originalDataSize"])
print("Output Data Size:", job_summary["outputDataSize"])
print("Data Size % Change:", job_summary["dataSizePercentChange"])

Number of granules: 44
Original Data Size: 759.11 MiB
Output Data Size: 5.06 MiB
Data Size % Change: 99.33% reduction

If you want to access a specific job that has already run, you can simply assign a known id to job_id and continue below.

Results are kept for 30 days in the Harmony S3 bucket.

6. Plot the subsetted data#

Let’s do a quick plot of chlor_a from the first granule:

plot = ds["chlor_a"].plot()

../../_images/8e885f39c9082dff14d3091ea171d2166486e4c825becd4a17dbd8ad6564fc30.png

Now let’s plot multiple subsetted granules:

fig, axes = plt.subplots(2, 5, figsize=(10, 4), constrained_layout=True)
axes = axes.ravel()

for ax, file in zip(axes, urls[:10]):

    dt = xr.open_datatree(file, **kwargs)
    ds = xr.merge(dt.to_dict().values())
    ds = ds.set_coords(("longitude", "latitude"))
    date = ds.attrs["time_coverage_start"]
    im = ds["chlor_a"].plot(ax=ax, cmap="viridis", add_colorbar=False, vmin=0, vmax=20)
    ax.set_title(date, fontsize=8)

fig.colorbar(im, ax=axes, orientation="vertical", shrink=0.8, label="Chl a (mg m-3)")
plt.show()

../../_images/8b2b42bd50706513d7b9ca027848b43997987846b9aec39083038461041a2a06.png

To plot using lat, lon coordinates, we need to project the data onto a defined grid with a given resolution. We will use code presented in the Projecting PACE Data onto a Predefined Grid tutorial.

def grid_data(src, resolution, dst_crs="epsg:4326", resampling=Resampling.nearest):
    """
    Reproject a L2 dataset to match an input grid. Makes sure 3D variables are
        in (Z, Y, X) dimension order, and all variables have spatial dims/crs 
        assigned.
    Args:
        src - an xarray dataset or dataarray to reproject
        resolution - resolution of the output grid, in dst_crs units
        dst_crs - CRS of the output data
        resampling - resampling method (see rasterio.enums)
    Returns:
        dst - projected xr dataset
    """
    if (len(list(src.dims)) == 3) and (list(src.dims)[0] != "wavelength_3d"):
        src = src.transpose("wavelength_3d", ...)
    src = src.rio.set_spatial_dims("pixels_per_line", "number_of_lines")
    src = src.rio.write_crs("epsg:4326")

    # Calculating the default affine transform
    defaults = rasterio.warp.calculate_default_transform(
        src.rio.crs,
        dst_crs,
        src.rio.width,
        src.rio.height,
        left=src.attrs["geospatial_lon_min"],
        bottom=src.attrs["geospatial_lat_min"],
        right=src.attrs["geospatial_lon_max"],
        top=src.attrs["geospatial_lat_max"],
    )
    # Aligning that transform to our desired resolution
    transform, width, height = rasterio.warp.aligned_target(*defaults, resolution)
    
    # Run projection
    dst = src.rio.reproject(
        dst_crs=dst_crs,
        shape=(height, width),
        transform=transform,
        src_geoloc_array=(
            src["longitude"],
            src["latitude"],
        ),
        nodata=np.nan,
        resample=resampling,
    )
    dst["x"] = dst["x"].round(9)
    dst["y"] = dst["y"].round(9)
    
    return dst.rename({"x": "longitude", "y": "latitude"})

We choose a 0.015 degree resolution, and the function above employs the plate carrée (lat, lon) projection.

resolution = (0.015, 0.015)

dt = xr.open_datatree(urls[0], **kwargs)
ds = xr.merge(dt.to_dict().values())
ds = ds.set_coords(("longitude", "latitude"))
ds_gridded = grid_data(ds, resolution)
plot = ds_gridded["chlor_a"].plot()

../../_images/2f8463117d10e217131c779df7737edac828f43ff738839c7fa762f20a03c5dc.png

Plotting the first 10 files as subplots, and keeping the gridded “chlor_a” data array for the next section:

fig, axes = plt.subplots(2, 5, figsize=(10, 4), constrained_layout=True)
axes = axes.ravel()
da = []

for ax, file in zip(axes, urls[:10]):

    dt = xr.open_datatree(file, **kwargs)
    ds = xr.merge(dt.to_dict().values())
    ds = ds.set_coords(("longitude", "latitude"))

    ds_gridded = grid_data(ds, resolution)
    da.append(ds_gridded["chlor_a"])
    date = ds_gridded.attrs["time_coverage_start"]
    im = ds_gridded["chlor_a"].plot(ax=ax, cmap="viridis", add_colorbar=False, vmin=0, vmax=20)
    ax.set_title(date, fontsize=8)

    ax.set_xlabel("")
    ax.set_ylabel("")

fig.colorbar(im, ax=axes, orientation="vertical", shrink=0.8, label="Chl a (mg m-3)")
plt.show()

../../_images/25807084000bd61f875b0c680045dc45e89b6bffba24f6570fc118c3765a5259.png

Now, we can make a 10-day Chl a composite:

da = xr.concat(da, dim="scene")
chlor_a_mean = da.mean(dim="scene")
plot = chlor_a_mean.plot(cmap="viridis", vmin=0, vmax=20)

../../_images/2dbfc11b52cdd8932c6305af514f7f8ce4375b36289e6f3b97d838546bb3c267.png

Subsetting PACE OCI data

Contents

Subsetting PACE OCI data#

Summary#

Learning Objectives#

1. Setup#

2. Earthdata authentication and Harmony client initalization#

3. Discover subsetting capabilities for Level-2 PACE OCI data#

4. Build and submit a request#

5. Access the subsetted data#

Download a single file#

Download all files#

Stream the files#

6. Plot the subsetted data#

7. Subsetting L3M data#