Visualize HARP2 L2 V4.0 aerosol over ocean product (FastMAPOL)

Visualize HARP2 L2 V4.0 aerosol over ocean product (FastMAPOL)#

Authors: Meng Gao (NASA/SSAI), Sean Foley (NASA/MSU), Kamal Aryal (NASA/UMBC)

Summary#

This notebook explores the HARP2 Level 2 (L2) aerosol product derived from the joint aerosol and surface retrieval algorithm: FastMAPOL (Fast Multi-Angle Polarimetric Ocean and Land algorithm). For more detailed information about the algorithm, please refer to the relevant documentation.

Similar to the SPEXone notebook, we will analyze a scene from the Los Angeles wildfire, which includes both smoke and dust events. The analysis will focus on aerosol optical depth, absorption, and size information over ocean surface (MAPOL_OCEAN).

This notebook updates the previous V3.0 example to use the new V4.0 data. Note that the HARP2 Level-1 (L1) data is still undergoing calibration improvements, which may affect the quality of the Level-2 (L2) products. Several quality metrics are provided to help assess the reliability of the L2 data, and these are reviewed at the end of this tutorial.

How to Cite#

If you use the PACE HARP2 MAPOL_OCEAN Version 4.0 data, please refer to the dataset information page for the latest citation and DOI:

https://www.earthdata.nasa.gov/data/catalog/ob-cloud-pace-harp2-l2-mapol-ocean-4.0

The dataset may be cited as:

NASA Ocean Biology Processing Group. (2026). PACE HARP2 Level-2 Regional Aerosol Over Ocean Optical Properties, FastMAPOL Algorithm, Version 4.0. NASA Ocean Biology Distributed Active Archive Center (OB.DAAC). DOI: 10.5067/PACE/HARP2/L2/MAPOL_OCEAN/4.0. Accessed on: YYYY-MM-DD.

Learning Objectives#

By the end of this notebook, you will understand:

How to acquire HARP2 L2 data
What aerosol products are available
How to visualize basic aerosol properties
How to evaluate data quality

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 requests
import earthaccess
import numpy as np
import xarray as xr
from pathlib import Path

import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from matplotlib.patches import Ellipse
import matplotlib.patheffects as path_effects
import cartopy.crs as ccrs
import cartopy.feature as cfeature

auth = earthaccess.login(persist=True)
fs = earthaccess.get_fsspec_https_session()

4. Visulize HARP2 L2 aerosol properties#

In this example, we visualize the aerosol properties for a scene during LA wild fire with both smoke and dust events. We read the total aerosol optical depth, single scattering albedo, and fine mode volume fraction as below:

aot = dataset["aot"].values
ssa = dataset["ssa"].values
fvf = dataset["fvf"].values
aot.shape, ssa.shape, fvf.shape

((395, 519, 4), (395, 519, 4), (395, 519))

We also need the spatial and angle dimensions as below:

lat = dataset["latitude"].values
lon = dataset["longitude"].values
plot_range = [lon.min(), lon.max(), lat.min(), lat.max()]
wavelength = dataset["wavelength"].values
print(wavelength)

[440.159 549.465 664.564 865.283]

