CODEWITHSUNDEEP
pythongis
user13644866
user13644866

Reputation:

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.

#Read the image file
image_file = "my_tif_image.tif"


with rasterio.open(image_file) as src:
    band_red = src.read(3)

with rasterio.open(image_file) as src:
    band_nir = src.read(4)
   

# 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?

Upvotes: 3

Views: 1072

Answers (1)

gx0113
gx0113

Reputation: 36

To use the Burned Area Index, the data must be calibrated to reflectance. There are other vegetation indices:

You can apply the different indices in python or R (which cointains pre-establish packages for the manipulation). Other options are QGIS, PCI Geomatica (license needed) and ArcMAP (license Needed).

Upvotes: 2

Related Questions