diff --git a/docs/notebooks/88_nasa_earth_data.ipynb b/docs/notebooks/88_nasa_earth_data.ipynb index 23f50348bb..0102630f55 100644 --- a/docs/notebooks/88_nasa_earth_data.ipynb +++ b/docs/notebooks/88_nasa_earth_data.ipynb @@ -124,7 +124,6 @@ { "cell_type": "code", "execution_count": null, - "id": "9e232fb9", "metadata": {}, "outputs": [], "source": [ diff --git a/docs/notebooks/89_image_array_viz.ipynb b/docs/notebooks/89_image_array_viz.ipynb new file mode 100644 index 0000000000..e73e9d7281 --- /dev/null +++ b/docs/notebooks/89_image_array_viz.ipynb @@ -0,0 +1,201 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[![image](https://jupyterlite.rtfd.io/en/latest/_static/badge.svg)](https://demo.leafmap.org/lab/index.html?path=notebooks/89_image_array_viz.ipynb)\n", + "[![image](https://studiolab.sagemaker.aws/studiolab.svg)](https://studiolab.sagemaker.aws/import/github/opengeos/leafmap/blob/master/examples/notebooks/89_image_array_viz.ipynb)\n", + "[![image](https://img.shields.io/badge/Open-Planetary%20Computer-black?style=flat&logo=microsoft)](https://pccompute.westeurope.cloudapp.azure.com/compute/hub/user-redirect/git-pull?repo=https://github.com/opengeos/leafmap&urlpath=lab/tree/leafmap/examples/notebooks/89_image_array_viz.ipynb&branch=master)\n", + "[![image](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/opengeos/leafmap/blob/master/examples/notebooks/89_image_array_viz.ipynb)\n", + "[![image](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/opengeos/leafmap/HEAD)\n", + "\n", + "**Visualizing in-memory raster datasets and image arrays**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# %pip install leafmap rasterio rioxarray" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import leafmap\n", + "import rasterio\n", + "import rioxarray\n", + "import xarray as xr" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Download two sample raster datasets." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "url1 = \"https://open.gishub.org/data/raster/landsat7.tif\"\n", + "url2 = \"https://open.gishub.org/data/raster/srtm90.tif\"\n", + "satellite = leafmap.download_file(url1, \"landsat7.tif\")\n", + "dem = leafmap.download_file(url2, \"srtm90.tif\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The Landsat image contains 3 bands: nir, red, and green. Let's calculate NDVI using the nir and red bands." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dataset = rasterio.open(satellite)\n", + "nir = dataset.read(1).astype(float)\n", + "red = dataset.read(2).astype(float)\n", + "ndvi = (nir - red) / (nir + red)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create an in-memory raster dataset from the NDVI array and use the projection and extent of the Landsat image." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ndvi_image = leafmap.array_to_image(ndvi, source=satellite)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Visualize the Landsat image and the NDVI image on the same map." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m = leafmap.Map()\n", + "m.add_raster(satellite, band=[1, 2, 3], nodata=-1, layer_name=\"Landsat 7\")\n", + "m.add_raster(ndvi_image, cmap=\"Greens\", layer_name=\"NDVI\")\n", + "m" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Use rioxarray to read raster datasets into xarray DataArrays." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds = rioxarray.open_rasterio(dem)\n", + "ds" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Classify the DEM into 2 elevation classes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "array = ds.sel(band=1)\n", + "masked_array = xr.where(array < 2000, 0, 1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create an in-memory raster dataset from the elevation class array and use the projection and extent of the DEM." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "image = leafmap.array_to_image(masked_array, source=dem)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Visualize the DEM and the elevation class image on the same map." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m = leafmap.Map()\n", + "m.add_raster(dem, cmap=\"terrain\", layer_name=\"DEM\")\n", + "m.add_raster(image, cmap=\"coolwarm\", layer_name=\"Classified DEM\")\n", + "m" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.11.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/tutorials.md b/docs/tutorials.md index e7b81e575d..dc5db866f2 100644 --- a/docs/tutorials.md +++ b/docs/tutorials.md @@ -102,6 +102,7 @@ 86. Adding markers to the map ([notebook](https://leafmap.org/notebooks/86_add_markers)) 87. Cloud-based geoprocessing with Actinia ([notebook](https://leafmap.org/notebooks/87_actinia)) 88. Searching and downloading NASA Earth science data products ([notebook](https://leafmap.org/notebooks/88_nasa_earth_data)) +89. Visualizing in-memory raster datasets and image arrays ([notebook](https://leafmap.org/notebooks/89_image_array_viz)) ## Demo diff --git a/examples/README.md b/examples/README.md index 90a9ec47e1..d6f9e1715c 100644 --- a/examples/README.md +++ b/examples/README.md @@ -109,6 +109,7 @@ 86. Adding markers to the map ([notebook](https://leafmap.org/notebooks/86_add_markers)) 87. Cloud-based geoprocessing with Actinia ([notebook](https://leafmap.org/notebooks/87_actinia)) 88. Searching and downloading NASA Earth science data products ([notebook](https://leafmap.org/notebooks/88_nasa_earth_data)) +89. Visualizing in-memory raster datasets and image arrays ([notebook](https://leafmap.org/notebooks/89_image_array_viz)) ## Demo diff --git a/examples/notebooks/88_nasa_earth_data.ipynb b/examples/notebooks/88_nasa_earth_data.ipynb index 23f50348bb..0102630f55 100644 --- a/examples/notebooks/88_nasa_earth_data.ipynb +++ b/examples/notebooks/88_nasa_earth_data.ipynb @@ -124,7 +124,6 @@ { "cell_type": "code", "execution_count": null, - "id": "9e232fb9", "metadata": {}, "outputs": [], "source": [ diff --git a/examples/notebooks/89_image_array_viz.ipynb b/examples/notebooks/89_image_array_viz.ipynb new file mode 100644 index 0000000000..e73e9d7281 --- /dev/null +++ b/examples/notebooks/89_image_array_viz.ipynb @@ -0,0 +1,201 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[![image](https://jupyterlite.rtfd.io/en/latest/_static/badge.svg)](https://demo.leafmap.org/lab/index.html?path=notebooks/89_image_array_viz.ipynb)\n", + "[![image](https://studiolab.sagemaker.aws/studiolab.svg)](https://studiolab.sagemaker.aws/import/github/opengeos/leafmap/blob/master/examples/notebooks/89_image_array_viz.ipynb)\n", + "[![image](https://img.shields.io/badge/Open-Planetary%20Computer-black?style=flat&logo=microsoft)](https://pccompute.westeurope.cloudapp.azure.com/compute/hub/user-redirect/git-pull?repo=https://github.com/opengeos/leafmap&urlpath=lab/tree/leafmap/examples/notebooks/89_image_array_viz.ipynb&branch=master)\n", + "[![image](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/opengeos/leafmap/blob/master/examples/notebooks/89_image_array_viz.ipynb)\n", + "[![image](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/opengeos/leafmap/HEAD)\n", + "\n", + "**Visualizing in-memory raster datasets and image arrays**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# %pip install leafmap rasterio rioxarray" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import leafmap\n", + "import rasterio\n", + "import rioxarray\n", + "import xarray as xr" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Download two sample raster datasets." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "url1 = \"https://open.gishub.org/data/raster/landsat7.tif\"\n", + "url2 = \"https://open.gishub.org/data/raster/srtm90.tif\"\n", + "satellite = leafmap.download_file(url1, \"landsat7.tif\")\n", + "dem = leafmap.download_file(url2, \"srtm90.tif\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The Landsat image contains 3 bands: nir, red, and green. Let's calculate NDVI using the nir and red bands." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dataset = rasterio.open(satellite)\n", + "nir = dataset.read(1).astype(float)\n", + "red = dataset.read(2).astype(float)\n", + "ndvi = (nir - red) / (nir + red)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create an in-memory raster dataset from the NDVI array and use the projection and extent of the Landsat image." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ndvi_image = leafmap.array_to_image(ndvi, source=satellite)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Visualize the Landsat image and the NDVI image on the same map." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m = leafmap.Map()\n", + "m.add_raster(satellite, band=[1, 2, 3], nodata=-1, layer_name=\"Landsat 7\")\n", + "m.add_raster(ndvi_image, cmap=\"Greens\", layer_name=\"NDVI\")\n", + "m" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Use rioxarray to read raster datasets into xarray DataArrays." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds = rioxarray.open_rasterio(dem)\n", + "ds" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Classify the DEM into 2 elevation classes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "array = ds.sel(band=1)\n", + "masked_array = xr.where(array < 2000, 0, 1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create an in-memory raster dataset from the elevation class array and use the projection and extent of the DEM." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "image = leafmap.array_to_image(masked_array, source=dem)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Visualize the DEM and the elevation class image on the same map." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m = leafmap.Map()\n", + "m.add_raster(dem, cmap=\"terrain\", layer_name=\"DEM\")\n", + "m.add_raster(image, cmap=\"coolwarm\", layer_name=\"Classified DEM\")\n", + "m" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.11.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/leafmap/common.py b/leafmap/common.py index adeaf949bc..bde5af4722 100644 --- a/leafmap/common.py +++ b/leafmap/common.py @@ -2848,6 +2848,7 @@ def get_local_tile_layer( """ from osgeo import gdal + import rasterio # ... and suppress errors gdal.PushErrorHandler("CPLQuietErrorHandler") @@ -2896,6 +2897,11 @@ def get_local_tile_layer( raise ValueError("The source path does not exist.") else: source = github_raw_url(source) + elif isinstance(source, TileClient) or isinstance( + source, rasterio.io.DatasetReader + ): + pass + else: raise ValueError("The source must either be a string or TileClient") @@ -2911,7 +2917,11 @@ def get_local_tile_layer( else: layer_name = "LocalTile_" + random_string(3) - tile_client = TileClient(source, port=port, debug=debug) + if isinstance(source, str) or isinstance(source, rasterio.io.DatasetReader): + tile_client = TileClient(source, port=port, debug=debug) + + else: + tile_client = source if "cmap" not in kwargs: kwargs["cmap"] = palette @@ -10234,6 +10244,131 @@ def install_package(package): process.wait() +def array_to_memory_file( + array, + source: str = None, + dtype: str = None, + compress: str = "deflate", + transpose: bool = True, + cellsize: float = None, + crs: str = None, + transform: tuple = None, + driver="COG", + **kwargs, +): + """Convert a NumPy array to a memory file. + + Args: + array (numpy.ndarray): The input NumPy array. + source (str, optional): Path to the source file to extract metadata from. Defaults to None. + dtype (str, optional): The desired data type of the array. Defaults to None. + compress (str, optional): The compression method for the output file. Defaults to "deflate". + transpose (bool, optional): Whether to transpose the array from (bands, rows, columns) to (rows, columns, bands). Defaults to True. + cellsize (float, optional): The cell size of the array if source is not provided. Defaults to None. + crs (str, optional): The coordinate reference system of the array if source is not provided. Defaults to None. + transform (tuple, optional): The affine transformation matrix if source is not provided. Defaults to None. + driver (str, optional): The driver to use for creating the output file, such as 'GTiff'. Defaults to "COG". + **kwargs: Additional keyword arguments to be passed to the rasterio.open() function. + + Returns: + rasterio.DatasetReader: The rasterio dataset reader object for the converted array. + """ + import rasterio + import numpy as np + import xarray as xr + + if isinstance(array, xr.DataArray): + array = array.values + + if array.ndim == 3 and transpose: + array = np.transpose(array, (1, 2, 0)) + + if source is not None: + with rasterio.open(source) as src: + crs = src.crs + transform = src.transform + if compress is None: + compress = src.compression + else: + if cellsize is None: + raise ValueError("cellsize must be provided if source is not provided") + if crs is None: + raise ValueError( + "crs must be provided if source is not provided, such as EPSG:3857" + ) + + if "transform" not in kwargs: + # Define the geotransformation parameters + xmin, ymin, xmax, ymax = ( + 0, + 0, + cellsize * array.shape[1], + cellsize * array.shape[0], + ) + transform = rasterio.transform.from_bounds( + xmin, ymin, xmax, ymax, array.shape[1], array.shape[0] + ) + else: + transform = kwargs["transform"] + + if dtype is None: + # Determine the minimum and maximum values in the array + min_value = np.min(array) + max_value = np.max(array) + # Determine the best dtype for the array + if min_value >= 0 and max_value <= 1: + dtype = np.float32 + elif min_value >= 0 and max_value <= 255: + dtype = np.uint8 + elif min_value >= -128 and max_value <= 127: + dtype = np.int8 + elif min_value >= 0 and max_value <= 65535: + dtype = np.uint16 + elif min_value >= -32768 and max_value <= 32767: + dtype = np.int16 + else: + dtype = np.float64 + + # Convert the array to the best dtype + array = array.astype(dtype) + + # Define the GeoTIFF metadata + metadata = { + "driver": driver, + "height": array.shape[0], + "width": array.shape[1], + "dtype": array.dtype, + "crs": crs, + "transform": transform, + } + + if array.ndim == 2: + metadata["count"] = 1 + elif array.ndim == 3: + metadata["count"] = array.shape[2] + if compress is not None: + metadata["compress"] = compress + + metadata.update(**kwargs) + + # Create a new memory file and write the array to it + memory_file = rasterio.MemoryFile() + dst = memory_file.open(**metadata) + + if array.ndim == 2: + dst.write(array, 1) + elif array.ndim == 3: + for i in range(array.shape[2]): + dst.write(array[:, :, i], i + 1) + + dst.close() + + # Read the dataset from memory + dataset_reader = rasterio.open(dst.name, mode="r") + + return dataset_reader + + def array_to_image( array, output: str = None, @@ -10241,8 +10376,9 @@ def array_to_image( dtype: str = None, compress: str = "deflate", transpose: bool = True, - resolution: float = None, + cellsize: float = None, crs: str = None, + driver: str = "COG", **kwargs, ) -> str: """Save a NumPy array as a GeoTIFF using the projection information from an existing GeoTIFF file. @@ -10254,20 +10390,23 @@ def array_to_image( dtype (np.dtype, optional): The data type of the output array. Defaults to None. compress (str, optional): The compression method. Can be one of the following: "deflate", "lzw", "packbits", "jpeg". Defaults to "deflate". transpose (bool, optional): Whether to transpose the array from (bands, rows, columns) to (rows, columns, bands). Defaults to True. - resolution (float, optional): The resolution of the output image in meters. Defaults to None. + cellsize (float, optional): The resolution of the output image in meters. Defaults to None. crs (str, optional): The CRS of the output image. Defaults to None. + driver (str, optional): The driver to use for creating the output file, such as 'GTiff'. Defaults to "COG". + **kwargs: Additional keyword arguments to be passed to the rasterio.open() function. """ import numpy as np import rasterio - import tempfile + + if output is None: + return array_to_memory_file( + array, source, dtype, compress, transpose, cellsize, crs, driver, **kwargs + ) if array.ndim == 3 and transpose: array = np.transpose(array, (1, 2, 0)) - if output is None: - output = tempfile.NamedTemporaryFile(suffix=".tif", delete=False).name - out_dir = os.path.dirname(os.path.abspath(output)) if not os.path.exists(out_dir): os.makedirs(out_dir) @@ -10282,7 +10421,7 @@ def array_to_image( if compress is None: compress = src.compression else: - if resolution is None: + if cellsize is None: raise ValueError("resolution must be provided if source is not provided") if crs is None: raise ValueError( @@ -10294,8 +10433,8 @@ def array_to_image( xmin, ymin, xmax, ymax = ( 0, 0, - resolution * array.shape[1], - resolution * array.shape[0], + cellsize * array.shape[1], + cellsize * array.shape[0], ) transform = rasterio.transform.from_bounds( xmin, ymin, xmax, ymax, array.shape[1], array.shape[0] @@ -10303,11 +10442,10 @@ def array_to_image( else: transform = kwargs["transform"] - # Determine the minimum and maximum values in the array - min_value = np.min(array) - max_value = np.max(array) - if dtype is None: + # Determine the minimum and maximum values in the array + min_value = np.min(array) + max_value = np.max(array) # Determine the best dtype for the array if min_value >= 0 and max_value <= 1: dtype = np.float32 @@ -10326,33 +10464,23 @@ def array_to_image( array = array.astype(dtype) # Define the GeoTIFF metadata + metadata = { + "driver": driver, + "height": array.shape[0], + "width": array.shape[1], + "dtype": array.dtype, + "crs": crs, + "transform": transform, + } + if array.ndim == 2: - metadata = { - "driver": "GTiff", - "height": array.shape[0], - "width": array.shape[1], - "count": 1, - "dtype": array.dtype, - "crs": crs, - "transform": transform, - **kwargs, - } + metadata["count"] = 1 elif array.ndim == 3: - metadata = { - "driver": "GTiff", - "height": array.shape[0], - "width": array.shape[1], - "count": array.shape[2], - "dtype": array.dtype, - "crs": crs, - "transform": transform, - **kwargs, - } - + metadata["count"] = array.shape[2] if compress is not None: metadata["compress"] = compress - else: - raise ValueError("Array must be 2D or 3D.") + + metadata.update(**kwargs) # Create a new GeoTIFF file and write the array to it with rasterio.open(output, "w", **metadata) as dst: diff --git a/leafmap/leafmap.py b/leafmap/leafmap.py index db178bdde5..fcd4bc5080 100644 --- a/leafmap/leafmap.py +++ b/leafmap/leafmap.py @@ -72,7 +72,7 @@ def __init__(self, **kwargs): self.sandbox_path = None if "height" not in kwargs: - self.layout.height = "500px" + self.layout.height = "600px" else: if isinstance(kwargs["height"], int): kwargs["height"] = str(kwargs["height"]) + "px" @@ -93,7 +93,7 @@ def __init__(self, **kwargs): self.add(ipyleaflet.FullScreenControl()) if "search_control" not in kwargs: - kwargs["search_control"] = True + kwargs["search_control"] = False if kwargs["search_control"]: url = "https://nominatim.openstreetmap.org/search?format=json&q={s}" search_control = ipyleaflet.SearchControl( @@ -180,7 +180,7 @@ def handle_draw(target, action, geo_json): draw_control.on_draw(handle_draw) if "measure_control" not in kwargs: - kwargs["measure_control"] = True + kwargs["measure_control"] = False if kwargs["measure_control"]: self.add(ipyleaflet.MeasureControl(position="topleft")) diff --git a/mkdocs.yml b/mkdocs.yml index fd4c191231..2f121163bb 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -224,3 +224,4 @@ nav: - notebooks/86_add_markers.ipynb - notebooks/87_actinia.ipynb - notebooks/88_nasa_earth_data.ipynb + - notebooks/89_image_array_viz.ipynb