Visualize PACE-PAX Data

Visualize PACE-PAX Data#

Authors: Anna Windle (NASA, SSAI)
Adapted from code developed by: Ivona Cetinić (NASA, MSU), Kirk Knobelspiesse (NASA)

Last updated: October 28, 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#

The PACE Postlaunch Airborne eXperiment (PACE-PAX) was a field campaign conducted in September 2024 in California and adjacent coastal areas. PACE-PAX brought together airborne and ship-based in situ observations to collect complementary atmospheric and oceanic measurements for the validation of observations from NASA’s PACE mission.

All PACE-PAX datasets are publicly available through the NASA Earthdata Cloud. This tutorial demonstrates how to access multiple PACE-PAX datasets, combine observations from different platforms, and create integrated visualizations to explore coupled atmosphere–ocean processes. This tutorial generates a multi-platform figure corresponding to Figure 9 in Knobelspiesse et al. (in review).

We will use the seabass-utils Python library to read in some of the PACE-PAX data. The SeaWiFS Bio-optical Archive and Storage System (SeaBASS) is a public archive of in situ oceanographic and atmospheric data maintained by the Ocean Biology Processing Group (OBPG) at NASA Goddard Space Flight Center (GSFC). seabass-utils is a collection of Python functions, also developed by the OBPG, to assist with SeaBASS data processing, file manipulation, and mapping tasks.

Learning Objectives#

At the end of this notebook you will be able to:

Access PACE-PAX datasets from the NASA Earthdata Cloud
Read and explore PACE-PAX datasets from SeaBASS
Combine observations from multiple platforms into a single data fusion visualization
Explore coupled atmosphere–ocean observations collected during PACE-PAX

1. Setup#

Begin by importing all of the packages used in this notebook. If you followed the guidance on the Getting-Started page, then the imports will be successful.

import cartopy.crs as ccrs
import earthaccess
import matplotlib.pyplot as plt
import pandas as pd
import plotly.graph_objects as go
import plotly.io as pio
import seabass_utils as sb
import xarray as xr

Set your Earthdata Login credentials. You can add persist=True to save your credentials in plain text in a netrc file, allowing you to authenticate automatically in the future sessions without re-entering your Earthdata credentials.

auth = earthaccess.login()

2. Access PACE-OCI data#

Let’s find a PACE-OCI Level-2 Inherent Optical Property (IOP) granule acquired during the PACE-PAX campaign over the Santa Barbara Channel, California. The IOP suite contains particulate backscattering data.

tspan = ("2024-09-26", "2024-09-26")
bbox = (-120, 33, -110, 35)
results = earthaccess.search_data(
    short_name="PACE_OCI_L2_IOP",
    temporal=tspan,
    bounding_box=bbox,
)
len(results)

oci_paths = earthaccess.open(results)

We can open this file up using xarray and retrieve particulate backscattering at 442 nm data (bbp_442):

oci_dt = xr.open_datatree(oci_paths[-1])
oci_bbp = oci_dt["geophysical_data"]["bbp_442"]
for item in ("longitude", "latitude"):
    oci_bbp[item] = oci_dt["navigation_data"][item]
oci_bbp

<xarray.DataArray 'bbp_442' (number_of_lines: 1710, pixels_per_line: 1272)> Size: 9MB
[2175120 values with dtype=float32]
Coordinates:
    longitude  (number_of_lines, pixels_per_line) float32 9MB ...
    latitude   (number_of_lines, pixels_per_line) float32 9MB ...
Dimensions without coordinates: number_of_lines, pixels_per_line
Attributes:
    long_name:   Particulate backscattering at 442 nm, GIOP model
    units:       m^-1
    valid_min:   -27999
    valid_max:   2000
    reference:   P.J. Werdell, B.A. Franz, S.W. Bailey, G.C. Feldman and 15 c...
    wavelength:  442.327

Let’s subset the PACE-OCI bbp_442 data to a smaller geographical region in the Santa Barbara Channel, retain only positive values, and plot it:

subset = (
    (oci_bbp["longitude"] >= -120.5)
    & (oci_bbp["longitude"] <= -119.0)
    & (oci_bbp["latitude"] >= 33.8)
    & (oci_bbp["latitude"] <= 34.5)
)
oci_bbp = oci_bbp.where(subset & (oci_bbp >= 0), drop=True)
oci_bbp.sizes

Frozen({'number_of_lines': 79, 'pixels_per_line': 108})

crs = ccrs.PlateCarree()
fig, ax = plt.subplots(figsize=(8, 6), subplot_kw={"projection": crs})

oci_bbp.plot(
    x="longitude",
    y="latitude",
    cmap="viridis",
    vmax=0.015,
    cbar_kwargs={"shrink": 0.65},
    ax=ax,
)

