Produce a new raster with index stacked in new bands

This example shows how to add to a raster spectral index. Here we add LChloC and ACORVI (a modified NDVI).

Import libraries

import numpy as np
from museopheno import sensors,datasets

# to add custom  creation of new raster, import rasterMath
from museotoolbox.processing import RasterMath

Import dataset

raster,dates = datasets.Sentinel2_3a_2018(return_dates=True)
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 :

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'])

Define a custom function

def createSITSplusIndex(X):
    X1 = S2.generate_index(X,S2.get_index_expression('LChloC'),multiply_by=100,dtype=np.int16)
    X2 = S2.generate_index(X,S2.get_index_expression('ACORVI'),multiply_by=100,dtype=np.int16)

    return np.hstack((X,X1,X2)).astype(np.int16)

Use rasterMath to read and write block per block the raster according to a function

rM = RasterMath(raster)

X = rM.get_random_block()
print('Block has {} pixels and {} bands'.format(X.shape[0],X.shape[1]))

Out:

Total number of blocks : 1
Block has 55842 pixels and 70 bands

Now we can try our function

XwithIndex = createSITSplusIndex(X)
print('Raster+index produced has {} pixels and {} bands'.format(XwithIndex.shape[0],XwithIndex.shape[1]))

Out:

Raster+index produced has 55842 pixels and 84 bands

Now we add our function as the test was a success

rM.add_function(createSITSplusIndex,'/tmp/SITSwithIndex.tif')

Out:

Using datatype from numpy table : int16.
Detected 84 bands for function createSITSplusIndex.

Produce raster

rM.run()

##################
# Plot image
from matplotlib import pyplot as plt
rM = RasterMath('/tmp/SITSwithIndex.tif')
X=rM.get_random_block() #randomly select a block
plt.plot(X[:20,:].T)
addIndextoRaster

Out:

Batch processing (1 blocks using 20Mo of ram)

rasterMath... [########################################]100%
Total number of blocks : 1

[<matplotlib.lines.Line2D object at 0x7f086040c1d0>, <matplotlib.lines.Line2D object at 0x7f085bddcdd0>, <matplotlib.lines.Line2D object at 0x7f086040c510>, <matplotlib.lines.Line2D object at 0x7f086040c990>, <matplotlib.lines.Line2D object at 0x7f086040c550>, <matplotlib.lines.Line2D object at 0x7f086040ca10>, <matplotlib.lines.Line2D object at 0x7f086040c750>, <matplotlib.lines.Line2D object at 0x7f086040c910>, <matplotlib.lines.Line2D object at 0x7f086040c350>, <matplotlib.lines.Line2D object at 0x7f086040c190>, <matplotlib.lines.Line2D object at 0x7f085bdc4250>, <matplotlib.lines.Line2D object at 0x7f0860066490>, <matplotlib.lines.Line2D object at 0x7f0860066710>, <matplotlib.lines.Line2D object at 0x7f0860066250>, <matplotlib.lines.Line2D object at 0x7f0860066110>, <matplotlib.lines.Line2D object at 0x7f0860066850>, <matplotlib.lines.Line2D object at 0x7f0860066190>, <matplotlib.lines.Line2D object at 0x7f0860066150>, <matplotlib.lines.Line2D object at 0x7f0860066210>, <matplotlib.lines.Line2D object at 0x7f0860066810>]

Total running time of the script: ( 0 minutes 6.887 seconds)

Gallery generated by Sphinx-Gallery