Merge pull request #39 from gregory-halverson/main

daylight upscaling

Author: gregory-halverson-jpl

PTJPLSM/constants.py

@@ -1,6 +1,6 @@
 from PTJPL.constants import *
 
-UPSCALE_TO_DAILY = False
+UPSCALE_TO_DAYLIGHT = False
 
 FIELD_CAPACITY_SCALE = 0.7
 RESAMPLING = "cubic"

PTJPLSM/model.py

@@ -42,8 +42,8 @@
 
 from carlson_leaf_area_index import carlson_leaf_area_index
 from sun_angles import calculate_daylight
-from verma_net_radiation import verma_net_radiation, daily_Rn_integration_verma
-from daily_evapotranspiration_upscaling import daily_ET_from_instantaneous, daily_ET_from_daily_LE
+from verma_net_radiation import verma_net_radiation, daylight_Rn_integration_verma
+from daylight_evapotranspiration import daylight_ET_from_instantaneous_LE, daylight_ET_from_daylight_LE
 
 from PTJPL import GAMMA_PA
 from PTJPL import BETA_PA
@@ -77,7 +77,7 @@
 def PTJPLSM(
         NDVI: Union[Raster, np.ndarray, float],
         Rn_Wm2: Union[Raster, np.ndarray, float] = None,
-        Rn_daily_Wm2: Union[Raster, np.ndarray, float] = None,
+        Rn_daylight_Wm2: Union[Raster, np.ndarray, float] = None,
         geometry: RasterGeometry = None,
         time_UTC: datetime = None,
         hour_of_day: np.ndarray = None,
@@ -107,7 +107,7 @@ def PTJPLSM(
         wilting_point_directory: str = SOIL_CAPACITY_DIRECTORY,
         canopy_height_directory: str = GEDI_DOWNLOAD_DIRECTORY,
         floor_Topt: bool = FLOOR_TOPT,
-        upscale_to_daily: bool = UPSCALE_TO_DAILY,
+        upscale_to_daylight: bool = UPSCALE_TO_DAYLIGHT,
         regenerate_net_radiation: bool = False,
         resampling: str = RESAMPLING) -> Dict[str, Union[Raster, np.ndarray]]:
     """
@@ -307,8 +307,8 @@ def PTJPLSM(
         elif SWin_Wm2 is not None:
             logger.info("using given shortwave radiation (SWin_Wm2)")
 
-        if upscale_to_daily:
-            logger.info("running Verma net radiation with daily upscaling")
+        if upscale_to_daylight:
+            logger.info("running Verma net radiation with daylight upscaling")
         else:
             logger.info("running instantaneous Verma net radiation")
 
@@ -319,13 +319,13 @@ def PTJPLSM(
             emissivity=emissivity,
             Ta_C=Ta_C,
             RH=RH,
-            upscale_to_daily=upscale_to_daily,
+            upscale_to_daylight=upscale_to_daylight,
         )
 
         Rn_Wm2 = Rn_results["Rn_Wm2"]
 
-        if "Rn_daily_Wm2" in Rn_results:
-            Rn_daily_Wm2 = Rn_results["Rn_daily_Wm2"]
+        if "Rn_daylight_Wm2" in Rn_results:
+            Rn_daylight_Wm2 = Rn_results["Rn_daylight_Wm2"]
 
     elif Rn_Wm2 is not None:
         logger.info("using given net radiation (Rn_Wm2) for PT-JPL-SM processing")
@@ -465,49 +465,25 @@ def PTJPLSM(
     check_distribution(LE_Wm2, "LE_Wm2")
     results["LE_Wm2"] = LE_Wm2
 
-    if upscale_to_daily and time_UTC is not None:
-        logger.info("started daily ET upscaling")
+    if upscale_to_daylight and time_UTC is not None:
+        logger.info("started daylight ET upscaling")
         t_et = TicToc()
         t_et.tic()
 
-        if Rn_daily_Wm2 is None:
-            logger.info("running daily net radiation integration")
-            Rn_daily_Wm2 = daily_Rn_integration_verma(
-                Rn_Wm2=Rn_Wm2,
-                time_UTC=time_UTC,
-                geometry=geometry
-            )
-
-        check_distribution(Rn_daily_Wm2, "Rn_daily_Rn_Wm2")
-        results["Rn_daily_Wm2"] = Rn_daily_Wm2
-
-        EF = rt.where((LE_Wm2 == 0) | ((Rn_Wm2 - G_Wm2) == 0), 0, LE_Wm2 / (Rn_Wm2 - G_Wm2))
-        check_distribution(EF, "EF")
-        results["EF"] = EF
-
-        # Calculate latent heat flux during daylight
-        LE_daylight_Wm2 = EF * Rn_daily_Wm2
-        check_distribution(LE_daylight_Wm2, "LE_daylight_Wm2")
-        results["LE_daylight_Wm2"] = LE_daylight_Wm2
-
-        # Calculate daily ET
-        # ET = daily_ET_from_daily_LE(LE_daylight_Wm2, datetime_UTC=time_UTC, geometry=geometry)
-
-        daylight_hours = calculate_daylight(day_of_year=day_of_year, time_UTC=time_UTC, geometry=geometry)
-
-        # convert length of day in hours to seconds
-        daylight_seconds = daylight_hours * 3600.0
-
-        LAMBDA_JKG_WATER_20C = 2450000.0
-
-        # factor seconds out of watts to get joules and divide by latent heat of vaporization to get kilograms
-        ET_daily_kg = rt.clip(LE_daylight_Wm2 * daylight_seconds / LAMBDA_JKG_WATER_20C, 0.0, None)
-
-        check_distribution(ET_daily_kg, "ET_daily_kg")
-        results["ET_daily_kg"] = ET_daily_kg
+        # Use new upscaling function from daylight_evapotranspiration
+        daylight_results = daylight_ET_from_instantaneous_LE(
+            LE_instantaneous_Wm2=LE_Wm2,
+            Rn_instantaneous_Wm2=Rn_Wm2,
+            G_instantaneous_Wm2=G_Wm2,
+            day_of_year=day_of_year,
+            time_UTC=time_UTC,
+            geometry=geometry
+        )
+        # Add all returned daylight results to output
+        results.update(daylight_results)
 
         elapsed_et = t_et.tocvalue()
-        logger.info(f"completed daily ET upscaling (elapsed: {elapsed_et:.2f} seconds)")
+        logger.info(f"completed daylight ET upscaling (elapsed: {elapsed_et:.2f} seconds)")
 
     elapsed = t.tocvalue()
     logger.info(f"PT-JPL-SM model run complete (elapsed: {elapsed:.2f} seconds)")

PTJPLSM/process_PTJPLSM_table.py

@@ -28,7 +28,7 @@
 
 def process_PTJPLSM_table(
         input_df: DataFrame,
-        upscale_to_daily: bool = UPSCALE_TO_DAILY,
+        upscale_to_daylight: bool = UPSCALE_TO_DAYLIGHT,
         regenerate_net_radiation: bool = False
         ) -> DataFrame:
     """
@@ -245,7 +245,7 @@ def parse_geom(s):
         RH=RH,
         soil_moisture=soil_moisture,
         Rn_Wm2=Rn_Wm2,
-        Rn_daily_Wm2=Rn_daily_Wm2,
+        Rn_daylight_Wm2=Rn_daily_Wm2,
         Topt_C=Topt_C,
         fAPARmax=fAPARmax,
         canopy_height_meters=canopy_height_meters,
@@ -256,7 +256,7 @@ def parse_geom(s):
         SWin_Wm2=SWin_Wm2,
         time_UTC=time_UTC,
         regenerate_net_radiation=regenerate_net_radiation,
-        upscale_to_daily=upscale_to_daily
+        upscale_to_daylight=upscale_to_daylight
     )
 
     output_df = input_df.copy()