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.
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¶
import sys
sys.path.append('../../sentipy')
import rasterio as rio
import numpy as np
import numpy.ma as ma
from rasterio.plot import show
import matplotlib.pyplot as plt
from sentipy import agri_indices
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.
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 should help to get a sense for the output we'd expect.
NDVI (Normalised Difference Vegetation Index) - as the name suggests - should correlate strongly with healthy, green vegetation.
Look 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...
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).
<matplotlib.image.AxesImage at 0x7f825efc2b80>
Multi-spectral bands¶
NDVI is computed using the red and the near-infrared bands. It can be helpful to consider how the intensity of each band varies across the scene to understand 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.
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-15-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¶
Finally, we use sentipy to calculate ndvi for us:
- We grab the bands that we need (red & NIR), and stack them in a numpy array
- Pass the array to
agri_indices.ndvi()
- Visualise the result (using a rasterio utility function)
img_arr = np.stack([band4, band8], axis=2)
ndvi = agri_indices.ndvi(img_arr)
_ = show(ndvi, cmap="RdYlGn")
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 from the imagery was easy!
img_arr = np.stack([band8, band11], axis=2)
ndmi = agri_indices.ndmi(img_arr)
_ = show(ndvi, cmap="RdYlBu")
Using insights from the sentinel-2 toolbox¶
from sentipy.s2_toolbox import Fapar, Fcover, CanopyWater
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
)
fapar_processor = Fapar()
fapar = np.apply_along_axis(fapar_processor.run, 2, toolbox_input, validate=False)
_ = show(fapar, cmap="RdYlGn")
fcover_processor = Fcover()
fcover = np.apply_along_axis(fcover_processor.run, 2, toolbox_input, validate=False)
_ = show(fcover, cmap="RdYlGn")
cwater_processor = CanopyWater()
cwater = np.apply_along_axis(cwater_processor.run, 2, toolbox_input, validate=False)
_ = show(cwater, cmap="RdYlBu")