Merge pull request #30 from gregory-halverson/main

checking distribution of intermediate values

Author: gregory-halverson

breathing_earth_system_simulator/BESS.py

@@ -31,9 +31,7 @@
 from .load_ball_berry_slope_C3 import load_ball_berry_slope_C3
 from .load_ball_berry_slope_C4 import load_ball_berry_slope_C4
 
-GEOS5FP_DOWNLOAD_DIRECTORY = "~/data/GEOS5FP_download"
-
-## TODO automatic acquisition of SRTM for default elevation
+from .check_distribution import check_distribution
 
 def BESS(
         ST_C: Union[Raster, np.ndarray],  # surface temperature in Celsius
@@ -78,7 +76,6 @@ def BESS(
         peakVCmax_C3: np.ndarray = None,  # peak maximum carboxylation rate for C3 plants
         peakVCmax_C4: np.ndarray = None,  # peak maximum carboxylation rate for C4 plants
         CI: Union[Raster, np.ndarray] = None,
-        GEOS5FP_download_directory: str = GEOS5FP_DOWNLOAD_DIRECTORY,
         resampling: str = RESAMPLING):  # clumping index
     if geometry is None and isinstance(ST_C, Raster):
         geometry = ST_C.geometry
@@ -244,6 +241,30 @@ def BESS(
         canopy_height_meters=canopy_height_meters
     )
 
+    meteorology_outputs = {
+        "Ps_Pa": Ps_Pa,
+        "VPD_Pa": VPD_Pa,
+        "RH": RH,
+        "desTa": desTa,
+        "ddesTa": ddesTa,
+        "gamma": gamma,
+        "Cp": Cp,
+        "rhoa": rhoa,
+        "epsa": epsa,
+        "R": R,
+        "Rc": Rc,
+        "Rs": Rs,
+        "SFd": SFd,
+        "SFd2": SFd2,
+        "DL": DL,
+        "Ra": Ra,
+        "fStress": fStress
+    }
+
+    # Check the distribution for each variable
+    for var_name, var_value in meteorology_outputs.items():
+        check_distribution(var_value, var_name, time_UTC)
+
     # convert NDVI to LAI
     LAI = LAI_from_NDVI(NDVI)
     LAI_minimum = LAI_from_NDVI(NDVI_minimum)
@@ -259,6 +280,18 @@ def BESS(
         kn=kn
     )
 
+    # List of variable names and their corresponding values
+    VCmax_outputs = {
+        "VCmax_C3_sunlit": VCmax_C3_sunlit,
+        "VCmax_C4_sunlit": VCmax_C4_sunlit,
+        "VCmax_C3_shaded": VCmax_C3_shaded,
+        "VCmax_C4_shaded": VCmax_C4_shaded
+    }
+
+    # Check the distribution for each variable
+    for var_name, var_value in VCmax_outputs.items():
+        check_distribution(var_value, var_name, time_UTC)
+
     sunlit_fraction, APAR_sunlit, APAR_shaded, ASW_sunlit, ASW_shaded, ASW_soil, G = canopy_shortwave_radiation(
         PARDiff=VISdiff,  # diffuse photosynthetically active radiation in W/m^2
         PARDir=VISdir,  # direct photosynthetically active radiation in W/m^2
@@ -272,6 +305,22 @@ def BESS(
         albedo_NIR=albedo_NIR  # surface albedo in near-infrared wavelengths
     )
 
+    # List of variable names and their corresponding values
+    canopy_radiation_outputs = {
+        "sunlit_fraction": sunlit_fraction,
+        "APAR_sunlit": APAR_sunlit,
+        "APAR_shaded": APAR_shaded,
+        "ASW_sunlit": ASW_sunlit,
+        "ASW_shaded": ASW_shaded,
+        "ASW_soil": ASW_soil,
+        "G": G
+    }
+
+    # Check the distribution for each variable
+    for var_name, var_value in canopy_radiation_outputs.items():
+        check_distribution(var_value, var_name, time_UTC)
+
+
     canopy_temperature_K = canopy_temperature_C + 273.15
     soil_temperature_K = soil_temperature_C + 273.15
 
@@ -309,6 +358,22 @@ def BESS(
         C4_photosynthesis=False  # C3 or C4 photosynthesis
     )
 
+    # List of variable names and their corresponding values
+    carbon_water_fluxes_outputs = {
+        "GPP_C3": GPP_C3,
+        "LE_C3": LE_C3,
+        "LE_soil_C3": LE_soil_C3,
+        "LE_canopy_C3": LE_canopy_C3,
+        "Rn_C3": Rn_C3,
+        "Rn_soil_C3": Rn_soil_C3,
+        "Rn_canopy_C3": Rn_canopy_C3
+    }
+
+    # Check the distribution for each variable
+    for var_name, var_value in carbon_water_fluxes_outputs.items():
+        check_distribution(var_value, var_name, time_UTC)
+
+
     GPP_C4, LE_C4, LE_soil_C4, LE_canopy_C4, Rn_C4, Rn_soil_C4, Rn_canopy_C4 = carbon_water_fluxes(
         canopy_temperature_K=canopy_temperature_K,  # canopy temperature in Kelvin
         soil_temperature_K=soil_temperature_K,  # soil temperature in Kelvin
@@ -343,6 +408,21 @@ def BESS(
         C4_photosynthesis=True  # C3 or C4 photosynthesis
     )
 
+    # List of variable names and their corresponding values
+    carbon_water_fluxes_C4_outputs = {
+        "GPP_C4": GPP_C4,
+        "LE_C4": LE_C4,
+        "LE_soil_C4": LE_soil_C4,
+        "LE_canopy_C4": LE_canopy_C4,
+        "Rn_C4": Rn_C4,
+        "Rn_soil_C4": Rn_soil_C4,
+        "Rn_canopy_C4": Rn_canopy_C4
+    }
+
+    # Check the distribution for each variable
+    for var_name, var_value in carbon_water_fluxes_C4_outputs.items():
+        check_distribution(var_value, var_name, time_UTC)
+
     # interpolate C3 and C4 GPP
     ST_K = ST_C + 273.15
     GPP = np.clip(interpolate_C3_C4(GPP_C3, GPP_C4, C4_fraction), 0, 50)

breathing_earth_system_simulator/check_distribution.py

@@ -0,0 +1,73 @@
+from typing import Union
+from rasters import Raster
+from datetime import date
+import numpy as np
+from logging import getLogger
+import colored_logging as cl
+
+logger = getLogger(__name__)
+
+class BlankOutputError(Exception):
+    pass
+
+def check_distribution(
+        image: Raster,
+        variable: str,
+        date_UTC: Union[date, str],
+        target: str = None):
+    unique = np.unique(image)
+    nan_proportion = np.count_nonzero(np.isnan(image)) / np.size(image)
+
+    target_message = f" at {cl.place(target)}" if target else ""
+
+    if len(unique) < 10:
+        logger.info(f"variable {cl.name(variable)} ({image.dtype}) on {cl.time(f'{date_UTC:%Y-%m-%d}')}{target_message} with {cl.val(unique)} unique values")
+
+        for value in unique:
+            if np.isnan(value):
+                count = np.count_nonzero(np.isnan(image))
+            else:
+                count = np.count_nonzero(image == value)
+
+            if value == 0 or np.isnan(value):
+                logger.info(f"* {cl.colored(value, 'red')}: {cl.colored(count, 'red')}")
+            else:
+                logger.info(f"* {cl.val(value)}: {cl.val(count)}")
+    else:
+        minimum = np.nanmin(image)
+
+        if minimum < 0:
+            minimum_string = cl.colored(f"{minimum:0.3f}", "red")
+        else:
+            minimum_string = cl.val(f"{minimum:0.3f}")
+
+        maximum = np.nanmax(image)
+
+        if maximum <= 0:
+            maximum_string = cl.colored(f"{maximum:0.3f}", "red")
+        else:
+            maximum_string = cl.val(f"{maximum:0.3f}")
+
+        if nan_proportion > 0.5:
+            nan_proportion_string = cl.colored(f"{(nan_proportion * 100):0.2f}%", "yellow")
+        elif nan_proportion == 1:
+            nan_proportion_string = cl.colored(f"{(nan_proportion * 100):0.2f}%", "red")
+        else:
+            nan_proportion_string = cl.val(f"{(nan_proportion * 100):0.2f}%")
+
+        message = "variable " + cl.name(variable) + \
+            " on " + cl.time(f"{date_UTC:%Y-%m-%d}") + \
+            target_message + \
+            " min: " + minimum_string + \
+            " mean: " + cl.val(f"{np.nanmean(image):0.3f}") + \
+            " max: " + maximum_string + \
+            " nan: " + nan_proportion_string + f" ({cl.val(image.nodata)})"
+
+        if np.all(image == 0):
+            message += " all zeros"
+            logger.warning(message)
+        else:
+            logger.info(message)
+
+    if nan_proportion == 1:
+        raise BlankOutputError(f"variable {variable} on {date_UTC:%Y-%m-%d}{target_message} is a blank image")

breathing_earth_system_simulator/version.txt

@@ -1 +1 @@
-1.5.0
+1.5.1

pyproject.toml

@@ -3,7 +3,7 @@ requires = ["setuptools>=60", "setuptools-scm>=8.0", "wheel"]
 
 [project]
 name = "breathing-earth-system-simulator"
-version = "1.5.0"
+version = "1.5.1"
 description = "Breathing Earth System Simulator (BESS) Gross Primary Production (GPP) and Evapotranspiration (ET) Model Python"
 readme = "README.md"
 authors = [