From 81dadfecbafa1100ebc3f3b2c2e01a5cbf694a65 Mon Sep 17 00:00:00 2001 From: Qiusheng Wu Date: Sun, 5 May 2024 14:02:58 -0400 Subject: [PATCH 1/5] Add read_pace function --- hypercoast/common.py | 22 +++++++++ hypercoast/hypercoast.py | 2 +- hypercoast/pace.py | 104 ++++++++++++++++++++++++++++++--------- 3 files changed, 105 insertions(+), 23 deletions(-) diff --git a/hypercoast/common.py b/hypercoast/common.py index e9cb47d..04e4379 100644 --- a/hypercoast/common.py +++ b/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 diff --git a/hypercoast/hypercoast.py b/hypercoast/hypercoast.py index a0b0de9..16ef4ed 100644 --- a/hypercoast/hypercoast.py +++ b/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 * diff --git a/hypercoast/pace.py b/hypercoast/pace.py index 7f17d5f..f0a286e 100644 --- a/hypercoast/pace.py +++ b/hypercoast/pace.py @@ -17,34 +17,94 @@ 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 plot_pace( + dataset, + wavelengths=None, + method="nearest", + cmap="jet", + vmin=0, + vmax=0.02, + ncols=1, + **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". + 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. + """ + + 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) + + fig, axes = plt.subplots( + nrows=nrows, ncols=ncols, figsize=(6.4 * ncols, 4.8 * 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() From 4466414acfc812f09a4179bd3c39ffec600375d1 Mon Sep 17 00:00:00 2001 From: Qiusheng Wu Date: Sun, 5 May 2024 14:03:38 -0400 Subject: [PATCH 2/5] Change to viz_pace --- hypercoast/pace.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hypercoast/pace.py b/hypercoast/pace.py index f0a286e..a99e44e 100644 --- a/hypercoast/pace.py +++ b/hypercoast/pace.py @@ -39,7 +39,7 @@ def read_pace(filepath, wavelengths=None, method="nearest", **kwargs): return dataset -def plot_pace( +def viz_pace( dataset, wavelengths=None, method="nearest", From 91d717ba9ef3ba231cf2505e7a685abcb057780e Mon Sep 17 00:00:00 2001 From: Qiusheng Wu Date: Sun, 5 May 2024 16:41:45 -0400 Subject: [PATCH 3/5] Add projection --- hypercoast/pace.py | 134 ++++++++++++++++++++++++++++++++++----------- 1 file changed, 101 insertions(+), 33 deletions(-) diff --git a/hypercoast/pace.py b/hypercoast/pace.py index a99e44e..0072444 100644 --- a/hypercoast/pace.py +++ b/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): @@ -40,13 +41,17 @@ def read_pace(filepath, wavelengths=None, method="nearest", **kwargs): def viz_pace( - dataset, - wavelengths=None, - method="nearest", - cmap="jet", - vmin=0, - vmax=0.02, - ncols=1, + 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, ): """ @@ -56,10 +61,16 @@ def viz_pace( 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 @@ -81,30 +92,87 @@ def viz_pace( nrows = math.ceil(len(wavelengths) / ncols) - fig, axes = plt.subplots( - nrows=nrows, ncols=ncols, figsize=(6.4 * ncols, 4.8 * nrows), **kwargs - ) + 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: - 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() + 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() From 1b39239908441817dbf32162e2f71e029012f7d6 Mon Sep 17 00:00:00 2001 From: Qiusheng Wu Date: Sun, 5 May 2024 18:00:52 -0400 Subject: [PATCH 4/5] Add PACE notebook --- docs/examples/pace.ipynb | 121 +++++++++++++++++++++++++++++++++++++++ mkdocs.yml | 1 + pyproject.toml | 2 +- 3 files changed, 123 insertions(+), 1 deletion(-) create mode 100644 docs/examples/pace.ipynb diff --git a/docs/examples/pace.ipynb b/docs/examples/pace.ipynb new file mode 100644 index 0000000..248d3ea --- /dev/null +++ b/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 +} diff --git a/mkdocs.yml b/mkdocs.yml index 01753ae..52cc829 100644 --- a/mkdocs.yml +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index 00787f5..35d0289 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -34,7 +34,7 @@ all = [ ] extra = [ - "pandas", + "cartopy", ] From 8737375927012a662dd1cf59d7e390ced5634933 Mon Sep 17 00:00:00 2001 From: Qiusheng Wu Date: Sun, 5 May 2024 18:04:16 -0400 Subject: [PATCH 5/5] Add cartopy --- requirements_dev.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/requirements_dev.txt b/requirements_dev.txt index 7b255af..51807e1 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -2,6 +2,8 @@ black black[jupyter] build bump-my-version + +cartopy click codespell flake8