Generate a NDVI time series

This example shows how to produce a NDVI time series from a Sentinel-2 time series.

Import libraries

from museopheno import sensors,datasets
import numpy as np

Import dataset

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

Create an instance of sensors.Sentinel2()

S2 = sensors.Sentinel2(n_bands=10)

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

Out:

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

List of available index :

print(S2.available_indices.keys())

Out:

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)

This is useful to show date in band number in Qgis

S2.set_description_metadata(raster,dates)

Produce NDVI time series from and to a raster

S2.generate_raster(input_raster=raster,output_raster='/tmp/index.tif',expression=S2.get_index_expression('NDVI'),dtype=np.float32)

Out:

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 NDVI index

from museotoolbox.processing import RasterMath
from matplotlib import pyplot as plt

rM = RasterMath('/tmp/index.tif')
NDVI = rM.get_random_block() #randomly select a block
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('NDVI')

Out:

Total number of blocks : 1

Text(0, 0.5, 'NDVI')