diff --git a/docs/examples/pace.ipynb b/docs/examples/pace.ipynb index 1a9c1aa9..1afc657c 100644 --- a/docs/examples/pace.ipynb +++ b/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": { diff --git a/hypercoast/emit.py b/hypercoast/emit.py index 18632c4a..22bde848 100644 --- a/hypercoast/emit.py +++ b/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: diff --git a/hypercoast/hypercoast.py b/hypercoast/hypercoast.py index 4e803f8e..b499e873 100644 --- a/hypercoast/hypercoast.py +++ b/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, diff --git a/hypercoast/pace.py b/hypercoast/pace.py index b764f5b9..5c22cad8 100644 --- a/hypercoast/pace.py +++ b/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,6 +237,13 @@ 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 @@ -244,24 +251,65 @@ def grid_pace(dataset, wavelengths, method="nearest", **kwargs): 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)