Add function to plot composites with ipyleaflet (#98)

Generates an image with matplotlib and overlays on a leaflet map.
Requires pyproj to get the image boundaries in lay, lon. Add ipyleaflet,
pyproj, and matplotlib as dependencies. Use the new function in a few
examples.

Fixes #94

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;
+}

doc/api/index.rst

@@ -26,6 +26,7 @@ Visualization
     composite
     equalize_histogram
     adjust_l1_colors
+    plot_composite_leaflet
 
 Indices
 -------

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
+<https://ipyleaflet.readthedocs.io/en/latest/>`__ 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.
 

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

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")
 
 
 ----

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

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

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]

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__

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