ax.coastlines(resolution="10m", linewidth=1)
ax.gridlines(
    draw_labels={"bottom": True, "left": True, "top": False, "right": False},
    xlocs=plt.MaxNLocator(5),
    ylocs=plt.MaxNLocator(5),
    linewidth=0,
)

plt.tight_layout()
plt.show()

../_images/d19281aea006622360a7876fc3fa1d9c7611ba59a7b943eea7e5fe55d822a9e9.png

Higher bbp_442 values are observed along the coast, indicating greater particulate backscattering and potentially higher particle concentrations.

5. Access ship cast data#

In situ oceanic observations were collected aboard the NOAA R/V Shearwater. Particulate backscattering was measured using an SC6 backscattering sensor.

Let’s open the SC6 data collected on Sep 26th:

results = earthaccess.search_data(
    short_name="PACE-PAX",
    temporal=tspan,
    granule_name="*Shearwater*SC6026*",
)
len(results)

These granules each reference multiple files containing data for that granule. At least one file for each granule is a plain .sb, rather than a .tgz.sb, indicating that it is a SeaBASS file.

urls = []
for item in results:
    urls += [i for i in item.data_links() if not i.endswith(".tgz.sb")]
urls

['https://oceandata.sci.gsfc.nasa.gov/cmr/getfile/80370ff860_PACE-PAX_RV_Shearwater_SC6026_st38_R1.sb',
 'https://oceandata.sci.gsfc.nasa.gov/cmr/getfile/80370ff860_PACE-PAX_RV_Shearwater_SC6026_st39_R1.sb',
 'https://oceandata.sci.gsfc.nasa.gov/cmr/getfile/80370ff860_PACE-PAX_RV_Shearwater_SC6026_st40_R1.sb',
 'https://oceandata.sci.gsfc.nasa.gov/cmr/getfile/80370ff860_PACE-PAX_RV_Shearwater_SC6026_st37_R1.sb']

sc6_paths = earthaccess.open(urls)

SeaBASS files can be read using the seabass_utils library, which we imported above as sb. The sb_read() function reads the SeaBASS file into a dictionary that pandas converts into a DataFrame.

Station information is extracted from header metadata and added as pandas.DataFrame columns:

sc6_bbp = []

for file in sc6_paths:
    data = sb.sb_read(file, no_warn=True)
    df = pd.DataFrame(data.data)
    headers = data.headers

    lat = float(headers["north_latitude"].split("[")[0])
    lon = float(headers["east_longitude"].split("[")[0])
    date = headers["start_date"]
    time = headers["start_time"].split("[")[0]

    bbp = df[["depth", "bbp440"]].assign(
        **{
            "date": date,
            "time": time,
            "latitude": lat,
            "longitude": lon,
        }
    )
    sc6_bbp.append(bbp)

sc6_bbp = pd.concat(sc6_bbp, ignore_index=True)
sc6_bbp.head()

	depth	bbp440	date	time	latitude	longitude
0	0.4405	0.006945	20240926	16:41:02	34.1889	-119.6281
1	0.9380	0.007183	20240926	16:41:02	34.1889	-119.6281
2	1.9430	0.007060	20240926	16:41:02	34.1889	-119.6281
3	2.8510	0.007014	20240926	16:41:02	34.1889	-119.6281
4	4.1280	0.007370	20240926	16:41:02	34.1889	-119.6281

And we can now plot the vertical profiles of particulate backscattering at 440 nm (bbp440):

fig, ax = plt.subplots(figsize=(5, 6))

sc = ax.scatter(
    sc6_bbp["longitude"],
    sc6_bbp["depth"],
    c=sc6_bbp["bbp440"],
    s=35,
    cmap="viridis",
)

ax.invert_yaxis()
ax.set_xlabel("Longitude")
ax.set_ylabel("Depth (m)")
ax.set_title("R/V Shearwater SC6 $b_{bp}(440)$ profiles")

cbar = plt.colorbar(sc, ax=ax)
cbar.set_label("$b_{bp}(440)$ (m$^{-1}$)")

plt.tight_layout()
plt.show()

../_images/f18342114de0f7d33dac228b0f540ce41424dd822ac82f3ef97d1a6ef16ee4df.png

6. Access ocean glider data#

An autonomous ocean glider was also deployed along a 10 km transect for 25 days. The glider was equipped with a Seabird ECO Triplet optical sensor that measured particulate backscatter at 532 nm.

results = earthaccess.search_data(
    short_name="PACE-PAX",
    temporal=tspan,
    granule_name="*glider*bbp*",
)
len(results)

eco_paths = earthaccess.open(results)

We will use sb.read() to open this data as well:

