From 97cfbbcb5d66f202b140da13e3d52d1fccc52e32 Mon Sep 17 00:00:00 2001 From: Varitius24 <149641921+Varitius24@users.noreply.github.com> Date: Sat, 24 Feb 2024 15:40:03 +0000 Subject: [PATCH] Update and rename __init__.py to health_monitor.py --- src/__init__.py | 23 ----- src/health_monitor.py | 193 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 193 insertions(+), 23 deletions(-) delete mode 100644 src/__init__.py create mode 100644 src/health_monitor.py diff --git a/src/__init__.py b/src/__init__.py deleted file mode 100644 index 4a8677e..0000000 --- a/src/__init__.py +++ /dev/null @@ -1,23 +0,0 @@ -import os - -import toml - - - -def get_version_from_pyproject(): - try: - # Determine the path to pyproject.toml relative to this file - base_dir = os.path.dirname(os.path.abspath(__file__)) - pyproject_path = os.path.join(base_dir, "..", "pyproject.toml") - - # Load the pyproject.toml file - pyproject_data = toml.load(pyproject_path) - - # Extract the version from the project section - version = pyproject_data["project"]["version"] - except FileNotFoundError: - version = "unknown" - return version - - -__version__ = get_version_from_pyproject() diff --git a/src/health_monitor.py b/src/health_monitor.py new file mode 100644 index 0000000..d14548b --- /dev/null +++ b/src/health_monitor.py @@ -0,0 +1,193 @@ +import numpy as np +from osgeo import gdal +import matplotlib.pyplot as plt +from settings import * + +def ndvi(nir: np.ndarray, red: np.ndarray)-> np.ndarray: + top = nir - red + bot = nir + red + non_zero_mask = bot != 0 + valid_mask = ~np.isnan(top) & ~np.isnan(bot) + final_mask = valid_mask & non_zero_mask + result = np.divide(top, bot, where=final_mask, out=np.full_like(top, np.nan)) + return result + +def ndre(nir: np.ndarray, re: np.ndarray)-> np.ndarray: + top = nir-re + bot = nir+re + non_zero_mask = bot != 0 + valid_mask = ~np.isnan(top) & ~np.isnan(bot) + final_mask = valid_mask & non_zero_mask + result = np.divide(top, bot, where=final_mask, out=np.full_like(top, np.nan)) + return result + +def sr(nir: np.ndarray, red: np.ndarray) -> np.ndarray: + non_zero_mask = red != 0 + valid_mask = ~np.isnan(nir) & ~np.isnan(red) + final_mask = valid_mask & non_zero_mask + result = np.divide(nir, red, where=final_mask, out=np.full_like(nir, np.nan)) + return result + +def srre(nir: np.ndarray, red_edge: np.ndarray) -> np.ndarray: + non_zero_mask = red_edge != 0 + valid_mask = ~np.isnan(nir) & ~np.isnan(red_edge) + final_mask = valid_mask & non_zero_mask + result = np.divide(nir, red_edge, where=final_mask, out=np.full_like(nir, np.nan)) + return result + +def evi(nir: np.ndarray, red: np.ndarray, blue: np.ndarray, g =2.5, c1 = 6, c2 = 7.5, l = 1, scale = 0.0001) -> np.ndarray: + + #for some reason you have to scale it down? See below: + #https://community.esri.com/t5/arcgis-spatial-analyst-questions/formula-for-enhanced-vegetation-index-evi-in/td-p/32453 + top = (nir*scale-red*scale) *g + bot = nir*scale+(c1*red*scale)-(c2*blue*scale)+l + non_zero_mask = bot != 0 + valid_mask = ~np.isnan(top) & ~np.isnan(bot) + final_mask = valid_mask & non_zero_mask + result = np.divide(top, bot, where=final_mask, out=np.full_like(top, np.nan)) + return result + +def msavi(nir: np.ndarray, red: np.ndarray) -> np.ndarray: + nir2_1 = 2*nir + 1 + nir_red = nir-red + sqr_root = np.sqrt((nir2_1)**2 - 8*nir_red) + top = nir2_1 - sqr_root + final_mask = ~np.isnan(top) + result = np.divide(top, 2, where=final_mask, out=np.full_like(top, np.nan)) + return result + + +def scaleMinMax(x : np.ndarray) -> np.ndarray: + return((x - np.nanmin(x))/(np.nanmax(x) - np.nanmin(x))) + +def get_bands(min : int ,max: int(),image): + band = 0 + for band in range(min,max+1): + band += np.array(image.GetRasterBand(band).ReadAsArray()) + return scaleMinMax(band) + +def make_and_save_graph(title, data): + plt.figure() + plt.title(title) + plt.imshow(data,cmap = "RdYlGn") + plt.colorbar(label='Intensity') + plt.savefig(title +'.png') + plt.close() + +def heatmaps(name, to_return = False): + indices = {"NDVI": None, "NDRE": None, "SR":None, "SRRE":None, "EVI" : None, "MSAVI" : None} + ds = gdal.Open(name) + if ds == None: + return indices + + #get the band values + if type(GREEN_BAND) == type(1): + green = scaleMinMax(np.array(ds.GetRasterBand(GREEN_BAND).ReadAsArray())) + else: + green = get_bands(GREEN_BAND[0],GREEN_BAND[-1],ds) + + if type(BLUE_BAND) == type(1): + blue = scaleMinMax(np.array(ds.GetRasterBand(BLUE_BAND).ReadAsArray())) + else: + blue = get_bands(BLUE_BAND[0],BLUE_BAND[-1],ds) + + if type(RED_BAND) == type(1): + red = scaleMinMax(np.array(ds.GetRasterBand(RED_BAND).ReadAsArray())) + else: + red = get_bands(RED_BAND[0],RED_BAND[-1],ds) + + if type(RED_EDGE_BAND) == type(1): + red_edge = scaleMinMax(np.array(ds.GetRasterBand(RED_EDGE_BAND).ReadAsArray())) + else: + red_edge = get_bands(RED_EDGE_BAND[0],RED_EDGE_BAND[-1],ds) + + if type(NIR_BAND) == type(1): + nir = scaleMinMax(np.array(ds.GetRasterBand(NIR_BAND).ReadAsArray())) + else: + nir = get_bands(NIR_BAND[0],NIR_BAND[-1],ds) + + indices["NDVI"] = ndvi(nir,red) + indices["NDRE"] = ndre(nir, red_edge) + indices["SR"] = sr(nir, red) + indices["SRRE"] = srre(nir, red_edge) + indices["EVI"] = evi(nir, red, blue) + indices["MSAVI"] = msavi(nir, red) + if to_return: + return indices + for key in indices: + if isinstance(indices[key], np.ndarray): + make_and_save_graph(key,indices[key]) + + +def heatmap_selector(name, selector, to_return = False): + ds = gdal.Open(name) + + """Press 1 for NDVI heatmap + Press 2 for NDRE heatmap, + Press 3 for SR heatmap, + Press 4 for SRRE heatmap """ + + if selector == 1: #NDVI from nir and red + if type(RED_BAND) == type(1): + red = scaleMinMax(np.array(ds.GetRasterBand(RED_BAND).ReadAsArray())) + else: + red = get_bands(RED_BAND[0],RED_BAND[-1],ds) + + if type(NIR_BAND) == type(1): + nir = scaleMinMax(np.array(ds.GetRasterBand(NIR_BAND).ReadAsArray())) + else: + nir = get_bands(NIR_BAND[0],NIR_BAND[-1],ds) + heatmap = ndvi(nir,red) + + elif selector == 2: #NDRE from nir and re + if type(RED_EDGE_BAND) == type(1): + red_edge = scaleMinMax(np.array(ds.GetRasterBand(RED_EDGE_BAND).ReadAsArray())) + else: + red_edge = get_bands(RED_EDGE_BAND[0],RED_EDGE_BAND[-1],ds) + + if type(NIR_BAND) == type(1): + nir = scaleMinMax(np.array(ds.GetRasterBand(NIR_BAND).ReadAsArray())) + else: + nir = get_bands(NIR_BAND[0],NIR_BAND[-1],ds) + heatmap = ndre(nir, red_edge) + + elif selector == 3: # SR from nir and red + if type(RED_BAND) == type(1): + red = scaleMinMax(np.array(ds.GetRasterBand(RED_BAND).ReadAsArray())) + else: + red = get_bands(RED_BAND[0],RED_BAND[-1],ds) + + if type(NIR_BAND) == type(1): + nir = scaleMinMax(np.array(ds.GetRasterBand(NIR_BAND).ReadAsArray())) + else: + nir = get_bands(NIR_BAND[0],NIR_BAND[-1],ds) + heatmap = sr(nir, red) + + elif selector == 4: # SRRE from nir and re + if type(RED_EDGE_BAND) == type(1): + red_edge = scaleMinMax(np.array(ds.GetRasterBand(RED_EDGE_BAND).ReadAsArray())) + else: + red_edge = get_bands(RED_EDGE_BAND[0],RED_EDGE_BAND[-1],ds) + + if type(NIR_BAND) == type(1): + nir = scaleMinMax(np.array(ds.GetRasterBand(NIR_BAND).ReadAsArray())) + else: + nir = get_bands(NIR_BAND[0],NIR_BAND[-1],ds) + heatmap = srre(nir, red_edge) + + if to_return: + return heatmap + + else: + plt.figure() + plt.imshow(heatmap) + plt.savefig() + plt.show() + + + + +if __name__ == '__main__': + heatmaps("example.tif") + arrays = heatmaps("example.tif", to_return=True) #returns the indices as a dictionary of np.arrays. Graphs wont be saved to avoid cluttering +