Add NDVI calculation (#92)

Add a function that calculates NDVI for a given scene. It's better than
having to remember the equations all the time. The function can also
take care of adding some metadata and copying the metadata from the
original scene. The function allows customizing the names of the red and
NIR bands.

doc/api/index.rst

@@ -11,22 +11,37 @@ Input and output
 ----------------
 
 .. autosummary::
-   :toctree: generated/
+    :toctree: generated/
 
     load_scene
     save_scene
     load_panchromatic
 
-Processing
-----------
+Visualization
+-------------
 
 .. autosummary::
-   :toctree: generated/
+    :toctree: generated/
 
     composite
-    pansharpen
     equalize_histogram
     adjust_l1_colors
+
+Indices
+-------
+
+.. autosummary::
+    :toctree: generated/
+
+    ndvi
+
+Processing
+----------
+
+.. autosummary::
+    :toctree: generated/
+
+    pansharpen
     interpolate_missing
 
 Sample datasets

doc/index.rst

@@ -163,10 +163,10 @@ Here's a quick example:
 
     io.rst
     composites.rst
+    equalize-histogram.rst
     indices.rst
-    pansharpen.rst
     missing-values.rst
-    equalize-histogram.rst
+    pansharpen.rst
     plot-overlay.rst
 
 .. toctree::

doc/indices.rst

@@ -10,9 +10,26 @@ differentiate between them.
 Many of these indices can be calculated with simple arithmetic operations.
 So now that our data are in :class:`xarray.Dataset`'s, it's fairly easy to
 calculate them.
-As an example, we'll use two example scenes from before and after the
+Some commonly used indices are provided as functions in :mod:`xlandsat` but
+any other index can be calculated using the power of :mod:`xarray`.
+
+Here, we'll see some examples of indices that can be calculated.
+First, we must import the required packages for the examples below.
+
+.. jupyter-execute::
+
+    import xlandsat as xls
+    import matplotlib.pyplot as plt
+    import numpy as np
+
+NDVI
+----
+
+As an example, we'll use two scenes from before and after the
 `Brumadinho tailings dam disaster <https://en.wikipedia.org/wiki/Brumadinho_dam_disaster>`__
-to try to image and quantify the total area flooded by the damn collapse.
+to try to image and quantify the total area flooded by the damn collapse using
+`NDVI <https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index>`__,
+which is available in :func:`xlandsat.ndvi`.
 
 .. admonition:: Trigger warning
     :class: warning
@@ -25,24 +42,18 @@ to try to image and quantify the total area flooded by the damn collapse.
     **Some readers may find this topic disturbing and may not wish to read
     futher.**
 
-First, we must import the required packages, download our two sample scenes,
-and load them with :func:`xlandsat.load_scene`:
+First, we must download the sample data and load the two scenes with
+:func:`xlandsat.load_scene`:
 
 .. jupyter-execute::
 
-    import xlandsat as xls
-    import matplotlib.pyplot as plt
-
-
     path_before = xls.datasets.fetch_brumadinho_before()
     path_after = xls.datasets.fetch_brumadinho_after()
-
     before = xls.load_scene(path_before)
     after = xls.load_scene(path_after)
     after
 
-Let's make RGB composites to get a sense of what these
-two scenes contain:
+Let's make RGB composites to get a sense of what these two scenes contain:
 
 .. jupyter-execute::
 
@@ -65,42 +76,23 @@ as it was contaminated by the mud flow.
 
      See :ref:`composites` for more information on what we did above.
 
-NDVI
-----
-
-We can calculate the
-`NDVI <https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index>`__
-for these scenes to see if we can isolate the effect of the flood following the
-dam collapse.
+We can calculate the NDVI for these scenes to see if we can isolate the effect
+of the flood following the dam collapse.
 NDVI highlights vegetation, which we assume will have decreased in the after
 scene due to the flood.
-NDVI is defined as:
-
-.. math::
-
-    NDVI = \dfrac{NIR - Red}{NIR + Red}
-
-which we can calculate with xarray as:
-
-.. jupyter-execute::
-
-    ndvi_before = (before.nir - before.red) / (before.nir + before.red)
-    ndvi_before
-
-Now we can do the same for the after scene:
+Calculating the NDVI is as simple as calling :func:`xlandsat.ndvi` with the
+scene as an argument:
 
 .. jupyter-execute::
 
-    ndvi_after = (after.nir - after.red) / (after.nir + after.red)
+    ndvi_before = xls.ndvi(before)
+    ndvi_after = xls.ndvi(after)
     ndvi_after
 
-And add some metadata for xarray to find when making plots:
+And add some extra metadata for xarray to find when making plots:
 
 .. jupyter-execute::
 
-    for ndvi in [ndvi_before, ndvi_after]:
-        ndvi.attrs["long_name"] = "normalized difference vegetation index"
-        ndvi.attrs["units"] = "dimensionless"
     ndvi_before.attrs["title"] = "NDVI before"
     ndvi_after.attrs["title"] = "NDVI after"
 
@@ -119,7 +111,7 @@ disaster:
 
 
 Tracking differences
---------------------
+++++++++++++++++++++
 
 An advantage of having our data in :class:`xarray.DataArray` format is that we
 can perform **coordinate-aware** calculations. This means that taking the
@@ -168,7 +160,7 @@ collapse in purple at the center.
 
 
 Estimating area
----------------
++++++++++++++++
 
 One things we can do with indices and their differences in time is calculated
 **area estimates**. If we know that the region of interest has index values
@@ -268,42 +260,35 @@ that many people will understand:
    `Silva Rotta et al. (2020) <https://doi.org/10.1016/j.jag.2020.102119>`__.
 
 
-Other indices
--------------
+Calculating your own indices
+----------------------------
 
-Calculating other indices will follow a very similar strategy to NDVI since
-most of them only involve arithmetic operations on different bands.
+Calculating other indices will is relatively straight forward since most of
+them only involve arithmetic operations on different bands.
 As an example, let's calculate and plot the
 `Modified Soil Adjusted Vegetation Index (MSAVI) <https://doi.org/10.1016/0034-4257(94)90134-1>`__
-for our two scenes:
+for our Manaus, Brazil, sample data:
 
 .. jupyter-execute::
 
-    import numpy as np
+    scene = xls.load_scene(xls.datasets.fetch_manaus())
 
-    # This time, use a loop and put them in a list to avoid repeated code
-    msavi_collection = []
-    for scene in [before, after]:
-        msavi = (
-            (
-                2 * scene.nir + 1 - np.sqrt(
-                    (2 * scene.nir + 1) * 2 - 8 * (scene.nir - scene.red)
-                )
-            ) / 2
-        )
-        msavi.name = "msavi"
-        msavi.attrs["long_name"] = "modified soil adjusted vegetation index"
-        msavi.attrs["units"] = "dimensionless"
-        msavi.attrs["title"] = scene.attrs["title"]
-        msavi_collection.append(msavi)
-
-    # Plotting is mostly the same
-    fig, axes = plt.subplots(2, 1, figsize=(10, 12), layout="tight")
-    for ax, msavi in zip(axes, msavi_collection):
-        msavi.plot(ax=ax, vmin=-0.5, vmax=0.5, cmap="RdBu_r")
-        ax.set_title(msavi.attrs["title"])
+    msavi = 0.5 * (
+        2 * scene.nir + 1
+        - np.sqrt((2 * scene.nir + 1) * 2 - 8 * (scene.nir - scene.red))
+    )
+    msavi.name = "msavi"
+    msavi.attrs["long_name"] = "modified soil adjusted vegetation index"
+
+    # Plot an RGB as well for comparison
+    rgb = xls.composite(scene, rescale_to=[0.02, 0.2])
+
+    fig, axes = plt.subplots(2, 1, figsize=(10, 9.75), layout="constrained")
+    rgb.plot.imshow(ax=axes[0])
+    msavi.plot(ax=axes[1], vmin=-0.5, vmax=0.5, cmap="RdBu_r")
+    axes[0].set_title("Manaus, Brazil")
+    for ax in axes:
         ax.set_aspect("equal")
     plt.show()
 
-**With this same logic, you could calculate NBR and dNBR, other variants of
-NDVI, NDSI, etc.**
+**With this same logic, you could calculate any other index.**

doc/plot-overlay.rst

@@ -28,10 +28,8 @@ version to get the thermal band as surface temperature:
 
     path_l1 = xls.datasets.fetch_momotombo(level=1)
     scene = xls.load_scene(path_l1)
-
     path_l2 = xls.datasets.fetch_momotombo(level=2)
     scene = scene.merge(xls.load_scene(path_l2, bands=["thermal"]))
-
     scene
 
 Now we can plot an RGB composite and thermal band separately to see that they
@@ -45,14 +43,11 @@ have to show:
     rgb = xls.adjust_l1_colors(rgb)
 
     # Plot the RGB and thermal separately
-    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 12))
-
-    rgb.plot.imshow(ax=ax1)
-    scene.thermal.plot.imshow(ax=ax2, cmap="magma")
-
-    ax1.set_aspect("equal")
-    ax2.set_aspect("equal")
-
+    fig, axes = plt.subplots(2, 1, figsize=(10, 12), layout="constrained")
+    for ax in axes:
+        ax.set_aspect("equal")
+    rgb.plot.imshow(ax=axes[0])
+    scene.thermal.plot.imshow(ax=axes[1], cmap="magma")
     plt.show()
 
 Notice that the lava flow is clearly visible as temperatures above 320 K in the
