I am attempting to calculate a vegetative index from a mask layer for vegetation in a GeoTIFF image. The vegetative index is calculated simply as (r - b)/r on all pixels such that r !=0 and that exceed a certain threshold to be included in the mask layer.
If my math is right, the range of this VI is [-1, 1], however when the vegetative index geotiff is viewed in QGIS, the values for this VI range from 0 to 65,534:
When inspecting individual pixels with the "Value Tool" plugin, the individual RGB bands correspond to the manually calculated VI, except for pixels with very high values, such as this example:
Does anyone have an idea as to what the issue could be here?
# Coloration index
def ci_ind(output_file, r, g, b, threshold, affine, projection):
np.seterr(divide='ignore', invalid='ignore')
mask = np.where(r != 0, 2 * (g / (r + g + b)) - r / (r + g + b) - b / (r + g + b) >= threshold, np.nan)
check = mask == 1
vi = np.where(check, (r - b) / r, np.nan)
geotiff_exp(vi, output_file, affine, projection)
# Geotiff export
def geotiff_exp(vi_array, output_file, affine, projection):
driver = gdal.GetDriverByName('GTiff')
driver.Register()
outds = driver.Create(output_file, xsize=vi_array.shape[1],
ysize=vi_array.shape[0], bands=1,
eType=gdal.GDT_Float32)
outds.SetGeoTransform(affine)
outds.SetProjection(projection)
outband = outds.GetRasterBand(1)
outband.WriteArray(vi_array)
outband.SetNoDataValue(0)
outband.FlushCache()
outband = None
outds = None
fp = r"TEST 10-20-21 Area H NNN Plant Rows.tif"
ds = gdal.Open(fp)
gt = ds.GetGeoTransform()
proj = ds.GetProjection()
exg_thresh = 0.05
r = ds.GetRasterBand(1).ReadAsArray().astype(np.uint16)
g = ds.GetRasterBand(2).ReadAsArray().astype(np.uint16)
b = ds.GetRasterBand(3).ReadAsArray().astype(np.uint16)
ci_ind('ci_test.tif', r, g, b, exg_thresh, gt, proj)