def plot_l2_product(
    lon, lat, data,
    plot_range, label, title,
    vmin=None, vmax=None,
    figsize=(12, 4),
    ellipses=[],
    cmap="viridis",
    log_scale=False,
    land_color="#f2efe9",
    ocean_color="#dbe9f6",
):
    """Make map + histogram with optional log color scaling and land/ocean background.

    Notes:
      - Assumes lon/lat are 2D or 1D arrays available in the outer scope
        (or change signature to pass them in).
      - For log_scale=True, only positive values are used for autoscaling and histogram.
    """

    # ------------------
    # Determine vmin / vmax if not given
    # ------------------
    valid = data[np.isfinite(data)]
    if valid.size == 0:
        raise ValueError("No finite values in `data`.")

    if log_scale:
        valid = valid[valid > 0]
        if valid.size == 0:
            raise ValueError("log_scale=True but `data` has no positive finite values.")

        if vmin is None:
            vmin = np.percentile(valid, 2)
        if vmax is None:
            vmax = np.percentile(valid, 98)

        # Safety: avoid invalid/degenerate bounds
        if (vmin is None) or (vmax is None) or (vmin <= 0) or (vmin >= vmax):
            vmin = float(np.min(valid))
            vmax = float(np.max(valid))

    else:
        if vmin is None:
            vmin = np.percentile(valid, 2)
        if vmax is None:
            vmax = np.percentile(valid, 98)

        if (vmin is None) or (vmax is None) or (vmin >= vmax):
            vmin = float(np.min(valid))
            vmax = float(np.max(valid))

    # ------------------
    # Figure layout
    # ------------------
    fig = plt.figure(figsize=figsize)
    gs = fig.add_gridspec(1, 2, width_ratios=[3, 1], wspace=0.3)

    # ------------------
    # Map subplot
    # ------------------
    ax_map = fig.add_subplot(gs[0], projection=ccrs.PlateCarree())
    ax_map.set_extent(plot_range, crs=ccrs.PlateCarree())

    # Land / ocean background (behind data)
    ax_map.add_feature(cfeature.OCEAN, facecolor=ocean_color, zorder=0)
    ax_map.add_feature(cfeature.LAND, facecolor=land_color, zorder=1)

    ax_map.coastlines(resolution="110m", color="black", linewidth=0.8)
    ax_map.gridlines(draw_labels=True)

    norm = LogNorm(vmin=vmin, vmax=vmax) if log_scale else None

    pm = ax_map.pcolormesh(
        lon, lat, data,
        norm=norm,
        vmin=None if log_scale else vmin,
        vmax=None if log_scale else vmax,
        transform=ccrs.PlateCarree(),
        cmap=cmap,
        zorder=2
    )

    cbar = plt.colorbar(pm, ax=ax_map, orientation="vertical", pad=0.1)
    cbar.set_label(label)

    ax_map.set_title(title, fontsize=12)

    # ------------------
    # Histogram subplot
    # ------------------
    ax_hist = fig.add_subplot(gs[1])

    hist_data = data[np.isfinite(data)]
    if log_scale:
        hist_data = hist_data[hist_data > 0]
        bins = np.logspace(np.log10(vmin), np.log10(vmax), 40)
        ax_hist.hist(hist_data, bins=bins, color="gray", edgecolor="black")
        ax_hist.set_xscale("log")
    else:
        ax_hist.hist(
            hist_data, bins=40, range=[vmin, vmax],
            color="gray", edgecolor="black"
        )

    ax_hist.set_xlabel(label)
    ax_hist.set_ylabel("Count")
    ax_hist.set_title(f"Histogram: N={hist_data.size}")
    for ellipse in ellipses:
        ell = ax_map.add_patch(Ellipse(facecolor=(0, 0, 0, 0), lw=2, zorder=3, linestyle="dashed", **ellipse))
        ell.set_path_effects([path_effects.Stroke(linewidth=4, foreground='white'), path_effects.Normal()])
    ax_map.legend(loc="upper right")

    plt.show()

wavelength_index = 1
title = "Aerosol Optical Depth (AOD): " + str(wavelength[wavelength_index]) + " nm"
label = "AOD"
data = aot[:, :, wavelength_index]
ellipses = [
    {"label": "Smoke", "xy": (-118.5, 31), "width": 3, "height": 8, "angle": 0, "edgecolor": (0.004, 0.451, 0.698)},
    {"label": "Dust", "xy": (-113.5, 24), "width": 10, "height": 11, "angle": 15, "edgecolor": (0.871, 0.561, 0.02)},
]
plot_l2_product(
    lon, lat, data, plot_range=plot_range, label=label, title=title, vmin=0, vmax=0.5, ellipses=ellipses, cmap="jet"
)

../../_images/99d016d0b4a926fb7c4d8f5dfe3c3ddaa3942888a0170a9e869cac1098f78965.png

wavelength_index = 1
title = "Single scattering albedo (SSA): " + str(wavelength[wavelength_index]) + " nm"
label = "SSA"
data = ssa[:, :, wavelength_index]
plot_l2_product(
    lon, lat, data, plot_range=plot_range, label=label, title=title, vmin=0.7, vmax=1, ellipses=ellipses, cmap="jet"
)

../../_images/b4966cb7159d68ad0f725eae1bfa9c348e5e1d65b39a5705c215b211908756b2.png

wavelength_index = 1
title = "Fine mode fraction"
label = "FVF"
data = fvf
plot_l2_product(
    lon, lat, data, plot_range=plot_range, label=label, title=title, vmin=0, vmax=1, ellipses=ellipses, cmap="jet"
)

../../_images/5df8c68f15476084012b3ba0c442507e338d7b1d2e338d72c376ab04fa610c07.png

5. Improve data quality: filter low AOD pixels#

Aerosol absorption and microphysics have larger uncertainties when aerosol loading is low. User can further remove low AOD cases when necessary. We can clearly see, two type of aerosol events with relatively high AOD, the upper one with high absorption (low SSA) and small size (high FVF), probably smoke due to fire; the lower one with less absorption (high SSA) and large size (low FVF), probably dust.

