Welcome to sentipy’s documentation page

Getting started

Installation

sentipy is hosted on pypy and can therefore be pip installed:

$ pip install sentipy

Usage

The intended responsibility of the sentipy project is to process pre-loaded Sentinel imagery into useful derivative indices such as those included in the S2 Toolbox (Fcover, Fapar, Canopy water). This is likely to expand into other spectral indices as the project grows.

We therefore do not involve ourselves in the loading of bands / pixels from imagery; our contract assumes that you will be able to provide pixel values for each band (as required) in the form of numpy (masked) arrays, and we take care of the processing logic.

Our preferred method for loading Sentinel imagery is to use the rasterio package.

from sentipy import s2_toolbox
import rasterio

with rasterio.open('example.tif') as dataset:
    # read your bands in & stack them into a single array here
    input_arr = np.stack(
        [
            dataset.read(1),
            dataset.read(2),
            dataset.read(3),
            # ...
            # and you'll need the image meta-data here too: view zenith, sun zenith, rel azimuth
        ]
        , axis=0
    )

# Initialise whichever calculator you wish to use, eg the Fapar calculator:
calculator = s2_toolbox.Fapar()
# All you have to do is call the calculator's `.run()` method on the input array and you'll get an output array of
# matching size (but not shape) with your calculated value(s)
fapar = calculator.run(input_arr)

And that’s really all there is to it!

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

API Reference

sentipy.settings

Settings and configuration that apply across the whole project

Sentinel 2 bands

https://sentinel.esa.int/web/sentinel/user-guides/sentinel-2-msi/resolutions/radiometric

Band label	Band code	Band description	Central wavelength	Resolution (m)
Band 1	“B01”	Coastal aerosol	0.443	60
Band 2	“B02”	Blue	0.49	10
Band 3	“B03”	Green	0.56	10
Band 4	“B04”	Red	0.665	10
Band 5	“B05”	Vegetation Red Edge	0.705	20
Band 6	“B06”	Vegetation Red Edge	0.74	20
Band 7	“B07”	Vegetation Red Edge	0.783	20
Band 8	“B08”	NIR	0.842	10
Band 8A	“B8a”	Vegetation Red Edge	0.865	20
Band 9	“B09”	Water vapour	0.945	60
Band 10	“B10”	SWIR - Cirrus	1.375	60
Band 11	“B11”	SWIR	1.61	20
Band 12	“B12”	SWIR	2.19	20

> The Band code is used to refer to Sentinel-2 bands throughout the project.

sentipy.settings.DEFAULT_BAND_SEQUENCE = ('B03', 'B04', 'B05', 'B06', 'B07', 'B8a', 'B11', 'B12', 'COS_VIEW_ZENITH', 'COS_SUN_ZENITH', 'COS_REL_AZIMUTH')

Default sequence of Sentinel-2 bands & meta-data assumed throughout project

sentipy.s2_toolbox

Python implementation of some Sentinel-2 toolbox processing tools.

Home for the objects & logic that deal with the calculation of FAPAR

class sentipy.s2_toolbox.CanopyWater

Bases: sentipy.s2_toolbox.S2BiophysicalCalculator

Calculator for the Canopy Water biophysical parameter.

VALIDATION_RANGES = {'B03': {'max': 0.263, 'min': 0.0}, 'B04': {'max': 0.3, 'min': 0.0}, 'B05': {'max': 0.315, 'min': 0.0}, 'B06': {'max': 0.619, 'min': 0.0}, 'B07': {'max': 0.764, 'min': 0.004}, 'B11': {'max': 0.503, 'min': 0.006}, 'B12': {'max': 0.503, 'min': 0.0}, 'B8a': {'max': 0.792, 'min': 0.017}, 'Cos(rel. azimuth)': {'max': 1.0, 'min': -1.0}, 'Cos(sun zenith)': {'max': 1.0, 'min': 0.0}, 'Cos(view zenith)': {'max': 1.0, 'min': 0.0}}

class sentipy.s2_toolbox.Fapar

Bases: sentipy.s2_toolbox.S2BiophysicalCalculator

