From 1d5ac6a2acf3b147c0f45df7361531a9214785fe Mon Sep 17 00:00:00 2001 From: Leonardo Uieda Date: Wed, 28 Aug 2024 15:52:15 -0300 Subject: [PATCH 1/2] Add function to plot composites with ipyleaflet --- doc/_static/custom.css | 5 ++ doc/api/index.rst | 1 + doc/composites.rst | 12 +++- doc/conf.py | 1 + doc/index.rst | 16 ++---- env/requirements-docs.txt | 1 - environment.yml | 4 +- pyproject.toml | 3 + xlandsat/__init__.py | 1 + xlandsat/_leaflet.py | 116 ++++++++++++++++++++++++++++++++++++++ 10 files changed, 146 insertions(+), 14 deletions(-) create mode 100644 xlandsat/_leaflet.py diff --git a/doc/_static/custom.css b/doc/_static/custom.css index ebfd14f..a71d1af 100644 --- a/doc/_static/custom.css +++ b/doc/_static/custom.css @@ -32,3 +32,8 @@ div.jupyter_container .cell_output { .jupyter_container div.code_cell pre { padding: 10px; } + +/* Disable borders in darkmode since it messes up the ipyleaflet plots */ +html[data-theme="dark"] .bd-content img:not(.only-dark):not(.dark-light) { + background: none; +} diff --git a/doc/api/index.rst b/doc/api/index.rst index 6ce698c..1bc576f 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -26,6 +26,7 @@ Visualization composite equalize_histogram adjust_l1_colors + plot_composite_leaflet Indices ------- diff --git a/doc/composites.rst b/doc/composites.rst index 55df601..9f02d0a 100644 --- a/doc/composites.rst +++ b/doc/composites.rst @@ -84,8 +84,16 @@ instead: plt.show() -Well, this looks bad because some very bright pixels in the city are making the -majority of the other pixels have only a small share of the full range of +It's also possible to add a composite to an interactive `ipyleaflet +`__ map using +:func:`xlandsat.plot_composite_leaflet`: + +.. jupyter-execute:: + + xls.plot_composite_leaflet(rgb) + +This composite looks bad because some very bright pixels in the city are making +the majority of the other pixels have only a small share of the full range of available values. This can be mitigated by rescaling the intensity of the image to a smaller range of reflectance values. diff --git a/doc/conf.py b/doc/conf.py index 087bb93..e122ec6 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -41,6 +41,7 @@ "pooch": ("https://www.fatiando.org/pooch/latest/", None), "matplotlib": ("https://matplotlib.org/stable/", None), "scipy": ("https://docs.scipy.org/doc/scipy/", None), + "ipyleaflet": ("https://ipyleaflet.readthedocs.io/en/latest/", None), } # Autosummary pages will be generated by sphinx-autogen instead of sphinx-build diff --git a/doc/index.rst b/doc/index.rst index bc9e791..0f240fa 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -36,24 +36,20 @@ Here's a quick example: .. jupyter-execute:: import xlandsat as xls - import matplotlib.pyplot as plt # Download a sample Landsat 9 scene in EarthExplorer format path_to_scene_file = xls.datasets.fetch_manaus() - # Load the data from the file into an xarray.Dataset scene = xls.load_scene(path_to_scene_file) + # Display the scene and included metadata + scene + +.. jupyter-execute:: # Make an RGB composite as an xarray.DataArray rgb = xls.composite(scene, rescale_to=[0.02, 0.2]) - - # Plot the composite using xarray's plotting machinery - rgb.plot.imshow() - - # Annotate the plot with the rich metadata xlandsat adds to the scene - plt.title(f"{rgb.attrs['title']}\n{rgb.attrs['long_name']}") - plt.axis("scaled") - plt.show() + # Plot the composite on an interactive Leaflet map + xls.plot_composite_leaflet(rgb, height="400px") ---- diff --git a/env/requirements-docs.txt b/env/requirements-docs.txt index 8933143..f181e39 100644 --- a/env/requirements-docs.txt +++ b/env/requirements-docs.txt @@ -4,5 +4,4 @@ sphinx-book-theme==1.1.* sphinx-copybutton==0.5.* sphinx-design==0.5.* jupyter-sphinx==0.5.* -matplotlib ipykernel diff --git a/environment.yml b/environment.yml index 63299a9..13d4f20 100644 --- a/environment.yml +++ b/environment.yml @@ -12,6 +12,9 @@ dependencies: - xarray - scikit-image - pooch + - pyproj + - ipyleaflet + - matplotlib # Build - build - twine @@ -25,7 +28,6 @@ dependencies: - sphinx-copybutton==0.5.* - sphinx-design==0.5.* - jupyter-sphinx==0.5.* - - matplotlib - ipykernel # Style - black diff --git a/pyproject.toml b/pyproject.toml index eec2097..4a0dffd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,6 +31,9 @@ dependencies = [ "xarray>=2022.6.0", "scikit-image>=0.20", "pooch>=1.3.0", + "pyproj>=3.3.0", + "ipyleaflet>=0.18", + "matplotlib>=3.5", ] [project.urls] diff --git a/xlandsat/__init__.py b/xlandsat/__init__.py index 87756d7..1764be7 100644 --- a/xlandsat/__init__.py +++ b/xlandsat/__init__.py @@ -7,5 +7,6 @@ from ._indices import nbr, ndvi from ._interpolation import interpolate_missing from ._io import load_panchromatic, load_scene, save_scene +from ._leaflet import plot_composite_leaflet from ._pansharpen import pansharpen from ._version import __version__ diff --git a/xlandsat/_leaflet.py b/xlandsat/_leaflet.py new file mode 100644 index 0000000..c5bec59 --- /dev/null +++ b/xlandsat/_leaflet.py @@ -0,0 +1,116 @@ +# Copyright (c) 2022 The xlandsat developers. +# Distributed under the terms of the MIT License. +# SPDX-License-Identifier: MIT +""" +Functions for plotting data with ipyleaflet +""" +import base64 +import io + +import ipyleaflet + +# Dependecy of ipyleaflet +import ipywidgets +import matplotlib.pyplot as plt +import pyproj + + +def _scene_boundaries_geodetic(scene): + """ + Determine the boundaries of the scene in geodetic longitude and latitude + + Uses pyproj to project the UTM coordinates of the scene and get the + geographic bounding box. + + Returns + ------- + [w, e, s, n] : list of floats + The west, east, south, north boundaries of the scene in degrees. + """ + projection = pyproj.Proj( + proj="utm", zone=scene.attrs["utm_zone"], ellps=scene.attrs["ellipsoid"] + ) + west, south = projection( + scene.easting.min().values, scene.northing.min().values, inverse=True + ) + east, north = projection( + scene.easting.max().values, scene.northing.max().values, inverse=True + ) + return (west, east, south, north) + + +def plot_composite_leaflet(composite, dpi=70, leaflet_map=None, height="600px"): + """ + Display a composite as an image overlay in an interactive HTML map + + Adds the composite to a Leaflet.js map, which can be displayed in a Jupyter + notebook or HTML page. By default, adds a control widget for the opacity of + the image overlay. + + Parameters + ---------- + composite : :class:`xarray.DataArray` + A composite, as generated by :func:`xlandsat.composite`. + dpi : int + The dots-per-inch resolution of the image. + leaflet_map : :class:`ipyleaflet.Map` + A Leaflet map instance to which the image overlay will be added. If + None (default), a new map will be created. Pass an existing map to add + the overlay to it. + height : str + The height of the map which is embedded in the HTML. Should contain the + proper CSS units (px, em, rem, etc). + + Returns + ------- + leaflet_map : :class:`ipyleaflet.Map` + The map with the image overlay and opacity controls added to it. + """ + west, east, south, north = _scene_boundaries_geodetic(composite) + center = (0.5 * (north + south), 0.5 * (east + west)) + bounds = ((south, west), (north, east)) + # Create a plot of the composite with no decoration + fig, ax = plt.subplots(1, 1, layout="constrained") + composite.plot.imshow(ax=ax, add_labels=False) + ax.axis("off") + ax.set_aspect("equal") + # Save the PNG to an in-memory buffer + png = io.BytesIO() + fig.savefig(png, bbox_inches="tight", dpi=dpi, pad_inches=0, transparent=True) + plt.close(fig) + # Create the image overlay with the figure as base64 encoded png + image_overlay = ipyleaflet.ImageOverlay( + url=f"data:image/png;base64,{base64.b64encode(png.getvalue()).decode()}", + bounds=bounds, + ) + image_overlay.name = ( + f"{composite.attrs["long_name"].title()} | {composite.attrs["title"]}" + ) + # Create a map if one wasn't given + if leaflet_map is None: + leaflet_map = ipyleaflet.Map( + center=center, + scroll_wheel_zoom=True, + layout={"height": height}, + ) + leaflet_map.add(ipyleaflet.ScaleControl(position="bottomleft")) + leaflet_map.add(ipyleaflet.LayersControl(position="bottomright")) + leaflet_map.add(ipyleaflet.FullScreenControl()) + leaflet_map.fit_bounds(bounds) + # Add a widget to control the opacity of the image + opacity_slider = ipywidgets.FloatSlider( + description="Opacity:", + min=0, + max=1, + step=0.1, + value=1, + readout_format=".1f", + style={"description_width": "initial"}, + layout={"margin": "0 0 0 0.5rem"}, + ) + ipywidgets.jslink((opacity_slider, "value"), (image_overlay, "opacity")) + leaflet_map.add( + ipyleaflet.WidgetControl(widget=opacity_slider, position="topright") + ) + leaflet_map.add(image_overlay) + return leaflet_map From 20494265c7c011995f142648780c84c6b3affc41 Mon Sep 17 00:00:00 2001 From: Leonardo Uieda Date: Wed, 28 Aug 2024 19:59:21 -0300 Subject: [PATCH 2/2] Use single quotes in docstring --- xlandsat/_leaflet.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xlandsat/_leaflet.py b/xlandsat/_leaflet.py index c5bec59..2075f9b 100644 --- a/xlandsat/_leaflet.py +++ b/xlandsat/_leaflet.py @@ -84,7 +84,7 @@ def plot_composite_leaflet(composite, dpi=70, leaflet_map=None, height="600px"): bounds=bounds, ) image_overlay.name = ( - f"{composite.attrs["long_name"].title()} | {composite.attrs["title"]}" + f"{composite.attrs['long_name'].title()} | {composite.attrs['title']}" ) # Create a map if one wasn't given if leaflet_map is None: