Examples

Processing a Sentinel-2 raster image into actionable info.

Processing a Sentinel-2 raster image into agri-relevant indices¶

In this example, we'll show how easy it is to use the sentipy package to convert Sentinel-2 imagery into useful indices such as NDVI and OSAVI, and compute indicators for bio-physical properties of the fields (such as fcover and fapar).

Overview¶

We'll assume that you have some imagery that you want to work with, and use the rasterio package to read an image (in GeoTIFF format) in and convert the values into NDVI values.

We import the sentipy agri_indices module to compute NDVI on the numpy image array that rasterio gives us.

Imports¶

In [1]:

import rasterio as rio
import numpy as np
import numpy.ma as ma
from rasterio.plot import show, show_hist
import matplotlib.pyplot as plt

Reading in the data using rasterio¶

We have a test fixture in the respository that we use to verify the output of our index calculations. We'll work with that in this example.

We begin by creating a rasterio Dataset object and then reading the individual bands that we're interested in into numpy arrays. You'll see that we work with masked arrays as a convention, but it isn't actually necessary for this imagery.

In [3]:

with rio.open('../tests/fixtures/s2_example_input.tif') as dataset:
    band2 = ma.masked_array(dataset.read(2))
    band3 = ma.masked_array(dataset.read(3))
    band4 = ma.masked_array(dataset.read(4))
    band5 = ma.masked_array(dataset.read(5))
    band6 = ma.masked_array(dataset.read(6))
    band7 = ma.masked_array(dataset.read(7))
    band8 = ma.masked_array(dataset.read(8))
    band8a = ma.masked_array(dataset.read(9))
    band11 = ma.masked_array(dataset.read(11))
    band12 = ma.masked_array(dataset.read(12))

Understanding the data¶

Visualising the imagery is an important first step and should help to get a sense for the output we'd expect.

For example, NDVI (Normalised Difference Vegetation Index) - as the name suggests - should correlate strongly with healthy, green vegetation.

Looking at an RGB rendering of the imagery below, we can see that the imagery shows a region with many fields; some of which are bare, some of which are densely vegetated, and some of which lie in between...

In [4]:

rgb = np.stack([band4, band3, band2], axis=2) / 0.3
plt.imshow(rgb)

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

Out[4]:

<matplotlib.image.AxesImage at 0x7ff1334bc220>

Multi-spectral bands¶

Chlorophyll has high reflectance in the near-infrared band (NIR). Many vegetative indices therefore use this band in their calculations to remotely assess the state of crops under monitoring. It can be helpful to consider how the intensity of each band varies across the scene to better intuit how these vegetative indices work - especially for the bands that we can't actually see.

It's interesting to note that the red band is typically higher (more red) on the field parcels that are not heavily vegetated and, conversely, the NIR band is most intense where the vegetation is most dense.

In [5]:

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15,10))
ax1 = show(band4, ax=ax1, cmap="Reds", title="red")
ax2 = show(band3, ax=ax2, cmap="Greens", title="green")
ax3 = show(band2, ax=ax3, cmap="Blues", title="blue")
ax4 = show(band8, ax=ax4, cmap="Purples", title="NIR")
fig.show()

<ipython-input-5-f58baa5666d3>:6: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()

Calculate Vegetative indices & canopy moisture indicators¶

We now use sentipy to calculate a few vegetative indices:

For each index, we:

• stack the individual bands that we need (eg.red & NIR) in a numpy array
• pass the array to the sentipy method, eg. agri_indices.ndvi()
• visualise the result both a spatial map in a divergent colourscheme, and a histogram of the index values for each pixel

In [6]:

import matplotlib
font = {'weight' : 'bold',
        'size'   : 22}

matplotlib.rc('font', **font)

In [7]:

img_arr = np.stack([band4, band8], axis=2)
ndvi = agri_indices.ndvi(img_arr)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(30,10))
ax1.set_axis_off()
_ = show(ndvi, ax=ax1, cmap="RdYlGn")
_ = show_hist(ndvi, ax=ax2, bins=30, stacked=False, alpha=0.6, histtype='stepfilled', label='frequency')
fig.show()

<ipython-input-7-aefde010eba6>:8: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()