wavelength_index = 1
aot_min = 0.1
title = (
    "Filtered single scattering albedo (SSA): "
    + str(wavelength[wavelength_index])
    + " nm (AOD 550>"
    + str(aot_min)
    + ")"
)
label = "SSA"
data = filtered_ssa = np.where(
    aot[:, :, wavelength_index] >= aot_min, ssa[:, :, wavelength_index], np.nan
)
plot_l2_product(
    lon, lat,data, plot_range=plot_range, label=label, title=title, vmin=0.7, vmax=1, ellipses=ellipses, cmap="jet"
)

../../_images/82f3218cab4448273ba68931ca59a265e23386a8695e5ef33431d6b724c01320.png

The difference in appearance (after matplotlib automatically normalizes the data) is negligible, but the difference in the physical meaning of the array values is quite important.

wavelength_index = 1
aot_min = 0.1
title = "Fine mode fraction (AOD 550>" + str(aot_min) + ")"
label = "FVF"
data = filtered_ssa = np.where(aot[:, :, wavelength_index] >= aot_min, fvf, np.nan)
plot_l2_product(
    lon, lat,data, plot_range=plot_range, label=label, title=title, vmin=0, vmax=1, ellipses=ellipses, cmap="jet"
)

../../_images/d7d92f497f1a8df458aa41cdde38bbb7df45849e4214ea7f20489266a1bb50ca.png

6. Advanced quality assessment#

Since the retrieval algorithm is based on optimal estimation by minimizing a \(\chi^2\) cost function defined as the difference between measurement (m) and forward model fitting (f), normalized by total uncertainties (\(\sigma\)).

\(\chi^2 = \frac{1}{N} \sum (f - m)^2/\sigma^2\)

Here N is the total number of measureents used in retreival. The algorithm also adaptively evalue fitting performance, if the fitting perform poor, it will be removed from the retreival process. Therefore, the \(\chi^2\) and \(N\) can be used to evaluate retrieval performance, the pixels with small \(\chi^2\) (good fitting) and large \(N\) (more pixels can be fitted) will better quality. A more quantitatively approach based on error propogation can be also used to compute retrieval uncertainty, which will be include in future product.

To support L3 data processing, a quality flag is also defined, which is usually based on \(\chi^2\) and \(N\). For the HARP2 test data, we choose

quality_flag = 0: when \(\chi^2<1.5\) and \(N_{ref}>70\) and \(N_{DoLP}>70\)
quality_flag = 1: when \(\chi^2<1.5\) and \(N_{ref}>60\) and \(N_{DoLP}>60\)
quality_flag = 2: when \(\chi^2<1.5\) and \(N_{ref}>40\) and \(N_{DoLP}>40\)
quality_flag = 3: when \(\chi^2<2\) and \(N_{ref}>30\) and \(N_{DoLP}>30\)
quality_flag > 4: for higher value \(\chi^2\) and lower values of \(N_{ref}\) and \(N_{DoLP}\)

The quality flag definition has been revised compared with V3_0 by adding a new quality_flag = 0 criterion. As a result, all existing quality flag values have been renumbered (e.g., 0→1, 1→2, 2→3). The quality flag criteria may be updated further as Level-1 calibration improves.

chi2 = dataset["chi2"].values
nv_ref = dataset["nv_ref"].values
nv_dolp = dataset["nv_dolp"].values
quality_flag = dataset["quality_flag"].values

title = r"Retrieval cost function: $\chi^2$"
label = r"$\chi^2$"
data = chi2
plot_l2_product(
    lon, lat,data, plot_range=plot_range, label=label, title=title, vmin=0, vmax=3, cmap="jet"
)

/tmp/ipykernel_1258/2197280911.py:111: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
  ax_map.legend(loc="upper right")

../../_images/177bcac58df05df697eb89c88f1c844df2fa82db897a0f54b21e7988c8e98006.png

title = r"Total number of reflectance measurements"
label = r"$N_{ref}$"
data = nv_ref
plot_l2_product(
    lon, lat, data, plot_range=plot_range, label=label, title=title, vmin=0, vmax=90, cmap="jet"
)

/tmp/ipykernel_1258/2197280911.py:111: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
  ax_map.legend(loc="upper right")

../../_images/39bd62115496c9f4474246c10000f7b6e6d0d4ee3f7c9046b354d622d81b6f7d.png

title = r"Total number of reflectance measurements"
label = r"$N_{dolp}$"
data = nv_dolp
plot_l2_product(
    lon, lat, data, plot_range=plot_range, label=label, title=title, vmin=0, vmax=90, cmap="jet"
)

/tmp/ipykernel_1258/2197280911.py:111: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
  ax_map.legend(loc="upper right")