@@ -72,10 +67,8 @@ this is with the :func:`xarray.where` function:
     # If the condition is true, use the thermal values. If it's false, use nan
     lava = xr.where(scene.thermal >= 320, scene.thermal, np.nan, keep_attrs=True)
 
-    fig, ax = plt.subplots(1, 1, figsize=(10, 6))
-
+    fig, ax = plt.subplots(1, 1, figsize=(10, 6), layout="constrained")
     lava.plot.imshow(ax=ax, cmap="magma")
-
     ax.set_aspect("equal")
     plt.show()
 
@@ -94,12 +87,10 @@ plot that on top of the RGB composite and add a bit of transparency using the
 
 .. jupyter-execute::
 
-    fig, ax = plt.subplots(1, 1, figsize=(10, 6))
-
+    fig, ax = plt.subplots(1, 1, figsize=(10, 6), layout="constrained")
     # RGB goes first so it's at the bottom
     rgb.plot.imshow(ax=ax)
     lava.plot.imshow(ax=ax, cmap="magma", alpha=0.6)
-
     ax.set_aspect("equal")
     plt.show()
 

xlandsat/__init__.py

@@ -4,6 +4,7 @@
 from . import datasets
 from ._composite import composite
 from ._enhancement import adjust_l1_colors, equalize_histogram