As you can hopefully see, the NDVI is highest (shown in deeper green) on the fields where vegetation was most dense in the input image & it is lowest (shown in light shades) where the fields were bare.

And hopefully you could also see that calculating the NDVI with sentipy was easy!

Looking at the histogram for NDVI values, we can see that a lot of pixels are bunched up in the RHS of the distribution. This is because NDVI saturates (ie. hits the values near its max of 1) faster than some other indices.

Next, we'll compare NDVI to another popular VI; OSAVI.

In [17]:

img_arr = np.stack([band4, band8], axis=2)
osavi = agri_indices.osavi(img_arr)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(30,10))
ax1.set_axis_off()
_ = show(osavi, ax=ax1, cmap="RdYlGn")
_ = show_hist(osavi, ax=ax2, bins=30, stacked=False, alpha=0.6, histtype='stepfilled', label='frequency')
fig.show()

<ipython-input-17-0f28efd7dc79>:8: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()

As we can see in the histogram, the values are less bunched at the extreme RHS - a sign that OSAVI does not 'suffer' from the same issue of early saturation.

Rather than vegetation, we can also use indices to indicate moisture conditions in the field (although total canopy will correlate strongly with total moisture). In the example below, we use another normalised-difference index called NDMI (darker blue indicates more moisture).

In [21]:

img_arr = np.stack([band8, band11], axis=2)
ndmi = np.clip(agri_indices.ndmi(img_arr), a_min=0, a_max=1.)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(30,10))
ax1.set_axis_off()
_ = show(ndmi, ax=ax1, cmap="RdYlBu")
_ = show_hist(ndmi, ax=ax2, bins=30, stacked=False, alpha=0.6, histtype='stepfilled', label='frequency')
fig.show()

<ipython-input-21-20066eab9022>:8: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()

Using insights from the sentinel-2 toolbox¶

In addition to these commonly-used indices, sentipy can also be used to reproduce some of the processing outputs available in the S2 toolbox, including Fapar, Fcover and Canopy water.

For the imagery here, we pick values within the accepted range for the metadata values (view zenith, sun zenith and relative azimuth), but users will need to retrieve these values with their imagery.

In [22]:

from sentipy.s2_toolbox import Fapar, Fcover, CanopyWater

for each toolbox output, we need all of the bands below:

In [10]:

cos_view_zenith = np.full(band4.shape, 0.5)
cos_sun_zenith = np.full(band4.shape, 0.5)
cos_rel_azimuth = np.full(band4.shape, 0.)

toolbox_input = np.stack(
    [
        band3,
        band4,
        band5,
        band6,
        band7,
        band8a,
        band11,
        band12,
        cos_view_zenith,
        cos_sun_zenith,
        cos_rel_azimuth
    ],
    axis=2
)

In [23]:

fapar_processor = Fapar()
fapar = np.clip(np.apply_along_axis(fapar_processor.run, 2, toolbox_input, validate=False), a_min=0., a_max=1.)

In [24]:

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(30,10))
ax1.set_axis_off()
_ = show(fapar, ax=ax1, cmap="RdYlGn")
_ = show_hist(fapar, ax=ax2, bins=30, stacked=False, alpha=0.6, histtype='stepfilled', label='frequency')
fig.show()

<ipython-input-24-37be08475231>:5: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()

In [13]:

fcover_processor = Fcover()
fcover = np.apply_along_axis(fcover_processor.run, 2, toolbox_input, validate=False)

In [14]:

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(30,10))
ax1.set_axis_off()
_ = show(fcover, ax=ax1, cmap="RdYlGn")
_ = show_hist(fcover, ax=ax2, bins=30, stacked=False, alpha=0.6, histtype='stepfilled', label='frequency')
fig.show()

<ipython-input-14-2b010e2d9c7f>:5: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()

In [25]:

cwater_processor = CanopyWater()
cwater = np.clip(np.apply_along_axis(cwater_processor.run, 2, toolbox_input, validate=False), a_min=0., a_max=1.)

In [26]:

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(30,10))
ax1.set_axis_off()
_ = show(cwater, ax=ax1, cmap="RdYlBu")
_ = show_hist(cwater, ax=ax2, bins=30, stacked=False, alpha=0.6, histtype='stepfilled', label='frequency')
fig.show()

<ipython-input-26-31ab7c124a89>:5: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()