Compute a spectral index from S2 Time Series

This example shows how to compute an index (here LChloC) from a S2 with 10 bands. The raster is order date per date (blue,green,red…date 1 then blue,green,red… date 2…)

Import libraries

import numpy as np
from museopheno import sensors,datasets
from museotoolbox.processing import RasterMath

Import dataset

raster,dates = datasets.Sentinel2_3a_2018(return_dates=True)

Define an instance of sensors.Sentinel2()

S2 = sensors.Sentinel2(n_bands=10)

# check default band_order
print('Default band order for 10 bands is : '+', '.join(S2.bands_order)+'.')

# List of available index :
S2.available_indices.keys()

Out:

Default band order for 10 bands is : 2, 3, 4, 8, 5, 6, 7, 8A, 11, 12.

dict_keys(['B2', 'B3', 'B4', 'B8', 'B5', 'B6', 'B7', 'B8A', 'B11', 'B12', 'ACORVI', 'ACORVInarrow', 'SAVI', 'EVI', 'EVI2', 'PSRI', 'ARI', 'ARI2', 'MARI', 'CHLRE', 'MCARI', 'MSI', 'MSIB12', 'NDrededgeSWIR', 'SIPI2', 'NDWI', 'LCaroC', 'LChloC', 'LAnthoC', 'Chlogreen', 'NDRESWIR', 'NDVI', 'NDVInarrow', 'NDVIre', 'RededgePeakArea', 'Rratio', 'MTCI', 'S2REP', 'IRECI', 'NBR', 'REP', 'WDRVI_B6B7'])

Write metadata in each band (date + band name)

S2.set_description_metadata(raster,dates)

Generate a raster with LChloC index

# show expression and condition of LChloC index
print(S2.get_index_expression('LChloC'))

# generate raster
S2.generate_raster(input_raster=raster,output_raster='/tmp/S2.tif',expression=S2.get_index_expression('LChloC'),dtype=np.float32)

Out:

{'expression': 'B7 / B5', 'condition': 'B5 != 0'}
Total number of blocks : 1
Using datatype from numpy table : float32.
Detected 7 bands for function expression_manager.
Batch processing (1 blocks using 11Mo of ram)

Computing index [########################################]100%

Plot image

rM = RasterMath(raster)
X=rM.get_random_block()
NDVI = S2.generate_index(X,S2.get_index_expression('LChloC'),dtype=np.float32)

from matplotlib import pyplot as plt
from datetime import datetime
dateToDatetime = [datetime.strptime(str(date),'%Y%m%d') for date in dates]
plt.plot_date(dateToDatetime,NDVI[:10,:].T,'-o')
plt.ylabel('Leaf Chlorophyll Content')

Out:

Total number of blocks : 1

Text(0, 0.5, 'Leaf Chlorophyll Content')