How to calculate Multispectral Vegetation Indices using Python

I am trying to calculate Multispectral Vegetation Indices from multi spectral photos, using Python.

So far, I've managed to calculate the NDVI, which is easy to implement as there are many examples online.

image_file = "my_tif_image.tif"

with rasterio.open(image_file) as src:

with rasterio.open(image_file) as src:

# Allow division by zero
numpy.seterr(divide='ignore', invalid='ignore')

band_nir = numpy.array(band_nir)
band_red = numpy.array(band_red)

NDVI = (band_nir.astype(float) - band_red.astype(float)) / (band_nir + band_red)

fig1, ax = plt.subplots(figsize=(12, 12))
fig1 = plt.gcf()
ep.plot_bands(NDVI,
cmap='PiYG',
scale=False,
cbar = False,  #Change this to TRUE to display the color bar
ax=ax,
vmin=-1, vmax=1,
title="NDVI")
plt.show()
fig1.savefig('NDVI.png',bbox_inches='tight')

However, apart from the above index, I couldn't manage to find anywhere examples on how to calculate other indexes (at least the ones I am interested in).

For example, instead of the NDVI index, I tried the formula of Burned Area Index:

BIR = 1.0 / (((0.1 + band_red) ** 2) + ((0.06 + band_nir) ** 2))

but I get an empty plot and the calculated array numbers are very small (in the numbers of 10^-9).

My question is, why is it so hard to find examples that calculate and plot indices other than the NDVI? And lastly, if I can't do it with python, is there any other programmatic way to do it?

I've stumbled across the Arcpy, but I wasn't able to import that module. Do I have to have the ArcGIS software in order to use that library?

• indexdatabase.de/db/i.php
– BERA
Nov 20, 2020 at 12:25
• Yes, you do need to have ArcGIS to use the arcpy module Nov 20, 2020 at 12:58
• You can easily download QGIS and install it. You can look at the multiple plugins available ie plugins.qgis.org/plugins/indices. You can also run your python from there. Nov 21, 2020 at 16:47