Add support for reading and visualizing PACE data (#20)

* Add read_pace function

* Change to viz_pace

* Add projection

* Add PACE notebook

* Add cartopy

docs/examples/pace.ipynb

@@ -0,0 +1,121 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "[![image](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/opengeos/HyperCoast/blob/main/docs/examples/pace.ipynb)\n",
+    "\n",
+    "# Visualizing PACE data interactively with HyperCoast\n",
+    "\n",
+    "This notebook demonstrates how to visualize [Plankton, Aerosol, Cloud, ocean Ecosystem (PACE)](https://pace.gsfc.nasa.gov) data interactively with HyperCoast."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %pip install \"hypercoast[extra]\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import hypercoast"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Download a sample PACE data file from [here](https://github.com/opengeos/datasets/releases/tag/netcdf)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "url = \"https://github.com/opengeos/datasets/releases/download/netcdf/PACE_OCI.20240423T184658.L2.OC_AOP.V1_0_0.NRT.nc\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "filepath = \"PACE_OCI.20240423T184658.L2.OC_AOP.V1_0_0.NRT.nc\"\n",
+    "hypercoast.download_file(url)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Load the dataset as a `xarray.Dataset` object."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dataset = hypercoast.read_pace(filepath)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hypercoast.viz_pace(dataset, wavelengths=[500, 510, 520, 530], ncols=2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Add projection."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hypercoast.viz_pace(dataset, wavelengths=[500, 510, 520, 530], ncols=2, crs=\"default\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "hyper",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.14"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

hypercoast/common.py

@@ -2,6 +2,7 @@
 """
 
 import os
+from typing import List
 
 
 def github_raw_url(url):
@@ -116,3 +117,24 @@ def download_file(
                     tar_ref.extractall(os.path.dirname(output))
 
     return os.path.abspath(output)
+
+
+def netcdf_groups(filepath: str) -> List[str]:
+    """
+    Get the list of groups in a NetCDF file.
+
+    Args:
+        filepath (str): The path to the NetCDF file.
+
+    Returns:
+        list: A list of group names in the NetCDF file.
+
+    Example:
+        >>> netcdf_groups('path/to/netcdf/file')
+        ['group1', 'group2', 'group3']
+    """
+    import h5netcdf
+
+    with h5netcdf.File(filepath) as file:
+        groups = list(file)
+    return groups

hypercoast/hypercoast.py

@@ -3,7 +3,7 @@
 import ipyleaflet
 import leafmap
 import xarray as xr
-from .common import download_file
+from .common import download_file, netcdf_groups
 from .emit import read_emit, plot_emit, viz_emit, emit_to_netcdf, emit_to_image
 from .pace import *
 

hypercoast/pace.py

@@ -2,6 +2,7 @@
 """
 
 import xarray as xr
+from typing import List, Tuple, Union, Optional
 
 
 def read_pace(filepath, wavelengths=None, method="nearest", **kwargs):
@@ -17,34 +18,161 @@ def read_pace(filepath, wavelengths=None, method="nearest", **kwargs):
     Returns:
         xr.Dataset: An xarray Dataset containing the PACE data.
     """
-    ds = xr.open_dataset(filepath, group="geophysical_data")
-    ds = ds.swap_dims(
+
+    rrs = xr.open_dataset(filepath, group="geophysical_data")["Rrs"]
+    wvl = xr.open_dataset(filepath, group="sensor_band_parameters")
+    dataset = xr.open_dataset(filepath, group="navigation_data")
+    dataset = dataset.set_coords(("longitude", "latitude"))
+    dataset = dataset.rename({"pixel_control_points": "pixels_per_line"})
+    dataset = xr.merge((rrs, dataset.coords))
+    dataset.coords["wavelength_3d"] = wvl.coords["wavelength_3d"]
+    dataset = dataset.rename(
         {
             "number_of_lines": "latitude",
             "pixels_per_line": "longitude",
+            "wavelength_3d": "wavelength",
         }
     )
-    wvl = xr.open_dataset(filepath, group="sensor_band_parameters")
-    loc = xr.open_dataset(filepath, group="navigation_data")
 
-    lat = loc.latitude
-    lat = lat.swap_dims(
-        {"number_of_lines": "latitude", "pixel_control_points": "longitude"}
-    )
+    if wavelengths is not None:
+        dataset = dataset.sel(wavelength=wavelengths, method=method, **kwargs)
 
-    lon = loc.longitude
-    wavelengths = wvl.wavelength_3d
-    Rrs = ds.Rrs
-
-    dataset = xr.Dataset(
-        {"Rrs": (("latitude", "longitude", "wavelengths"), Rrs.data)},
-        coords={
-            "latitude": (("latitude", "longitude"), lat.data),
-            "longitude": (("latitude", "longitude"), lon.data),
-            "wavelengths": ("wavelengths", wavelengths.data),
-        },
-    )
+    return dataset
+
+
+def viz_pace(
+    dataset: Union[xr.Dataset, str],
+    wavelengths: Optional[Union[List[float], float]] = None,
+    method: str = "nearest",
+    figsize: Tuple[float, float] = (6.4, 4.8),
+    cmap: str = "jet",
+    vmin: float = 0,
+    vmax: float = 0.02,
+    ncols: int = 1,
+    crs: Optional[str] = None,
+    xlim: Optional[List[float]] = None,
+    ylim: Optional[List[float]] = None,
+    **kwargs,
+):
+    """
+    Plots PACE data from a given xarray Dataset.
+
+    Args:
+        dataset (xr.Dataset): An xarray Dataset containing the PACE data.
+        wavelengths (array-like, optional): Specific wavelengths to select. If None, all wavelengths are selected.
+        method (str, optional): Method to use for selection when wavelengths is not None. Defaults to "nearest".
+        figsize (tuple, optional): Figure size. Defaults to (6.4, 4.8).
+        cmap (str, optional): Colormap to use. Defaults to "jet".
+        vmin (float, optional): Minimum value for the colormap. Defaults to 0.
+        vmax (float, optional): Maximum value for the colormap. Defaults to 0.02.
+        ncols (int, optional): Number of columns in the plot. Defaults to 1.
+        crs (str or cartopy.crs.CRS, optional): Coordinate reference system to use. If None, a simple plot is created. Defaults to None.
+            See https://scitools.org.uk/cartopy/docs/latest/reference/projections.html
+        xlim (array-like, optional): Limits for the x-axis. Defaults to None.
+        ylim (array-like, optional): Limits for the y-axis. Defaults to None.
+        **kwargs: Additional keyword arguments to pass to the `plt.subplots` function.
+    """
+
+    import matplotlib.pyplot as plt
+    import numpy as np
+    import math
+
+    if isinstance(dataset, str):
+        dataset = read_pace(dataset, wavelengths, method)
 
     if wavelengths is not None:
-        dataset = dataset.sel(wavelengths=wavelengths, method=method, **kwargs)
-    return dataset
+        if not isinstance(wavelengths, list):
+            wavelengths = [wavelengths]
+        dataset = dataset.sel(wavelength=wavelengths, method=method)
+    else:
+        wavelengths = dataset.coords["wavelength"][0].values.tolist()
+
+    lat = dataset.coords["latitude"]
+    lon = dataset.coords["longitude"]
+
+    nrows = math.ceil(len(wavelengths) / ncols)
+
+    if crs is None:
+
+        fig, axes = plt.subplots(
+            nrows=nrows,
+            ncols=ncols,
+            figsize=(figsize[0] * ncols, figsize[1] * nrows),
+            **kwargs,
+        )
+
+        for i in range(nrows):
+            for j in range(ncols):
+                index = i * ncols + j
+                if index < len(wavelengths):
+                    wavelength = wavelengths[index]
+                    data = dataset.sel(wavelength=wavelength, method=method)["Rrs"]
+
+                    if min(nrows, ncols) == 1:
+                        ax = axes[index]
+                    else:
+                        ax = axes[i, j]
+                    im = ax.pcolormesh(
+                        lon, lat, np.squeeze(data), cmap=cmap, vmin=vmin, vmax=vmax
+                    )
+                    ax.set_xlabel("Longitude")
+                    ax.set_ylabel("Latitude")
+                    ax.set_title(
+                        f"wavelength = {dataset.coords['wavelength'].values[index]} [nm]"
+                    )
+                    fig.colorbar(im, ax=ax, label="Reflectance")
+
+        plt.tight_layout()
+        plt.show()
+
+    else:
+
+        import cartopy
+        from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
+
+        if crs == "default":
+            crs = cartopy.crs.PlateCarree()
+
+        if xlim is None:
+            xlim = [math.floor(lon.min()), math.ceil(lon.max())]
+
+        if ylim is None:
+            ylim = [math.floor(lat.min()), math.ceil(lat.max())]
+
+        fig, axes = plt.subplots(
+            nrows=nrows,
+            ncols=ncols,
+            figsize=(figsize[0] * ncols, figsize[1] * nrows),
+            subplot_kw={"projection": cartopy.crs.PlateCarree()},
+            **kwargs,
+        )
+
+        for i in range(nrows):
+            for j in range(ncols):
+                index = i * ncols + j
+                if index < len(wavelengths):
+                    wavelength = wavelengths[index]
+                    data = dataset.sel(wavelength=wavelength, method=method)["Rrs"]
+
+                    if min(nrows, ncols) == 1:
+                        ax = axes[index]
+                    else:
+                        ax = axes[i, j]
+                    im = ax.pcolormesh(lon, lat, data, cmap="jet", vmin=0, vmax=0.02)
+                    ax.coastlines()
+                    ax.add_feature(cartopy.feature.STATES, linewidth=0.5)
+                    ax.set_xticks(np.linspace(xlim[0], xlim[1], 5), crs=crs)
+                    ax.set_yticks(np.linspace(ylim[0], ylim[1], 5), crs=crs)
+                    lon_formatter = LongitudeFormatter(zero_direction_label=True)
+                    lat_formatter = LatitudeFormatter()
+                    ax.xaxis.set_major_formatter(lon_formatter)
+                    ax.yaxis.set_major_formatter(lat_formatter)
+                    ax.set_xlabel("Longitude")
+                    ax.set_ylabel("Latitude")
+                    ax.set_title(
+                        f"wavelength = {dataset.coords['wavelength'].values[index]} [nm]"
+                    )
+                    plt.colorbar(im, label="Reflectance")
+
+        plt.tight_layout()
+        plt.show()

mkdocs.yml

@@ -81,6 +81,7 @@ nav:
     - Report Issues: https://github.com/opengeos/HyperCoast/issues
     - Examples:
           - examples/emit.ipynb
+          - examples/pace.ipynb
     - API Reference:
           - common module: common.md
           - emit module: emit.md

pyproject.toml

@@ -34,7 +34,7 @@ all = [
 ]
 
 extra = [
-    "pandas",
+    "cartopy",
 ]
 
 

requirements_dev.txt

@@ -2,6 +2,8 @@ black
 black[jupyter]
 build
 bump-my-version
+
+cartopy
 click
 codespell
 flake8