eco_data = sb.sb_read(eco_paths[0], no_warn=True)
eco_df = pd.DataFrame(eco_data.data)
eco_df.head()

	date	time	lat	lon	depth	wt	wt_sd	wt_bincount	cond	cond_sd	...	bbp532_bincount	bbp700	bbp700_sd	bbp700_bincount	bb532	bb532_sd	bb532_bincount	bb700	bb700_sd	bb700_bincount
0	20240907	21:23:10	34.26693	-119.86806	0.155	17.3445	4.2924	106	45.0698	0.0107	...	104	0.005235	0.001075	107	0.003669	0.000506	104	0.005540	0.001088	107
1	20240907	21:35:46	34.26700	-119.87155	1.477	19.0290	0.0403	7	45.0329	0.0529	...	7	0.004983	0.000242	7	0.003775	0.000131	7	0.005294	0.000242	7
2	20240907	21:35:52	34.26699	-119.87159	2.503	18.4357	0.3851	6	44.3475	0.3927	...	6	0.004860	0.000210	6	0.004545	0.000637	6	0.005171	0.000210	6
3	20240907	21:35:57	34.26699	-119.87161	3.459	17.7523	0.0431	4	43.7660	0.0982	...	4	0.004834	0.000109	4	0.004766	0.000540	4	0.005145	0.000109	4
4	20240907	21:36:02	34.26700	-119.87162	4.482	17.4256	0.1180	5	43.4607	0.0818	...	5	0.004660	0.000094	5	0.004438	0.000191	5	0.004971	0.000094	5

5 rows × 35 columns

Let’s filter out glider data for a specified date and location and keep only positive bbp532 values:

subset = (
    (eco_df["date"] == 20240926)
    & eco_df["lon"].between(-120.5, -119.0)
    & eco_df["lat"].between(33.8, 34.5)
)

eco_bbp = eco_df.loc[subset & (eco_df["bbp532"] > 0)]
eco_bbp.shape

(12400, 35)

And we can plot bbp(532) along the glider’s track as it descended and ascended throughout the day:

fig, ax = plt.subplots(figsize=(8, 5))

sc = ax.scatter(
    eco_bbp["lon"],
    eco_bbp["depth"],
    c=eco_bbp["bbp532"],
    s=15,
    cmap="viridis",
    vmax=0.004,
)

ax.invert_yaxis()
ax.set_xlabel("Longitude")
ax.set_ylabel("Depth [m]")
ax.set_title("PACE-PAX Glider $b_{bp}(532)$ on Sept 26 2024")

cbar = plt.colorbar(sc, ax=ax)
cbar.set_label("$b_{bp}(532)$ [m$^{-1}$]")

plt.tight_layout()
plt.show()

../_images/80f14af3c7e68433ac818204c4ea0aef7510489b036cdd5a1e5cd1f76101d1b0.png

7. Muli-platform Visualization#

In order to make all these data visually comparable, there are couple units discrepancies to smooth out. Ocean data use depth, which we’ll standardize to a “z” variable consistent with the altitude as used by airborne platforms. Both are in meters, just different conventions for positive.

sc6_bbp["z"] = -sc6_bbp["depth"]
eco_bbp["z"] = -eco_bbp["depth"]

The backscattering coefficients from the Twin Otter and ER-2 platforms are not the same units. We can convert the ER-2 platform’s data to units of Mm\(^{-1}\) sr\(^{-1}\) for plotting together.

hsrl2_bbp = hsrl2_bbp * 1000
hsrl2_bbp.attrs["units"] = "Mm-1.sr-1"

Now, we are able to show all this data together in one 3-dimensional plot:

pio.renderers.default = "notebook"

fig = go.Figure()

ocean_cmin = 0
ocean_cmax = 0.015
atmos_cmin = 0
atmos_cmax = 4

# --- PACE ---
fig.add_trace(
    go.Surface(
        x=oci_bbp["longitude"],
        y=oci_bbp["latitude"],
        z=oci_bbp.where(oci_bbp.isnull(), 0),
        surfacecolor=oci_bbp,
        colorscale="viridis",
        showscale=True,
        opacity=0.8,
        cmin=ocean_cmin,
        cmax=ocean_cmax,
        colorbar={
            "title": {
                "text": "Ocean b<sub>bp</sub>(λ) [m⁻¹]",
                "side": "right",
                "font": {"size": 10},
            },
            "tickfont": {"size": 9},
            "x": 0.8,
            "y": 0.2,
            "len": 0.4,
            "tickvals": [0, 0.002, 0.004, 0.006, 0.008, 0.01],
            "ticktext": ["0", "0.002", "0.004", "0.006", "0.008", "0.01"],
        },
    )
)