+from ._indices import ndvi
 from ._interpolation import interpolate_missing
 from ._io import load_panchromatic, load_scene, save_scene
 from ._pansharpen import pansharpen

xlandsat/_indices.py

@@ -0,0 +1,41 @@
+# Copyright (c) 2022 The xlandsat developers.
+# Distributed under the terms of the MIT License.
+# SPDX-License-Identifier: MIT
+"""
+Calculate indices based on data in xarray.Datasets
+"""
+
+
+def ndvi(scene, red_band="red", nir_band="nir"):
+    r"""
+    Normalized Difference Vegetation Index
+
+    Calculate the NDVI for the given scene, defined as:
+
+    .. math::
+
+        NDVI = \dfrac{NIR - Red}{NIR + Red}
+
+    Parameters
+    ----------
+    scene : :class:`xarray.Dataset`
+        A Landsat scene, as read with :func:`xlandsat.load_scene`.
+    red_band : str
+        The name of the variable in ``scene`` that corresponds to the red band.
+    nir_band : str
+        The name of the variable in ``scene`` that corresponds to the NIR band.
+
+    Returns
+    -------
+    ndvi : :class:`xarray.DataArray`
+        The calculated NDVI, with the metadata attributes from the original
+        scene.
+    """
+    red = scene[red_band]
+    nir = scene[nir_band]
+    result = (nir - red) / (nir + red)
+    result.name = "ndvi"
+    attrs = {"long_name": "normalized difference vegetation index"}
+    attrs.update(scene.attrs)
+    result = result.assign_attrs(attrs)
+    return result