Calculator for the Fapar biophysical parameter; the fraction of photosynthetically-active radiation absorbed (by vegetation).

VALIDATION_RANGES = {'B03': {'max': 0.263, 'min': 0.0}, 'B04': {'max': 0.3, 'min': 0.0}, 'B05': {'max': 0.315, 'min': 0.0}, 'B06': {'max': 0.619, 'min': 0.0}, 'B07': {'max': 0.764, 'min': 0.004}, 'B11': {'max': 0.503, 'min': 0.006}, 'B12': {'max': 0.503, 'min': 0.0}, 'B8a': {'max': 0.792, 'min': 0.017}, 'Cos(rel. azimuth)': {'max': 1.0, 'min': -1.0}, 'Cos(sun zenith)': {'max': 1.0, 'min': 0.0}, 'Cos(view zenith)': {'max': 1.0, 'min': 0.0}}

class sentipy.s2_toolbox.Fcover

Bases: sentipy.s2_toolbox.S2BiophysicalCalculator

Calculator for the Fcover biophysical parameter; the fraction of ground covered by vegetation.

class sentipy.s2_toolbox.S2BiophysicalCalculator

Bases: abc.ABC

run(input_arr: <sphinx.ext.autodoc.importer._MockObject object at 0x7fa369017950>, band_sequence: List[str] = ('B03', 'B04', 'B05', 'B06', 'B07', 'B8a', 'B11', 'B12', 'COS_VIEW_ZENITH', 'COS_SUN_ZENITH', 'COS_REL_AZIMUTH'), validate: bool = True) → Union[float, <sphinx.ext.autodoc.importer._MockObject object at 0x7fa368eeacd0>]

Run the calculator on an input array

By default, the calculator expects only the following bands to be passed in the sequence below:

• B03
• B04
• B05
• B06
• B07
• B8a
• B11
• B12
• COS_VIEW_ZENITH
• COS_SUN_ZENITH
• COS_REL_AZIMUTH

If band values are to be passed in a different set or sequence, the band_sequence parameter must be passed with band names (matching those above) for each element in the input array. eg. [“extra_band_1”, “COS_SUN_ZENITH”, “B03”, “B04”, …, “COS_REL_AZIMUTH”, “extra_band_2”]

Parameters:

• input_arr – Input values for the calculator to use
• band_sequence – Names of bands included in the input array (names must match those used above for the required bands)
• validate – Flag for whether or not to apply validation ranges to the inputs

Returns:

Scalar estimate of the biophysical property

sentipy.agri_indices

A simple python API for processing raster arrays into common indices used in the agricultural domain.

There are a wide range of agriculturally-relevant indices to be processed from (multi-spectral) optical satellite imagery in general.

This module provides implementations of a number of indices, largely inspired by the Index Database and the SentinelHub custom scripts repository. They are implemented principally with Sentinel-2 in mind, but can be applied to imagery from other constellations where equivalent bands are available (it is left to the user to decide where this is the case).

Available implementations:

• NDVI
• NDWI
• OSAVI

sentipy.agri_indices.ndmi(img_arr, band_sequence: tuple = ('B08', 'B11'))

Normalized Difference Moisture Index

NDMI = (NIR - SWIR)/(NIR + SWIR)

= (B08 - B11)/(B08 + B11)

sentipy.agri_indices.ndvi(img_arr, band_sequence: tuple = ('B04', 'B08'))

Normalized Difference Vegetation Index

NDVI = (NIR - RED)/(NIR + RED)

= (B08 - B04)/ (B08 + B04)

sentipy.agri_indices.ndwi(img_arr, band_sequence: tuple = ('B03', 'B08'))

Normalized Difference Water Index

NDWI = -(GREEN - NIR)/(GREEN + NIR)

= -(B03 - B08)/(B03 + B08)

sentipy.agri_indices.osavi(img_arr, band_sequence: tuple = ('B04', 'B08'))

Optimized Soil-Adjusted Vegetation Index Index DB link

OSAVI = (1+0.16)*(NIR - RED)/(NIR + RED + 0.16)

= (1.16)*(B08 - B04)/(B08 + B04 + 0.16)

Indices and tables

• Index
• Module Index
• Search Page