File Structure at Three Processing Levels for the Ocean Color Instrument (OCI)

File Structure at Three Processing Levels for the Ocean Color Instrument (OCI)#

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

Last updated: August 3, 2025

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#

In this example we will use the earthaccess package to access OCI Level-1B (L1B), Level-2 (L2), and Level-3-Mapped (L3M) NetCDF files and open them using xarray. While you can learn alot exploring the datasets in this way, be ready to refer to the full dataset documentation for details about the data products and processing.

NetCDF (Network Common Data Format) is a binary file format for storing multidimensional scientific data (variables). It is optimized for array-oriented data access and support a machine-independent format for representing scientific data. Files ending in .nc are NetCDF files.

XArray is a package that supports the use of multi-dimensional arrays in Python. It is widely used to handle Earth observation data, which often involves multiple dimensions — for instance, longitude, latitude, time, and channels/bands.

Learning Objectives#

At the end of this notebok you will know:

How to find groups in a NetCDF file
How to use xarray to open OCI data
What key variables are present in the groups within OCI L1B, L2, and L3M files

1. Setup#

Begin by importing all of the packages used in this notebook. If your kernel uses an environment defined following the guidance on the tutorials page, then the imports will be successful.

import cartopy.crs as ccrs
import earthaccess
import h5netcdf
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr

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

auth = earthaccess.login(persist=True)

4. Explore L3M File Structure#

At L3M there are binned (B) and mapped (M) products available. The L3M remote sensing reflectance (Rrs) files contain global maps of the full spectra, so they can be big! We’ll use the same earthaccess method to find the data.

results = earthaccess.search_data(
    short_name="PACE_OCI_L3M_RRS",
    temporal=tspan,
)
len(results)

The L3M collections contain granules at different temporal and spatial resolutions. Let’s take a closer look at one of the larger sized granules first, using the metadata available on each result.

for item in results:
    if item.size() > 2000:
        display(item)
        break

Data: PACE_OCI.20240321_20240620.L3m.SNSP.RRS.V3_1.Rrs.4km.nc

Size: 3150.73 MB

Cloud Hosted: True

Notice some special parts in the granule name:

Right after L3m you have a period indicator. The value DAY means a daily aggregate, 8D means an 8-day aggregate, MO means a monthly aggregate, and SN(SP|SU|AU|WI) indicates one of four seasonal aggregates.
Right after chlor_a you have a spatial resolution marker, as horizontal cell size at the equator given either in km or deg. Read p as a decimal point.

Within a specific collection, earthaccess can retrieve granules whose name matches a pattern given in the granule_name argument. Let’s use that to narrow our search to values aggregated to a month and 0.1 degree resolutions.

results = earthaccess.search_data(
    short_name="PACE_OCI_L3M_RRS",
    granule_name="*.MO.*.4km.*",
    temporal=tspan,
)
len(results)

paths = earthaccess.open(results)

L3M data do not have any groups, so we can open as a dataset rather than a datatree.

Notice that this L3M granule has wavelength, lat, and lon under “Indexes”, so it’s easy to slice out a bounding box and map the “Rrs” variable at a given wavelength.

rrs = dataset["Rrs"]
rrs_440_bbox = rrs.sel({
    "wavelength": 440,
    "lat": slice(bbox[3], bbox[1]),
    "lon": slice(bbox[0], bbox[2]),
})
rrs_440_bbox

<xarray.DataArray 'Rrs' (lat: 49, lon: 24)> Size: 5kB
[1176 values with dtype=float32]
Coordinates:
  * lat         (lat) float32 196B 38.98 38.94 38.9 38.85 ... 37.06 37.02 36.98
  * lon         (lon) float32 96B -76.73 -76.69 -76.65 ... -75.85 -75.81 -75.77
    wavelength  float64 8B 440.0
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
    display_scale:  linear
    display_min:    0.0
    display_max:    0.025

fig, ax = plt.subplots(subplot_kw={"projection": ccrs.PlateCarree()})
im = rrs_440_bbox.plot(x="lon", y="lat", cmap="viridis", robust=True, ax=ax)
ax.gridlines(draw_labels={"left": "y", "bottom": "x"})
ax.coastlines()
plt.show()

../../_images/feb7450078aff9eb19caf14a9bcfa744d76098b67c019bf9ccb74213c98ad3a2.png

