in geospatial/visualization/imagery_visualizer.py [0:0]
def get_landsat8_ndvi(tile: Union[rasterio.DatasetReader, np.ndarray],
window: Union[rasterio.windows.Window, Tuple] = None) -> np.ndarray:
"""Computes the NDVI (Normalized Difference Vegetation Index) over a tile or a section
on the tile specified by the window, for Landsat 8 tiles.
NDVI values are between -1.0 and 1.0 (hence "normalized"), mostly representing greenness, where any
negative values are mainly generated from clouds, water, and snow, and values
near zero are mainly generated from rock and bare soil. Very low values (0.1 and below)
of NDVI correspond to barren areas of rock, sand, or snow. Moderate values (0.2 to 0.3)
represent shrub and grassland, while high values (0.6 to 0.8) indicate temperate
and tropical rainforests. Source: https://desktop.arcgis.com/
For surface reflectance products, there are some out of range values in the red and near infrared bands
for areas around water/clouds. We should get rid of these before NDVI calculation,
otherwise the NDVI will also be out of range.
Source: https://www.researchgate.net/post/Why_Landsat_8_NDVI_Values_are_out_of_Range_Not_in_between-1_to_1
If we decide not to get rid of invalid values in the red and NIR bands and see some out-of-range NDVI values,
we can check to see if these out-of-range NDVI values are a small minority in all pixels in the dataset.
Args:
tile:
- a rasterio.io.DatasetReader object returned by rasterio.open(), or
- a numpy array of dims (height, width, bands)
window: a tuple of four (col_off x, row_off y, width delta_x, height delta_y)
to specify the section of the tile to compute NDVI over, or a rasterio Window object
Returns:
2D numpy array of dtype float32 of the NDVI values at each pixel, of dims (height, width)
Pixel value is set to 0 if the sum of the red and NIR value there is 0 (empty).
"""
if isinstance(tile, rasterio.io.DatasetReader):
if window and isinstance(window, Tuple):
window = rasterio.windows.Window(window[0], window[1], window[2], window[3])
band_red = tile.read(4, window=window, boundless=True, fill_value=0).squeeze()
band_nir = tile.read(5, window=window, boundless=True, fill_value=0).squeeze()
elif isinstance(tile, np.ndarray):
if window and isinstance(window, rasterio.windows.Window):
window = [window.col_off, window.row_off, window.width, window.height]
# get the red and NIR bands
tile_bands = []
for numpy_band_index in (3, 4):
if window:
b = tile[
window[1]: window[1] + window[3],
window[0]: window[0] + window[2],
numpy_band_index
]
else:
b = tile[:, :, numpy_band_index]
tile_bands.append(b)
band_red = tile_bands[0]
band_nir = tile_bands[1]
sum_red_nir = band_nir + band_red
# sum of the NIR and red bands being zero is most likely because this section is empty
# this workaround means that the final NDVI at such pixels are 0.
sum_red_nir[sum_red_nir == 0.0] = 1
ndvi = (band_nir - band_red) / sum_red_nir
return ndvi