--- jupytext: text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.19.5 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- # Visualize HARP2 L2 V4.0 aerosol over ocean product (FastMAPOL) **Authors:** Meng Gao (NASA/SSAI), Sean Foley (NASA/MSU), Kamal Aryal (NASA/UMBC) [edl]: https://urs.earthdata.nasa.gov/ [oci-data-access]: https://oceancolor.gsfc.nasa.gov/resources/docs/tutorials/notebooks/oci_data_access/ ## 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. [tutorials]: https://oceancolor.gsfc.nasa.gov/resources/docs/tutorials/ ```{code-cell} ipython3 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 ``` ```{code-cell} ipython3 auth = earthaccess.login(persist=True) fs = earthaccess.get_fsspec_https_session() ``` ## 2. Get Level-2 Data HARP2 L2 data is available on both OB.DAAC and earth data cloud. Please refer L1C notebook on the access of cloud. The following block retrieves a single HARP2 L2 granule at the time stamp 2025/01/09 20:00:19 UTC. V4.0 data is specified in the search. ```{code-cell} ipython3 results = earthaccess.search_data( short_name="PACE_HARP2_L2_MAPOL_OCEAN", temporal=("2025-01-09T20:00:20", "2025-01-09T20:00:21"), granule_name='*V4_0*', count=1, ) paths = earthaccess.open(results) ``` ```{code-cell} ipython3 # this cell is tagged to be removed from HTML renders, # but we currently want to download when we don't have direct access if not earthaccess.__store__.in_region: paths = earthaccess.download(results, "./") ``` PACE polarimeter L2 products for both HARP2 and SPEXone include four data groups - geolocation_data - geophysical_data - diagnostic_data - sensor_band_parameters ```{code-cell} ipython3 datatree = xr.open_datatree(paths[0]) datatree ``` Here we merge all the data group together for convenience in data manipulations. ```{code-cell} ipython3 dataset = xr.merge(datatree.to_dict().values()) dataset ``` ## 3. Understanding HARP2 L2 product structure The HARP2 FastMAPOL L2 product suite includes a long list of aerosol optical properties for both fine and coarse modes (defined in the same format as SPEXone L2 products): - Aerosol optical depth (aot and aot_fine/coarse) - Aerosol single scattering albedo (ssa and ssa_fine/coarse) - Ångström coefficient (angstrom_440_870 and angstrom_440_670) - Aerosol fine mode optical depth fraction (fmf) - etc As well as aerosol microphysical properties: - Aerosol effective radius (reff_fine/coarse) and variance (veff_fine/coarse) - Aerosol refractive index: real part (mr and mr_fine/coarse), imaginary part (mi and mi_fine/coarse) - Aerosol spherical fraction (sph and sph_fine/coarse) - Aerosol volume density (vd_fine/coarse) - Aerosol fine mode volume fraction (fvf) - Aerosol layer height (alh) - etc And a set of other products: - Wind speed (wind_speed) - Chlorophyll-a (chla) - Remote sensing reflectance (Rrs*) The remote sensing reflectance characterizes ocean-leaving reflectance. It is derived via atmospheric correction based on the retrieved aerosol properties at all HARP2 viewing angles. Therefore, it includes an angle dimension, as in the L1C data. There are two versions of remote sensing reflectance: Rrs_angular (before BRDF correction, renamed from Rrs1 in V3 data) and Rrs_nadir (after BRDF correction, renamed from Rrs2 in V3 data). Rrs_angular and Rrs_nadir are both standard output in L2 file. Meanwhile, their angular means and standard deviations are also included as Rrs_nadir_mean/std and Rrs_angular_mean/std. ```{code-cell} ipython3 :scrolled: true datatree["geophysical_data"] ``` ## 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: ```{code-cell} ipython3 aot = dataset["aot"].values ssa = dataset["ssa"].values fvf = dataset["fvf"].values aot.shape, ssa.shape, fvf.shape ``` We also need the spatial and angle dimensions as below: ```{code-cell} ipython3 lat = dataset["latitude"].values lon = dataset["longitude"].values plot_range = [lon.min(), lon.max(), lat.min(), lat.max()] wavelength = dataset["wavelength"].values print(wavelength) ``` ```{code-cell} ipython3 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() ``` ```{code-cell} ipython3 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" ) ``` ```{code-cell} ipython3 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" ) ``` ```{code-cell} ipython3 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" ) ``` ## 5. Improve data quality: filter low AOD pixels +++ {"lines_to_next_cell": 2} 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. ```{code-cell} ipython3 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" ) ``` 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. ```{code-cell} ipython3 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" ) ``` ## 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. ```{code-cell} ipython3 chi2 = dataset["chi2"].values nv_ref = dataset["nv_ref"].values nv_dolp = dataset["nv_dolp"].values quality_flag = dataset["quality_flag"].values ``` ```{code-cell} ipython3 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" ) ``` ```{code-cell} ipython3 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" ) ``` ```{code-cell} ipython3 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" ) ``` ```{code-cell} ipython3 np.nanmean(chi2), np.nanmean(nv_ref), np.nanmean(nv_dolp) ``` 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. ```{code-cell} ipython3 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" ) ``` 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. ```{code-cell} ipython3 mask_ref = dataset["mask_ref"].values mask_dolp = dataset["mask_dolp"].values mask_ref.shape, mask_dolp.shape ``` ```{code-cell} ipython3 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" ) ``` ```{code-cell} ipython3 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" ) ``` ## 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.