# --- Shearwater ---
fig.add_trace(
    go.Scatter3d(
        x=sc6_bbp["longitude"],
        y=sc6_bbp["latitude"],
        z=sc6_bbp["z"],
        mode="markers",
        marker={
            "size": 3,
            "color": sc6_bbp["bbp440"],
            "colorscale": "viridis",
            "showscale": False,
            "cmin": ocean_cmin,
            "cmax": ocean_cmax,
        },
    )
)

# --- Glider ---
fig.add_trace(
    go.Scatter3d(
        x=eco_bbp["lon"],
        y=eco_bbp["lat"],
        z=eco_bbp["z"],
        mode="markers",
        marker={
            "symbol": "square",
            "size": 3,
            "color": eco_bbp["bbp532"],
            "colorscale": "viridis",
            "showscale": False,
            "cmin": ocean_cmin,
            "cmax": ocean_cmax,
        },
        line={"color": "lightgray", "width": 2},
    )
)

# --- ER-2 ---
x, y, z, c = xr.broadcast(hsrl2_bbp["lon"], hsrl2_bbp["lat"], hsrl2_bbp["z"], hsrl2_bbp)
fig.add_trace(
    go.Surface(
        x=x,
        y=y,
        z=z.where(c.notnull()),
        surfacecolor=c,
        colorscale="inferno_r",
        showscale=True,
        cmin=atmos_cmin,
        cmax=atmos_cmax,
        colorbar={
            "title": {
                "text": "Aerosol b<sub>bp</sub>(532) [Mm⁻¹sr⁻¹]",
                "side": "right",
                "font": {"size": 10},
            },
            "tickfont": {"size": 9},
            "x": 0.8,
            "y": 0.6,
            "len": 0.4,
        },
    )
)


# --- Twin Otter ---
fig.add_trace(
    go.Scatter3d(
        x=isara_bbp["longitude"],
        y=isara_bbp["latitude"],
        z=isara_bbp["altitude"],
        mode="markers",
        marker={
            "size": 3,
            "color": isara_bbp,
            "colorscale": "plasma_r",
            "showscale": False,
            "cmin": atmos_cmin,
            "cmax": atmos_cmax,
            "line": {
                "color": "black",
                "width": 1,
            },
        },
    )
)

fig.update_layout(
    title={
        "text": "Multi-Platform 3D Ocean/Atmosphere Visualization (Santa Barbara)",
        "y": 0.995,
        "yanchor": "top",
        "x": 0.4,
        "xanchor": "center",
    },
    scene={
        "aspectmode": "manual",
        "aspectratio": {"x": 1, "y": 1, "z": 0.8},
        "xaxis": {"title": "Longitude", "range": [-120.2, -119.2]},
        "yaxis": {"title": "Latitude", "range": [34.15, 34.45]},
        "zaxis": {"title": "Height [m]", "range": [-400, 1000]},
        "camera": {
            "eye": {"x": -1.5, "y": -1.5, "z": 0.1}
        },
    },
    width=1100,
    height=800,
    margin={"l": 0, "r": 100, "b": 200, "t": 0},
    showlegend=False,
)

annotations_3d = [
    {
        "x": -120.1,
        "y": 34.4,
        "z": 5,
        "text": "<b>PACE</b>",
        "showarrow": False,
        "font": {"size": 12, "color": "white"},
        "bgcolor": "rgba(0,0,0,0.7)",
        "borderpad": 4,
    },
    {
        "x": -119.6,
        "y": 34.22,
        "z": -100,
        "text": "<b>Shearwater</b>",
        "showarrow": False,
        "font": {"size": 12, "color": "blue"},
        "bgcolor": "rgba(255,255,255,0.85)",
        "borderpad": 4,
    },
    {
        "x": -120.05,
        "y": 34.3,
        "z": -300,
        "text": "<b>Glider</b>",
        "showarrow": False,
        "font": {"size": 12, "color": "green"},
        "bgcolor": "rgba(255,255,255,0.85)",
        "borderpad": 4,
    },
    {
        "x": -119.6,
        "y": 34.2,
        "z": 800,
        "text": "<b>Twin Otter</b>",
        "showarrow": False,
        "font": {"size": 12, "color": "purple"},
        "bgcolor": "rgba(255,255,255,0.85)",
        "borderpad": 4,
    },
    {
        "x": -120.07,
        "y": 34.41,
        "z": 900,
        "text": "<b>ER-2</b>",
        "showarrow": False,
        "font": {"size": 12, "color": "orange"},
        "bgcolor": "rgba(255,255,255,0.85)",
        "borderpad": 4,
    },
]

fig.update_layout(scene={"annotations": annotations_3d})
fig.show()

You have completed the notebook on accessing and analyzing PACE-PAX data!