I've followed this tutorial to calculate the dNBR Index, using the given sample tif files which contain a burned area. My results were as shown on the tutorial. So far so good.
Then I tried to calculate the Burned Area Index (BAI) after making a few modifications on the above tutorial's code. The result is an empty plot.
This is my implementation:
(Note: The code downloads tif images from here)
from glob import glob
import matplotlib.pyplot as plt
import xarray as xr
import rioxarray as rxr
import earthpy.plot as ep
import earthpy as et
import os
def combine_tifs(tif_list):
"""A function that combines a list of tifs in the same CRS
and of the same extent into an xarray object
Parameters
----------
tif_list : list
A list of paths to the tif files that you wish to combine.
Returns
-------
An xarray object with all of the tif files in the listmerged into
a single object.
"""
out_xr = []
for i, tif_path in enumerate(tif_list):
out_xr.append(rxr.open_rasterio(tif_path, masked=True).squeeze())
out_xr[i]["band"] = i+1
return xr.concat(out_xr, dim="band")
path_to_downloaded_data = et.data.get_data('cold-springs-fire', replace=True)
all_landsat_bands_path = glob(os.path.join(path_to_downloaded_data,
"landsat_collect",
"LC080340322016072301T1-SC20180214145802",
"crop",
"*band[1-7]*.tif"))
all_landsat_bands_path.sort()
all_bands = combine_tifs(all_landsat_bands_path)
red_square = (0.1 - all_bands[3]) * (0.1 - all_bands[3])
nir_square = (0.06 - all_bands[4]) * (0.06 - all_bands[4])
burned_area_index = 1 / (red_square + nir_square)
fig1, ax = plt.subplots(figsize=(12, 12))
fig1 = plt.gcf()
ep.plot_bands(burned_area_index,
cmap='PiYG',
scale=False,
cbar=False, # Change this to TRUE to display the color bar
ax=ax,
vmin=-1, vmax=1,
title="BAI")
plt.show()
The values of the burned_area_index array are very low, in the magnitude of 10^-7.
My question is, if the formula is correct, what might be the problem that results in an empty plot?