../../_images/b7fcf123a56e6ebded361bd686bdecc24a41e32e12eb4208a7ea2f7f8bcee371.png

np.nanmean(chi2), np.nanmean(nv_ref), np.nanmean(nv_dolp)

(np.float32(1.6170229),
 np.float64(57.435484896661364),
 np.float64(41.33418124006359))

Note that \(\chi^2\) converges reasonably well with slight under fitting (averaged around 1.6, peaked around 1.3). Since HARP2 measures 90 angles across 4 bands, the average number of measurement satisfied good fitting are only 52 for reflectance and 42 for polarization, which indicate potential discrepany between forward model and measurements, due to forward model assumptions or likely measurement calibrations.

title = "Retrieval quality flag"
label = "quality_flag"
data = quality_flag
plot_l2_product(
    lon, lat, data, plot_range=plot_range, label=label, title=title, vmin=0, vmax=4, cmap="viridis" 
)

/tmp/ipykernel_1258/2197280911.py:111: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
  ax_map.legend(loc="upper right")

../../_images/fa552d8e71c84d38b73fbec2bc7d9fa6deaabe42b5af2a4fdb57e81782715f24.png

We can evaluate quality flag based on the \(\chi^2\) and \(N\), and only a small portion of data near center of swath reach best quality as we defined as quality_flag=0. With the future improvement of data calibration, more data with better quality will be vailable.

7. Multi-angle data mask for cloud and data screening#

As mentioned previously, FastMAPOL algorithm conducted internal adaptive data screening on each HARP2 angle, the data mask are provided for both reflectance and DoLP. value 0 means the measurements are used in the retrievals, value 1 or NAN means the measurements are removed from retrieval. Therefore, the adaptive data mask can be also used to evaluate fitting quality and measurement quality at each angle. In the example below, please note the difference pattern for reflectance and polarization, which may indicates different calibration perforance.

mask_ref = dataset["mask_ref"].values
mask_dolp = dataset["mask_dolp"].values
mask_ref.shape, mask_dolp.shape

((395, 519, 90, 1), (395, 519, 90, 1))

angle_index = 5
title = "Adaptive data mask on reflectance: angle index " + str(angle_index)
label = "mask_ref"
data = mask_ref[:, :, angle_index, 0]
plot_l2_product(
    lon, lat, data, plot_range=plot_range, label=label, title=title, vmin=0, vmax=1, cmap="viridis"
)

/tmp/ipykernel_1258/2197280911.py:111: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
  ax_map.legend(loc="upper right")

../../_images/05e8e02cf664c225b008669877825e677d4ff6b739229d061e943b1b7d74c7c3.png

angle_index = 5
title = "Adaptive data mask on DoLP: angle index" + str(angle_index)
label = "mask_DOLP"
data = mask_dolp[:, :, angle_index, 0]
plot_l2_product(
    lon, lat, data, plot_range=plot_range, label=label, title=title, vmin=0, vmax=1, cmap="viridis"
)

/tmp/ipykernel_1258/2197280911.py:111: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
  ax_map.legend(loc="upper right")

../../_images/4f51f35e94575586263c08771b849e4704fe0a970d2d2cd43c61538721ea51ee.png

8. Optional: pixel level uncertainty estimation.#

As mentioned previously, pixel level uncertainty can be evalated through error propagation, which propgation measurement uncertainty through Jacobian of the forward model. The estimated uncertainties are discussed in (Gao et al 2021) for HARP2 and AirHARP, but currently not included in the HARP2 L2 products.

9. Reference#

[Validation] Gao, M., Aryal, K., Zhai, P.-W., Knobelspiesse, K., Franz, B. A., Cairns, B., Cetinić, I., Fu, G., Hasekamp, O., Ibrahim, A., Sayer, A. M., and Werdell, P. J.: Where sky meets sea: Integrated aerosol and ocean color retrieval from PACE SPEXone multi-angle polarimetry, Remote Sens. Environ., 345, 115603, https://doi.org/10.1016/j.rse.2026.115603, 2026.
[Algorithm] Gao, M., Franz, B. A., Zhai, P.-W., Knobelspiesse, K., Sayer, A. M., Xu, X., Martins, J. V., Cairns, B., Castellanos, P., Fu, G., Hannadige, N., Hasekamp, O., Hu, Y., Ibrahim, A., Patt, F., Puthukkudy, A., and Werdell, P. J.: Simultaneous retrieval of aerosol and ocean properties from PACE HARP2 with uncertainty assessment using cascading neural network radiative transfer models, Atmos. Meas. Tech., 16, 5863–5881, https://doi.org/10.5194/amt-16-5863-2023, 2023.