--- 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 SPEXone L2 V4.0 aerosol 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 There are two Level-2 (L2) processing suites available for PACE SPEXone data. RemoTAP (Remote Sensing of Trace Gases and Aerosol Products) retrieves aerosol properties over both land and ocean, while FastMAPOL (Fast Multi-Angle Polarimetric Ocean and Land) retrieves both aerosol and ocean color properties over the ocean. This notebook introduces the SPEXone Level-2 (L2) aerosol products generated by the **FastMAPOL** algorithm. For a detailed description of the retrieval algorithm and data products, please refer to the associated documentation. Similar to the HARP2 notebook, we analyze a scene from the Los Angeles wildfire, which includes both smoke and dust events (as observed by OCI and HARP2). However, due to the narrow swath of SPEXone data, the dust event will be the main focus of this tutorial. We will evaluate aerosol optical depth, aerosol absorption, and particle size information. ## 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-spexone-l2-mapol-ocean-4.0 The dataset may be cited as: > NASA Ocean Biology Processing Group. (2026). *PACE SPEXone 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/SPEXone/L2/MAPOL_OCEAN/4.0. Accessed on: YYYY-MM-DD. > ## Learning Objectives By the end of this notebook, you will understand: - How to acquire SPEXone 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 import matplotlib.gridspec as gridspec import cartopy.crs as ccrs import cartopy.feature as cfeature ``` ```{code-cell} ipython3 auth = earthaccess.login(persist=True) ``` ## 2. Get Level-2 Data SPEXone 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 SPEXone 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_SPEXONE_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 SPEXone L2 product structure The SPEXone FastMAPOL L2 product suite includes a long list of aerosol optical properties for both fine and coarse modes (defined in the same format as HARP2 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) ```{code-cell} ipython3 datatree["geophysical_data"] ``` ## 4. Visulize SPEXone 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, label, title, plot_range=None, vmin=None, vmax=None, figsize=(12, 4), 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 # ------------------ mask_valid = np.isfinite(data) valid = data[mask_valid] lat_valid, lon_valid = lat[mask_valid], lon[mask_valid] if plot_range is None: plot_range = [ lon_valid.min(), lon_valid.max(), lat_valid.min(), lat_valid.max() ] 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 = gridspec.GridSpec(16, 48, figure=fig) ax_map = fig.add_subplot(gs[:, :22], projection=ccrs.PlateCarree()) ax_cbar = fig.add_subplot(gs[:, 14:15]) ax_hist = fig.add_subplot(gs[:, 22:]) # ------------------ # Map subplot # ------------------ 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) gl = ax_map.gridlines(draw_labels={"bottom": "x", "left": "y"}) 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, cax=ax_cbar, orientation="vertical", pad=0.2) cbar.set_label(label) ax_map.set_title(title, fontsize=12) # ------------------ # Histogram subplot # ------------------ 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}") plt.show() ``` ```{code-cell} ipython3 wavelength_index = 7 title = "Aerosol Optical Depth (AOD): " + str(wavelength[wavelength_index]) + " nm" label = "AOD" data = aot[:, :, wavelength_index] plot_l2_product( lon, lat, data, label=label, title=title, vmin=0, vmax=0.3, cmap="jet" ) ``` ```{code-cell} ipython3 wavelength_index = 7 title = "Single scattering albedo (SSA): " + str(wavelength[wavelength_index]) + " nm" label = "SSA" data = filtered_ssa = np.where( aot[:, :, wavelength_index] > 0.1, ssa[:, :, wavelength_index], np.nan ) data = ssa[:, :, wavelength_index] plot_l2_product( lon, lat, data, label=label, title=title, vmin=0.7, vmax=1, cmap="jet" ) ``` ```{code-cell} ipython3 wavelength_index = 7 title = "Fine mode fraction" label = "FVF" data = fvf plot_l2_product( lon, lat, data, label=label, title=title, vmin=0, vmax=1, cmap="jet" ) ``` We can clearly see the aerosol event with less absorption (high SSA) and large size (low FVF), probably dust. +++ ## 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. ```{code-cell} ipython3 wavelength_index = 7 aot_min = 0.05 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, label=label, title=title, vmin=0.7, vmax=1, 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 = 7 aot_min = 0.05 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, label=label, title=title, vmin=0, vmax=1, 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 SPEXone L2 data, we choose - quality_flag = 0: when $\chi^2<1.5$ and $N_{ref}>160$ and $N_{DoLP}>160$ - quality_flag = 1: when $\chi^2<1.5$ and $N_{ref}>140$ and $N_{DoLP}>140$ - quality_flag = 2: when $\chi^2<2$ and $N_{ref}>140$ and $N_{DoLP}>140$ - quality_flag = 3: when $\chi^2<2$ and $N_{ref}>100$ and $N_{DoLP}>100$ - quality_flag > 4: for higher value $\chi^2$ and lower values of $N_{ref}$ and $N_{DoLP}$ The quality flag definition remained the same with V3_0 data. ```{code-cell} ipython3 chi2 = dataset["chi2"].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, label=label, title=title, vmin=0, vmax=3, cmap="jet" ) ``` Note that $\chi^2$ converges reasonably well with peak at 1. There are also pixels which do not converge well with relatively high cost function. These pixels need to be removed for more detailed analysis. ```{code-cell} ipython3 title = "Retrieval quality flag" label = "quality_flag" data = quality_flag plot_l2_product( lon, lat, data, label=label, title=title, vmin=0, vmax=3, cmap="viridis" ) ``` We can evaluate quality flag based on the $\chi^2$ and $N$, and most pixels have reach to the best quality with quality_flag=0. +++ ## 7. Multi-angle data mask for cloud and data screening +++ As mentioned previously, FastMAPOL algorithm conducted internal adaptive data screening on each SPEXone 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 wavelength_index = 7 angle_index = 0 title = "Adaptive data mask on reflectance: angle index " + str(angle_index) label = "mask_ref" data = mask_ref[:, :, angle_index, wavelength_index] plot_l2_product( lon, lat, data, label=label, title=title, vmin=0, vmax=1, cmap="viridis" ) ``` ```{code-cell} ipython3 wavelength_index = 7 angle_index = 0 title = "Adaptive data mask on DoLP: angle index " + str(angle_index) label = "mask_dolp" data = mask_dolp[:, :, angle_index, wavelength_index] plot_l2_product( lon, lat, data, label=label, title=title, vmin=0, vmax=1, cmap="viridis" ) ``` ## 8. [TBD]: 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.