Add support for multi-band visualization (#24)

* Add support for multi-band visualization

* Add viz examples to pace notebook

docs/examples/pace.ipynb

@@ -131,6 +131,63 @@
     "longitude = (-92, -91.055)\n",
     "hypercoast.filter_pace(dataset, latitude, longitude, return_plot=True)"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Single-band visualization."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m = hypercoast.Map()\n",
+    "m.add_basemap(\"Hybrid\")\n",
+    "wavelengths = [450]\n",
+    "m.add_pace(dataset, wavelengths, colormap=\"jet\", vmin=0, vmax=0.02, layer_name=\"PACE\")\n",
+    "m.add_colormap(cmap=\"jet\", vmin=0, vmax=0.02, label=\"Reflectance\")\n",
+    "m"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "![](https://i.imgur.com/2rBzam8.png)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Multiple-band visualization."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m = hypercoast.Map()\n",
+    "m.add_basemap(\"Hybrid\")\n",
+    "wavelengths = [450, 550, 650]\n",
+    "m.add_pace(\n",
+    "    dataset, wavelengths, indexes=[3, 2, 1], vmin=0, vmax=0.02, layer_name=\"PACE\"\n",
+    ")\n",
+    "m"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "![](https://i.imgur.com/40gX4og.png)"
+   ]
   }
  ],
  "metadata": {

hypercoast/emit.py

@@ -226,15 +226,15 @@ def emit_to_netcdf(data, output, **kwargs):
     ds_geo.to_netcdf(output, **kwargs)
 
 
-def emit_to_image(data, output=None, wavelengths=None, method="nearest", **kwargs):
+def emit_to_image(data, wavelengths=None, method="nearest", output=None, **kwargs):
     """
     Converts an EMIT dataset to an image.
 
     Args:
         data (xarray.Dataset or str): The dataset containing the EMIT data or the file path to the dataset.
-        output (str, optional): The file path where the image will be saved. If None, the image will be returned as a PIL Image object. Defaults to None.
         wavelengths (array-like, optional): The specific wavelengths to select. If None, all wavelengths are selected. Defaults to None.
         method (str, optional): The method to use for data selection. Defaults to "nearest".
+        output (str, optional): The file path where the image will be saved. If None, the image will be returned as a PIL Image object. Defaults to None.
         **kwargs: Additional keyword arguments to be passed to `leafmap.array_to_image`.
 
     Returns:

hypercoast/hypercoast.py

@@ -192,14 +192,15 @@ def add_pace(
         wavelengths=None,
         indexes=None,
         colormap="jet",
-        vmin=0,
-        vmax=0.02,
+        vmin=None,
+        vmax=None,
         nodata=np.nan,
         attribution=None,
         layer_name="PACE",
         zoom_to_layer=True,
         visible=True,
         method="nearest",
+        gridded=False,
         array_args={},
         **kwargs,
     ):
@@ -229,10 +230,15 @@ def add_pace(
 
             source = read_pace(source)
 
-        source = grid_pace(source, wavelengths, method=method)
+        image = pace_to_image(
+            source, wavelengths=wavelengths, method=method, gridded=gridded
+        )
+
+        if isinstance(wavelengths, list) and len(wavelengths) > 1:
+            colormap = None
 
         self.add_raster(
-            source,
+            image,
             indexes=indexes,
             colormap=colormap,
             vmin=vmin,

hypercoast/pace.py

@@ -221,7 +221,7 @@ def filter_pace(dataset, latitude, longitude, drop=True, return_plot=False, **kw
         return da_filtered
 
 
-def grid_pace(dataset, wavelengths, method="nearest", **kwargs):
+def grid_pace(dataset, wavelengths=None, method="nearest", **kwargs):
     """
     Grids a PACE dataset based on latitude and longitude.
 
@@ -237,31 +237,79 @@ def grid_pace(dataset, wavelengths, method="nearest", **kwargs):
     """
     from scipy.interpolate import griddata
 
+    if wavelengths is None:
+        wavelengths = dataset.coords["wavelength"].values[0]
+
+    # Ensure wavelengths is a list
+    if not isinstance(wavelengths, list):
+        wavelengths = [wavelengths]
+
     lat = dataset.latitude
     lon = dataset.longitude
 
     grid_lat = np.linspace(lat.min(), lat.max(), lat.shape[0])
     grid_lon = np.linspace(lon.min(), lon.max(), lon.shape[1])
     grid_lon_2d, grid_lat_2d = np.meshgrid(grid_lon, grid_lat)
 
-    data = dataset.sel(wavelength=wavelengths)["Rrs"]
-    gridded_data = griddata(
-        (lat.data.flatten(), lon.data.flatten()),
-        data.data.flatten(),
-        (grid_lat_2d, grid_lon_2d),
-        method=method,
-    )
+    gridded_data_dict = {}
+    for wavelength in wavelengths:
+        data = dataset.sel(wavelength=wavelength)["Rrs"]
+        gridded_data = griddata(
+            (lat.data.flatten(), lon.data.flatten()),
+            data.data.flatten(),
+            (grid_lat_2d, grid_lon_2d),
+            method=method,
+        )
+        gridded_data_dict[wavelength] = gridded_data
+
+    # Create a 3D array with dimensions latitude, longitude, and wavelength
+    gridded_data_3d = np.dstack(list(gridded_data_dict.values()))
 
     dataset2 = xr.Dataset(
-        {"Rrs": (("latitude", "longitude"), gridded_data)},
+        {"Rrs": (("latitude", "longitude", "wavelength"), gridded_data_3d)},
         coords={
             "latitude": ("latitude", grid_lat),
             "longitude": ("longitude", grid_lon),
+            "wavelength": ("wavelength", list(gridded_data_dict.keys())),
         },
         **kwargs,
     )
 
-    data = dataset2["Rrs"]
+    dataset2["Rrs"].rio.write_crs("EPSG:4326", inplace=True)
+
+    return dataset2
+
+
+def pace_to_image(
+    dataset, wavelengths=None, method="nearest", gridded=False, output=None, **kwargs
+):
+    """
+    Converts an PACE dataset to an image.
+
+    Args:
+        dataset (xarray.Dataset or str): The dataset containing the EMIT data or the file path to the dataset.
+        wavelengths (array-like, optional): The specific wavelengths to select. If None, all wavelengths are selected. Defaults to None.
+        method (str, optional): The method to use for data interpolation. Defaults to "nearest".
+        gridded (bool, optional): Whether the dataset is a gridded dataset. Defaults to False,
+        output (str, optional): The file path where the image will be saved. If None, the image will be returned as a PIL Image object. Defaults to None.
+        **kwargs: Additional keyword arguments to be passed to `leafmap.array_to_image`.
+
+    Returns:
+        rasterio.Dataset or None: The image converted from the dataset. If `output` is provided, the image will be saved to the specified file and the function will return None.
+    """
+    from leafmap import array_to_image
+
+    if isinstance(dataset, str):
+        dataset = read_pace(dataset, wavelengths=wavelengths, method="nearest")
+
+    if wavelengths is not None:
+        dataset = dataset.sel(wavelength=wavelengths, method="nearest")
+
+    if not gridded:
+        grid = grid_pace(dataset, wavelengths=wavelengths, method=method)
+    else:
+        grid = dataset
+    data = grid["Rrs"]
     data.rio.write_crs("EPSG:4326", inplace=True)
 
-    return data
+    return array_to_image(data, transpose=False, output=output, **kwargs)