Also becuase the L3M variables have gridded lat and lon coordinates, it’s possible to stack multiple granules along a new dimension that corresponds to time. Using the granule_name filter, let’s get low spatial-resolution chlorophyll-a data to see daily changes.

results = earthaccess.search_data(
    short_name="PACE_OCI_L3M_CHL",
    temporal=tspan,
    granule_name="*.DAY.*.0p1deg.*",
)
paths = earthaccess.open(results)

Instead of xr.open_dataset, we use xr.open_mfdataset to create a single xarray.Dataset (the “mf” in open_mfdataset stands for multiple files) from an array of paths. The paths list is sorted temporally by default, which means the shape of the paths array specifies the way we need to tile the files together into larger arrays. We specify combine="nested" to combine the files according to the shape of the array of files (or file-like objects), even though paths is not a “nested” list in this case. The concat_dim="date" argument generates a new dimension in the combined dataset, because “date” is not an existing dimension in the individual files.

dataset = xr.open_mfdataset(
    paths,
    combine="nested",
    concat_dim="date",
)

Add a date variable as a coordinate using the dates from the netCDF files.

dates = [xr.open_dataset(i).attrs["time_coverage_end"] for i in paths]
dt = pd.to_datetime(dates)
dataset = dataset.assign_coords(date=dt.values)
dataset

<xarray.Dataset> Size: 181MB
Dimensions:  (date: 7, lat: 1800, lon: 3600, rgb: 3, eightbitcolor: 256)
Coordinates:
  * date     (date) datetime64[ns] 56B 2024-05-02T02:28:09.020000 ... 2024-05...
  * lat      (lat) float32 7kB 89.95 89.85 89.75 89.65 ... -89.75 -89.85 -89.95
  * lon      (lon) float32 14kB -179.9 -179.9 -179.8 ... 179.8 179.9 180.0
Dimensions without coordinates: rgb, eightbitcolor
Data variables:
    chlor_a  (date, lat, lon) float32 181MB dask.array<chunksize=(1, 512, 1024), meta=np.ndarray>
    palette  (date, rgb, eightbitcolor) uint8 5kB dask.array<chunksize=(1, 3, 256), meta=np.ndarray>
Attributes: (12/64)
    product_name:                      PACE_OCI.20240501.L3m.DAY.CHL.V3_1.chl...
    instrument:                        OCI
    title:                             OCI Level-3 Standard Mapped Image
    project:                           Ocean Biology Processing Group (NASA/G...
    platform:                          PACE
    source:                            satellite observations from OCI-PACE
    ...                                ...
    identifier_product_doi:            10.5067/PACE/OCI/L3M/CHL/3.1
    keywords:                          Earth Science > Oceans > Ocean Chemist...
    keywords_vocabulary:               NASA Global Change Master Directory (G...
    data_bins:                         736046
    data_minimum:                      0.005030712
    data_maximum:                      98.41369

A common reason to generate a single dataset from multiple, daily images is to create a composite. Compare the map from a single day …

chla = np.log10(dataset["chlor_a"])
chla.attrs.update(
    {
        "units": f'log({dataset["chlor_a"].attrs["units"]})',
    }
)
im = chla.sel(date="2024-05-02").plot(aspect=2, size=4, cmap="GnBu_r")

../../_images/faa899fd05ec7fef677870b27cc50b351b1e74f7635ce0cd312754648b8ee0f2.png

… to a map of average values, which ignore missing values due to clouds.

chla_avg = chla.mean("date", keep_attrs=True)
im = chla_avg.plot(aspect=2, size=4, cmap="GnBu_r")

../../_images/d0fb0b60a0bc38c0e863b8fe90b77a02ae3d44d3f5cb817190ce41a9a7839358.png

We can also create a time series of mean values over the whole region.

chla_avg = chla.mean(dim=["lon", "lat"], keep_attrs=True)
im = chla_avg.plot(linestyle="-", marker="o", color="b")

../../_images/7261bf454e96a543135fe68d21ecde224a764d5692f06307f672dfbed33416a5.png

You have completed the notebook on OCI file structure.

File Structure at Three Processing Levels for the Ocean Color Instrument (OCI)

Contents

File Structure at Three Processing Levels for the Ocean Color Instrument (OCI)#

Summary#

Learning Objectives#

1. Setup#

2. Explore L1B File Structure#

3. Explore L2 File Structure#

4. Explore L3M File Structure#