From 659ed79a16b27894985b8420de1e0f34d176cc2b Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 8 Oct 2025 11:46:23 -0600 Subject: [PATCH 01/30] Envphys - Lint albedo code; fix typos --- smrf/envphys/albedo.py | 78 +++++++++++++++++++----------------------- 1 file changed, 36 insertions(+), 42 deletions(-) diff --git a/smrf/envphys/albedo.py b/smrf/envphys/albedo.py index 62efae8d..97493d42 100644 --- a/smrf/envphys/albedo.py +++ b/smrf/envphys/albedo.py @@ -1,4 +1,5 @@ import math +from typing import Tuple import numpy as np @@ -29,50 +30,42 @@ def growth(t): """ a = 4.0 - b = 3. + b = 3.0 c = 2.0 d = 1.0 - factor = (a+(b*t)+(t*t))/(c+(d*t)+(t*t)) - 1.0 + factor = (a + (b * t) + (t * t)) / (c + (d * t) + (t * t)) - 1.0 - return(1.0 - factor) + return 1.0 - factor -def albedo(telapsed, cosz, gsize, maxgsz, dirt=2): +def albedo(telapsed, cosz, gsize, maxgsz, dirt=2) -> Tuple[np.ndarray, np.ndarray]: """ - Calculate the abedo, adapted from IPW function albedo + Calculate the albedo, adapted from IPW function albedo Args: - telapsed - time since last snow storm (decimal days) - cosz - cosine local solar illumination angle matrix - gsize - gsize is effective grain radius of snow after last storm (mu m) - maxgsz - maxgsz is maximum grain radius expected from grain growth - (mu m) - dirt - dirt is effective contamination for adjustment to visible - albedo (usually between 1.5-3.0) + telapsed: Time since last snow storm (decimal days) + cosz: Cosine local solar illumination angle matrix + gsize: Grain size is effective grain radius of snow after last storm (mu m) + maxgsz: Max grain size is maximum grain radius expected from grain growth (mu m) + dirt: Dirt is effective contamination for adjustment to visible albedo + (usually between 1.5-3.0) Returns: - tuple: - Returns a tuple containing the visible and IR spectral albedo - - - **alb_v** (*numpy.array*) - albedo for visible specturm + Tuple + alb_v - albedo for visible spectrum + alb_ir - albedo for ir spectrum - - **alb_ir** (*numpy.array*) - albedo for ir spectrum - - Created April 17, 2015 Modified July 23, 2015 - take image of cosz and calculate albedo for one time step Scott Havens """ - -# telapsed = np.array(telapsed) - # check inputs if gsize <= 0 or gsize > 500: raise Exception("unrealistic input: gsize=%i", gsize) - if (maxgsz <= gsize or maxgsz > 2000): + if maxgsz <= gsize or maxgsz > 2000: raise Exception("unrealistic input: maxgsz=%i", maxgsz) if 1 >= dirt >= 10: @@ -87,11 +80,11 @@ def albedo(telapsed, cosz, gsize, maxgsz, dirt=2): # calc grain growth decay factor growth_factor = growth(telapsed + 1.0) - # calc effective gsizes for vis & ir + # calc effective grain size for vis & ir gv = radius_v + (range_v * growth_factor) gir = radius_ir + (range_ir * growth_factor) - # calc albedos for cos(z)=1 + # calc albedo for cos(z)=1 alb_v_1 = MAXV - (gv / VFAC) alb_ir_1 = MAXIR * np.exp(IRFAC * gir) @@ -113,10 +106,9 @@ def albedo(telapsed, cosz, gsize, maxgsz, dirt=2): return alb_v, alb_ir -def decay_alb_power(veg, veg_type, start_decay, end_decay, - t_curr, pwr, alb_v, alb_ir): +def decay_alb_power(veg, veg_type, start_decay, end_decay, t_curr, pwr, alb_v, alb_ir): """ - Find a decrease in albedo due to litter acccumulation. Decay is based on + Find a decrease in albedo due to litter accumulation. Decay is based on max decay, decay power, and start and end dates. No litter decay occurs before start_date. Fore times between start and end of decay, @@ -133,7 +125,7 @@ def decay_alb_power(veg, veg_type, start_decay, end_decay, end_decay: date at which to end albedo decay curve (datetime) t_curr: datetime object of current timestep pwr: power for power law decay - alb_v: numpy array of albedo for visibile spectrum + alb_v: numpy array of albedo for visible spectrum alb_ir: numpy array of albedo for IR spectrum Returns: @@ -151,13 +143,15 @@ def decay_alb_power(veg, veg_type, start_decay, end_decay, """ # Calculate hour past start of decay t_diff_hr = t_curr - start_decay - t_diff_hr = t_diff_hr.days*24.0 + \ - t_diff_hr.seconds/3600.0 # only need hours here + t_diff_hr = ( + t_diff_hr.days * 24.0 + t_diff_hr.seconds / 3600.0 + ) # only need hours here # Calculate total time of decay - t_decay_hr = (end_decay - start_decay) - t_decay_hr = t_decay_hr.days*24.0 + \ - t_decay_hr.seconds/3600.0 # only need hours here + t_decay_hr = end_decay - start_decay + t_decay_hr = ( + t_decay_hr.days * 24.0 + t_decay_hr.seconds / 3600.0 + ) # only need hours here # correct for veg alb_dec = np.zeros_like(alb_v) @@ -169,7 +163,7 @@ def decay_alb_power(veg, veg_type, start_decay, end_decay, # Use max decay if after start elif t_diff_hr > t_decay_hr: # Use default - alb_dec = alb_dec + veg['default'] + alb_dec = alb_dec + veg["default"] # Decay based on veg type for k, v in veg.items(): if isint(k): @@ -178,20 +172,20 @@ def decay_alb_power(veg, veg_type, start_decay, end_decay, # Power function decay if during decay period else: # Use defaults - max_dec = veg['default'] - tao = (t_decay_hr) / (max_dec**(1.0/pwr)) + max_dec = veg["default"] + tao = (t_decay_hr) / (max_dec ** (1.0 / pwr)) # Add default decay to array of zeros - alb_dec = alb_dec + ((t_diff_hr) / tao)**pwr + alb_dec = alb_dec + ((t_diff_hr) / tao) ** pwr # Decay based on veg type for k, v in veg.items(): max_dec = v - tao = (t_decay_hr) / (max_dec**(1.0/pwr)) + tao = (t_decay_hr) / (max_dec ** (1.0 / pwr)) # Set albedo decay at correct veg types if isint(k): - alb_dec[veg_type == int(k)] = ((t_diff_hr) / tao)**pwr + alb_dec[veg_type == int(k)] = ((t_diff_hr) / tao) ** pwr alb_v_d = alb_v - alb_dec alb_ir_d = alb_ir - alb_dec @@ -201,7 +195,7 @@ def decay_alb_power(veg, veg_type, start_decay, end_decay, def decay_alb_hardy(litter, veg_type, storm_day, alb_v, alb_ir): """ - Find a decrease in albedo due to litter acccumulation + Find a decrease in albedo due to litter accumulation using method from :cite:`Hardy:2000` with storm_day as input. .. math:: @@ -220,7 +214,7 @@ def decay_alb_hardy(litter, veg_type, storm_day, alb_v, alb_ir): litter: A dictionary of values for default,albedo,41,42,43 veg types veg_type: An image of the basin's NLCD veg type storm_day: numpy array of decimal day since last storm - alb_v: numpy array of albedo for visibile spectrum + alb_v: numpy array of albedo for visible spectrum alb_ir: numpy array of albedo for IR spectrum alb_litter: albedo of pure litter From 013432fae4e6a22e9476a5591fc1d91ceb78e52f Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 8 Oct 2025 17:20:07 -0600 Subject: [PATCH 02/30] Albedo - Rename decay start and end date config and break out method. Rename the parameters to define a decay window and move the function to determine whether the window is present to the distribute class for re-use. This prepares add a new post fire based decay method that should also be triggered within the same time window. The config parameters are now: * start_decay (formerly date_method_start_decay) * end_decay (formerly date_method_end_decay) --- smrf/distribute/albedo.py | 24 ++++++++ smrf/envphys/albedo.py | 57 +++++++------------ smrf/framework/CoreConfig.ini | 10 ++-- smrf/framework/changelog.ini | 4 +- smrf/framework/model_framework.py | 10 ++-- .../basins/Lakes/gold_hrrr/gold_config.ini | 4 +- smrf/tests/basins/RME/gold/gold_config.ini | 4 +- .../basins/RME/gold_hrrr/gold_config.ini | 4 +- 8 files changed, 61 insertions(+), 56 deletions(-) diff --git a/smrf/distribute/albedo.py b/smrf/distribute/albedo.py index 7efdf0fd..d82bafc5 100755 --- a/smrf/distribute/albedo.py +++ b/smrf/distribute/albedo.py @@ -1,3 +1,4 @@ +from typing import Tuple from datetime import datetime import numpy as np @@ -194,3 +195,26 @@ def distribute( self.albedo_vis = utils.set_min_max(alb_v, self.min, self.max) self.albedo_ir = utils.set_min_max(alb_ir, self.min, self.max) + + else: + self.albedo_vis = np.zeros(storm_day.shape) + self.albedo_ir = np.zeros(storm_day.shape) + + def decay_window(self, current_timestep: datetime) -> Tuple[float, float]: + # Calculate hour past start of decay + current_difference = current_timestep - self.config["start_decay"] + current_hours = ( + current_difference.days * 24.0 + current_difference.seconds / 3600.0 + ) + + # Exit if we are before the window starts + if current_hours < 0: + return -1, 0 + + # Calculate total time of decay + decay_difference = self.config["end_decay"] - self.config["start_decay"] + decay_hours = ( + decay_difference.days * 24.0 + decay_difference.seconds / 3600.0 + ) + + return current_hours, decay_hours diff --git a/smrf/envphys/albedo.py b/smrf/envphys/albedo.py index 97493d42..f6000de7 100644 --- a/smrf/envphys/albedo.py +++ b/smrf/envphys/albedo.py @@ -106,7 +106,15 @@ def albedo(telapsed, cosz, gsize, maxgsz, dirt=2) -> Tuple[np.ndarray, np.ndarra return alb_v, alb_ir -def decay_alb_power(veg, veg_type, start_decay, end_decay, t_curr, pwr, alb_v, alb_ir): +def decay_alb_power( + veg: dict, + veg_type: np.ndarray, + current_hours: float, + decay_hours: float, + pwr: float, + alb_v: np.ndarray, + alb_ir: np.ndarray, +) -> Tuple[np.ndarray, np.ndarray]: """ Find a decrease in albedo due to litter accumulation. Decay is based on max decay, decay power, and start and end dates. No litter decay occurs @@ -121,47 +129,22 @@ def decay_alb_power(veg, veg_type, start_decay, end_decay, t_curr, pwr, alb_v, a and :math:`end` are the current, start, and end times for the litter decay. Args: - start_decay: date to start albedo decay (datetime) - end_decay: date at which to end albedo decay curve (datetime) - t_curr: datetime object of current timestep + veg: Vegetation specific decay values + veg_type: Array of vegetation type from the topo + current_hours: Time delta in hours of current time minus decay start time + decay_hours: Time delta in hours of decay start time to end time pwr: power for power law decay alb_v: numpy array of albedo for visible spectrum alb_ir: numpy array of albedo for IR spectrum - Returns: - tuple: - Returns a tuple containing the corrected albedo arrays - based on date, veg type - - **alb_v** (*numpy.array*) - albedo for visible specturm - - - **alb_ir** (*numpy.array*) - albedo for ir spectrum - - - Created July 18, 2017 - Micah Sandusky + Returns: Tuple + alb_v_d, alb_ir_d : numpy arrays of decayed albedo """ - # Calculate hour past start of decay - t_diff_hr = t_curr - start_decay - t_diff_hr = ( - t_diff_hr.days * 24.0 + t_diff_hr.seconds / 3600.0 - ) # only need hours here - - # Calculate total time of decay - t_decay_hr = end_decay - start_decay - t_decay_hr = ( - t_decay_hr.days * 24.0 + t_decay_hr.seconds / 3600.0 - ) # only need hours here - - # correct for veg alb_dec = np.zeros_like(alb_v) - # Don't decay if before start - if t_diff_hr <= 0.0: - alb_dec = alb_dec * 0.0 - # Use max decay if after start - elif t_diff_hr > t_decay_hr: + if current_hours > decay_hours: # Use default alb_dec = alb_dec + veg["default"] # Decay based on veg type @@ -173,19 +156,19 @@ def decay_alb_power(veg, veg_type, start_decay, end_decay, t_curr, pwr, alb_v, a else: # Use defaults max_dec = veg["default"] - tao = (t_decay_hr) / (max_dec ** (1.0 / pwr)) + tao = current_hours / (max_dec ** (1.0 / pwr)) # Add default decay to array of zeros - alb_dec = alb_dec + ((t_diff_hr) / tao) ** pwr + alb_dec = alb_dec + (current_hours / tao) ** pwr # Decay based on veg type for k, v in veg.items(): max_dec = v - tao = (t_decay_hr) / (max_dec ** (1.0 / pwr)) + tao = decay_hours / (max_dec ** (1.0 / pwr)) # Set albedo decay at correct veg types if isint(k): - alb_dec[veg_type == int(k)] = ((t_diff_hr) / tao) ** pwr + alb_dec[veg_type == int(k)] = (current_hours / tao) ** pwr alb_v_d = alb_v - alb_dec alb_ir_d = alb_ir - alb_dec diff --git a/smrf/framework/CoreConfig.ini b/smrf/framework/CoreConfig.ini index f0a4f438..be706faf 100755 --- a/smrf/framework/CoreConfig.ini +++ b/smrf/framework/CoreConfig.ini @@ -794,17 +794,15 @@ default = None, options = [ hardy2000 date_method None], description = Describe how the albedo decays in the late season -date_method_start_decay : +start_decay : default = None, type = DatetimeOrderedPair, -description = Starting date for applying the decay method described by - date_method +description = Starting date for applying decay based methods -date_method_end_decay : +end_decay : default = None, type = DatetimeOrderedPair, -description = Starting date for applying the decay method described by - date_method +description = Starting date for applying decay based methods date_method_decay_power : default = 0.714, diff --git a/smrf/framework/changelog.ini b/smrf/framework/changelog.ini index f5705096..482e1571 100644 --- a/smrf/framework/changelog.ini +++ b/smrf/framework/changelog.ini @@ -107,8 +107,8 @@ precip/min_drift -> precip/winstral_min_drift precip/max_drift -> precip/winstral_max_drift # ALBEDO -albedo/start_decay -> albedo/date_method_start_decay -albedo/end_decay -> albedo/date_method_end_decay +albedo/date_method_start_decay -> albedo/start_decay +albedo/date_method_end_decay -> albedo/end_decay albedo/decay_power -> albedo/date_method_decay_power albedo/veg_default -> albedo/date_method_veg_default albedo/veg_41 -> albedo/date_method_veg_41 diff --git a/smrf/framework/model_framework.py b/smrf/framework/model_framework.py index 1ef68ab8..b4338424 100755 --- a/smrf/framework/model_framework.py +++ b/smrf/framework/model_framework.py @@ -117,13 +117,13 @@ def __init__(self, config, external_logger=None): self._setup_date_and_time() # need to align date time - if "date_method_start_decay" in self.config[Albedo.DISTRIBUTION_KEY].keys(): - self.config[Albedo.DISTRIBUTION_KEY]["date_method_start_decay"] = self.config[ + if self.config[Albedo.DISTRIBUTION_KEY].get("start_decay", None): + self.config[Albedo.DISTRIBUTION_KEY]["start_decay"] = self.config[ Albedo.DISTRIBUTION_KEY - ]["date_method_start_decay"].replace(tzinfo=self.time_zone) - self.config[Albedo.DISTRIBUTION_KEY]["date_method_end_decay"] = self.config[ + ]["start_decay"].replace(tzinfo=self.time_zone) + self.config[Albedo.DISTRIBUTION_KEY]["end_decay"] = self.config[ Albedo.DISTRIBUTION_KEY - ]["date_method_end_decay"].replace(tzinfo=self.time_zone) + ]["end_decay"].replace(tzinfo=self.time_zone) # Gridded dataset self.forecast_flag = False diff --git a/smrf/tests/basins/Lakes/gold_hrrr/gold_config.ini b/smrf/tests/basins/Lakes/gold_hrrr/gold_config.ini index a7c3435b..3ceaa36d 100644 --- a/smrf/tests/basins/Lakes/gold_hrrr/gold_config.ini +++ b/smrf/tests/basins/Lakes/gold_hrrr/gold_config.ini @@ -135,8 +135,8 @@ max_grain: 700.0 dirt: 2.0 decay_method: None grid_mask: True -date_method_start_decay: None -date_method_end_decay: None +start_decay: None +end_decay: None date_method_decay_power: 0.714 date_method_veg_default: 0.25 date_method_veg_41: 0.36 diff --git a/smrf/tests/basins/RME/gold/gold_config.ini b/smrf/tests/basins/RME/gold/gold_config.ini index bdfddb1f..9d4dda38 100644 --- a/smrf/tests/basins/RME/gold/gold_config.ini +++ b/smrf/tests/basins/RME/gold/gold_config.ini @@ -149,8 +149,8 @@ max_grain: 700.0 dirt: 2.0 decay_method: None grid_mask: True -date_method_start_decay: None -date_method_end_decay: None +start_decay: None +end_decay: None date_method_decay_power: 0.714 date_method_veg_default: 0.25 date_method_veg_41: 0.36 diff --git a/smrf/tests/basins/RME/gold_hrrr/gold_config.ini b/smrf/tests/basins/RME/gold_hrrr/gold_config.ini index b59948bf..b4c9db33 100644 --- a/smrf/tests/basins/RME/gold_hrrr/gold_config.ini +++ b/smrf/tests/basins/RME/gold_hrrr/gold_config.ini @@ -145,8 +145,8 @@ max_grain: 2000.0 dirt: 2.0 decay_method: None grid_mask: True -date_method_start_decay: None -date_method_end_decay: None +start_decay: None +end_decay: None date_method_decay_power: 0.714 date_method_veg_default: 0.25 date_method_veg_41: 0.36 From f5a99de7996bc84b5dd167ee2d4bc3b05baccb66 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 8 Oct 2025 17:24:28 -0600 Subject: [PATCH 03/30] Topo - Read in a burn mask if present. Adds reading in a new layer that defines a burn mask on the topo file. --- smrf/data/load_topo.py | 33 ++++++++++++++++----------------- smrf/tests/data/test_topo.py | 30 ++++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+), 17 deletions(-) diff --git a/smrf/data/load_topo.py b/smrf/data/load_topo.py index 1864233a..58f9ac21 100755 --- a/smrf/data/load_topo.py +++ b/smrf/data/load_topo.py @@ -23,19 +23,18 @@ class GdalAttributes: class Topo: """ - Class for topo images and processing those images. Images are: - - DEM - - Mask - - veg type - - veg height - - veg k - - veg tau - - Inputs to topo are the topo section of the config file - + Class reading the configured metadata for a model domain. """ - IMAGES = ["dem", "mask", "veg_type", "veg_height", "veg_k", "veg_tau"] + IMAGES = [ + "burn_mask", + "dem", + "mask", + "veg_height", + "veg_k", + "veg_tau", + "veg_type", + ] def __init__(self, topoConfig): self.topoConfig = topoConfig @@ -151,16 +150,16 @@ def readImages(self, f): f: netcdf dataset object """ # netCDF files are stored typically as 32-bit float, so convert to double or int - for v_smrf in self.IMAGES: + for image in self.IMAGES: result = None - if v_smrf in f.variables.keys(): - if v_smrf == "veg_type": - result = f.variables[v_smrf][:].astype(int) + if image in f.variables.keys(): + if image == "veg_type": + result = f.variables[image][:].astype(int) else: - result = f.variables[v_smrf][:].astype(np.float64) + result = f.variables[image][:].astype(np.float64) - setattr(self, v_smrf, result) + setattr(self, image, result) def get_center(self, ds, mask_name=None): """ diff --git a/smrf/tests/data/test_topo.py b/smrf/tests/data/test_topo.py index 6f9a837c..8eed6f11 100644 --- a/smrf/tests/data/test_topo.py +++ b/smrf/tests/data/test_topo.py @@ -29,11 +29,41 @@ def setUp(cls): def tearDown(cls): cls.ds.close() +<<<<<<< HEAD @mock.patch.object(Topo, "gradient") @mock.patch.object(Topo, "readNetCDF") def test_init(self, read_nc, gradient): topo = Topo(TOPO_CONFIG) self.assertEqual(TOPO_CONFIG, topo.topoConfig) +======= + def test_read_topo_images(self): + self.assertEqual( + [ + "burn_mask", + "dem", + "mask", + "veg_height", + "veg_k", + "veg_tau", + "veg_type", + ], + Topo.IMAGES + ) + + def test_topo_image_as_variables(self): + for variable in Topo.IMAGES: + self.assertTrue(hasattr(self.topo, variable)) + + def test_no_burn_mask_required(self): + self.assertIsNone(self.topo.burn_mask) + + @mock.patch.object(Topo, 'readNetCDF') + @mock.patch.object(Topo, 'gradient') + def test_init(self, gradient, read_nc): + topo = Topo(self.TOPO_CONFIG) + self.assertEqual(self.TOPO_CONFIG, topo.topoConfig) + gradient.assert_called_once() +>>>>>>> 190ae1b (Topo - Read in a burn mask if present.) read_nc.assert_called_once() gradient.assert_called_once() From cf4c41333d3aa8203df0185eb4cd33c0dfedd17b Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 8 Oct 2025 17:26:31 -0600 Subject: [PATCH 04/30] Distribute - Albedo - Lint code --- smrf/distribute/albedo.py | 25 +++++++++---------------- 1 file changed, 9 insertions(+), 16 deletions(-) diff --git a/smrf/distribute/albedo.py b/smrf/distribute/albedo.py index d82bafc5..6d3eb860 100755 --- a/smrf/distribute/albedo.py +++ b/smrf/distribute/albedo.py @@ -1,4 +1,6 @@ +from datetime import datetime from typing import Tuple + from datetime import datetime import numpy as np @@ -117,11 +119,8 @@ def distribute( Args: current_time_step: Current time step in datetime object - cosz: numpy array of the illumination angle for the current time - step - storm_day: numpy array of the decimal days since it last - snowed at a grid cell - + cosz: Llumination angle for the current time step + storm_day: Decimal days since it last snowed at a grid cell """ self._logger.debug("%s Distributing albedo" % current_time_step) @@ -171,28 +170,22 @@ def distribute( # Perform litter decay if self.config["decay_method"] == "date_method": - alb_v_d, alb_ir_d = albedo.decay_alb_power( + current_hours, decay_hours = self.decay_window(current_time_step) + alb_v, alb_ir = albedo.decay_alb_power( self.veg, self.veg_type, - self.config["date_method_start_decay"], - self.config["date_method_end_decay"], - current_time_step, + current_hours, + decay_hours, self.config["date_method_decay_power"], alb_v, alb_ir, ) - alb_v = alb_v_d - alb_ir = alb_ir_d - elif self.config["decay_method"] == "hardy2000": - alb_v_d, alb_ir_d = albedo.decay_alb_hardy( + alb_v, alb_ir = albedo.decay_alb_hardy( self.litter, self.veg_type, storm_day, alb_v, alb_ir ) - alb_v = alb_v_d - alb_ir = alb_ir_d - self.albedo_vis = utils.set_min_max(alb_v, self.min, self.max) self.albedo_ir = utils.set_min_max(alb_ir, self.min, self.max) From 9e358b0482075ac85724fc745e7fd3ebfb0e7ab0 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 8 Oct 2025 17:31:56 -0600 Subject: [PATCH 05/30] Albedo - Add new decay function based on binary burn mask Add a new option to decay albedo based on a burn mask in the topo file and configured values in the .ini via `post_fire_k_burned` and `post_fire_k_unburned`. The method itself will be triggered via the decay method `post_fire`. Co-authored-by: caldwellng --- smrf/distribute/albedo.py | 15 ++++- smrf/envphys/albedo.py | 44 +++++++++++++ smrf/framework/CoreConfig.ini | 16 ++++- smrf/tests/distribute/test_albedo.py | 92 ++++++++++++++++++++++++++++ 4 files changed, 164 insertions(+), 3 deletions(-) diff --git a/smrf/distribute/albedo.py b/smrf/distribute/albedo.py index 6d3eb860..8237d34d 100755 --- a/smrf/distribute/albedo.py +++ b/smrf/distribute/albedo.py @@ -1,8 +1,6 @@ from datetime import datetime from typing import Tuple -from datetime import datetime - import numpy as np import numpy.typing as npt @@ -92,6 +90,8 @@ def initialize(self, metadata): """ super().initialize(metadata) + self.burn_mask = self.topo.burn_mask + if ( self.config.get("decay_method", None) is None and self.config["source_files"] is None @@ -180,6 +180,17 @@ def distribute( alb_v, alb_ir, ) + elif self.config["decay_method"] == "post_fire": + current_hours, decay_hours = self.decay_window(current_time_step) + if current_hours > 0: + alb_v, alb_ir = albedo.decay_burned( + alb_v, + alb_ir, + storm_day, + self.burn_mask, + self.config["post_fire_k_burned"], + self.config["post_fire_k_unburned"], + ) elif self.config["decay_method"] == "hardy2000": alb_v, alb_ir = albedo.decay_alb_hardy( diff --git a/smrf/envphys/albedo.py b/smrf/envphys/albedo.py index f6000de7..eb511ce0 100644 --- a/smrf/envphys/albedo.py +++ b/smrf/envphys/albedo.py @@ -237,3 +237,47 @@ def decay_alb_hardy(litter, veg_type, storm_day, alb_v, alb_ir): alb_ir_d = alb_ir*sc + alb_litter*lc return alb_v_d, alb_ir_d + +def decay_burned( + alb_v: np.ndarray, + alb_ir: np.ndarray, + last_snow: float, + burn_mask: np.ndarray, + k_burned: float, + k_unburned: float, +) -> Tuple[np.ndarray, np.ndarray]: + """ + Apply exponential albedo decay as a function of time since last snowfall. + Different decay rates are applied for burned and unburned pixels. + + Args: + alb_v: Visible albedo + alb_ir: Infrared albedo + last_snow: Time since last snow storm (decimal days) + burn_mask: Mask of burned area + k_burned: Decay rate for burned area + k_unburned: Decay rate unburned area + + Returns: + alb_v_d, alb_ir_d : numpy arrays of decayed albedo + """ + # initialize output + alb_v_d = np.copy(alb_v) + alb_ir_d = np.copy(alb_ir) + + # masks + burned = burn_mask == 1 + unburned = burn_mask == 0 + + # exponential decay factors depending on burn condition + decay_burned = np.exp(-k_burned * last_snow) + decay_unburned = np.exp(-k_unburned * last_snow) + + # apply decay rates to vis and infrared albedo + alb_v_d[burned] = alb_v[burned] * decay_burned[burned] + alb_ir_d[burned] = alb_ir[burned] * decay_burned[burned] + + alb_v_d[unburned] = alb_v[unburned] * decay_unburned[unburned] + alb_ir_d[unburned] = alb_ir[unburned] * decay_unburned[unburned] + + return alb_v_d, alb_ir_d diff --git a/smrf/framework/CoreConfig.ini b/smrf/framework/CoreConfig.ini index be706faf..664a8d85 100755 --- a/smrf/framework/CoreConfig.ini +++ b/smrf/framework/CoreConfig.ini @@ -791,7 +791,7 @@ description = Effective contamination for adjustment to visible albedo (usually decay_method : default = None, -options = [ hardy2000 date_method None], +options = [hardy2000 date_method post_fire None], description = Describe how the albedo decays in the late season start_decay : @@ -857,12 +857,26 @@ type = float, description = Litter rate for places where vegetation not specified for Hardy et al. 2000 decay method for vegetation classes NLCD 43 +<<<<<<< HEAD source_files : default = None, type = Directory, description = Full path to netcdf files that contain albedo values. Once set none of the other albedo variables will be effective and the raw values from the files will be used. +======= +post_fire_k_burned: +default = None, +type = float, +description = Decay power for areas masked as 1 (burned) in the burned mask of the + topo file + +post_fire_k_unburned: +default = None, +type = float, +description = Decay power for areas masked as 0 (unburned) in the burned mask of the + topo file +>>>>>>> 2da0b90 (Albedo - Add new decay function based on binary burn mask) ################################################################################ # Cloud Factor - Fraction used to limit solar radiation Cloudy (0) - Sunny (1) diff --git a/smrf/tests/distribute/test_albedo.py b/smrf/tests/distribute/test_albedo.py index e51258cc..6fb6ab95 100644 --- a/smrf/tests/distribute/test_albedo.py +++ b/smrf/tests/distribute/test_albedo.py @@ -1,4 +1,5 @@ import unittest +<<<<<<< HEAD from unittest.mock import MagicMock, patch import numpy as np @@ -135,3 +136,94 @@ def test_distribute_calculated(self, mock_albedo_calc, mock_albedo_decay): mock_albedo_decay.assert_called_once() np.testing.assert_array_equal(subject.albedo_vis, values_vis_decay) np.testing.assert_array_equal(subject.albedo_ir, values_ir_decay) +======= +from unittest.mock import patch, MagicMock + +import numpy as np +import pandas as pd + +from smrf.distribute.albedo import Albedo + +CONFIG = { + "decay_method": "date_method", + "start_decay": pd.to_datetime("2025-04-01"), + "end_decay": pd.to_datetime("2025-07-01"), + "grain_size": 100.0, + "max_grain": 800.0, + "max": 1.0, + "min": 0.0, + "dirt": 2, + "date_method_veg_default": 0.2, +} +DECAY_TIME = pd.to_datetime("2025-04-30") +NO_DECAY_TIME = pd.to_datetime("2025-03-01") +STORM_DAYS = np.array([[0.0, 1.0], [1.0, 0.0]]) +TOPO = MagicMock( + veg_type=MagicMock("veg_type"), +) +DATA = MagicMock() + +ALBEDO_VIS = np.array([[0.0, 1.0], [1.0, 0.0]]) +ALBEDO_IR = np.array([[0.0, 1.0], [1.0, 0.0]]) + + +class TestAlbedo(unittest.TestCase): + def setUp(self): + self.subject = Albedo(CONFIG) + self.subject.initialize(TOPO, DATA) + + def test_before_decay_window(self): + current, decay = self.subject.decay_window(NO_DECAY_TIME) + self.assertEqual(-1, current) + + def test_in_decay_window(self): + current, decay = self.subject.decay_window(DECAY_TIME) + self.assertEqual(696.0, current) + self.assertEqual(2184.0, decay) + + @patch('smrf.distribute.albedo.albedo.albedo', return_value=(ALBEDO_VIS, ALBEDO_IR)) + @patch('smrf.distribute.albedo.albedo.decay_alb_power') + def test_distribute_date_method_not_in_window(self, decay_alb_power, envphys_albedo): + self.subject.distribute(NO_DECAY_TIME, 1, STORM_DAYS) + + envphys_albedo.assert_called() + envphys_albedo.assert_called_with( + STORM_DAYS, 1, CONFIG["grain_size"], CONFIG["max_grain"], CONFIG["dirt"] + ) + decay_alb_power.assert_not_called() + + @patch('smrf.distribute.albedo.albedo.albedo', return_value=(ALBEDO_VIS, ALBEDO_IR)) + @patch('smrf.distribute.albedo.albedo.decay_alb_power', return_value=(ALBEDO_VIS, ALBEDO_IR)) + def test_distribute_date_method_in_window(self, decay_alb_power, _envphys_albedo): + self.subject.config["date_method_decay_power"] = 0.7 + self.subject.config["date_method_veg_default"] = 0.2 + self.subject.distribute(DECAY_TIME, 1, STORM_DAYS) + + decay_alb_power.assert_called_with( + {"default": 0.2}, TOPO.veg_type, 696.0, 2184.0, 0.7, ALBEDO_VIS, ALBEDO_IR + ) + + @patch('smrf.distribute.albedo.albedo.albedo', return_value=(ALBEDO_VIS, ALBEDO_IR)) + @patch('smrf.distribute.albedo.albedo.decay_burned') + def test_distribute_post_fire_not_in_window(self, decay_burned, envphys_albedo): + self.subject.config["decay_method"] = "post_fire" + self.subject.distribute(NO_DECAY_TIME, 1, STORM_DAYS) + + envphys_albedo.assert_called() + envphys_albedo.assert_called_with( + STORM_DAYS, 1, CONFIG["grain_size"], CONFIG["max_grain"], CONFIG["dirt"] + ) + decay_burned.assert_not_called() + + @patch('smrf.distribute.albedo.albedo.albedo', return_value=(ALBEDO_VIS, ALBEDO_IR)) + @patch('smrf.distribute.albedo.albedo.decay_burned', return_value=(ALBEDO_VIS, ALBEDO_IR)) + def test_distribute_post_fire_in_window(self, decay_burned, _envphys_albedo): + self.subject.config["decay_method"] = "post_fire" + self.subject.config["post_fire_k_burned"] = 0.06 + self.subject.config["post_fire_k_unburned"] = 0.02 + self.subject.distribute(DECAY_TIME, 1, STORM_DAYS) + + decay_burned.assert_called_with( + ALBEDO_VIS, ALBEDO_IR, STORM_DAYS, TOPO.burn_mask, 0.06, 0.02 + ) +>>>>>>> 2da0b90 (Albedo - Add new decay function based on binary burn mask) From 61da0ee4bb95a38f1e2297e64aaed257a9946a59 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 8 Oct 2025 17:32:31 -0600 Subject: [PATCH 06/30] Precip - Use failsafe access to config key for storm_days_restart --- smrf/distribute/precipitation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/smrf/distribute/precipitation.py b/smrf/distribute/precipitation.py index fd9b0614..ae28c2b8 100755 --- a/smrf/distribute/precipitation.py +++ b/smrf/distribute/precipitation.py @@ -106,7 +106,7 @@ def initialize(self, data): self.storm_total = np.zeros((self.topo.ny, self.topo.nx)) # Assign storm_days array if given - if self.config["storm_days_restart"] is not None: + if self.config.get("storm_days_restart", None): self._logger.debug( "Reading {} from {}".format( "storm_days", self.config["storm_days_restart"] From 1b93edb1240fb16ea31d7df76362e0bcd62e1331 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Thu, 9 Oct 2025 19:33:04 -0600 Subject: [PATCH 07/30] Envphys - Small tweaks to decay_burned function * Fix type hint for `last_snow` parameter; it's an array * Improve output initialization and use empty values * Rename exponential variables that conflicted with method name --- smrf/envphys/albedo.py | 18 +++++----- smrf/tests/envphys/test_albedo.py | 55 +++++++++++++++++++++++++++++++ 2 files changed, 64 insertions(+), 9 deletions(-) create mode 100644 smrf/tests/envphys/test_albedo.py diff --git a/smrf/envphys/albedo.py b/smrf/envphys/albedo.py index eb511ce0..02aabdd2 100644 --- a/smrf/envphys/albedo.py +++ b/smrf/envphys/albedo.py @@ -241,7 +241,7 @@ def decay_alb_hardy(litter, veg_type, storm_day, alb_v, alb_ir): def decay_burned( alb_v: np.ndarray, alb_ir: np.ndarray, - last_snow: float, + last_snow: np.ndarray, burn_mask: np.ndarray, k_burned: float, k_unburned: float, @@ -262,22 +262,22 @@ def decay_burned( alb_v_d, alb_ir_d : numpy arrays of decayed albedo """ # initialize output - alb_v_d = np.copy(alb_v) - alb_ir_d = np.copy(alb_ir) + alb_v_d = np.empty_like(alb_v) + alb_ir_d = np.empty_like(alb_ir) # masks burned = burn_mask == 1 unburned = burn_mask == 0 # exponential decay factors depending on burn condition - decay_burned = np.exp(-k_burned * last_snow) - decay_unburned = np.exp(-k_unburned * last_snow) + burned_exp = np.exp(-k_burned * last_snow) + unburned_exp = np.exp(-k_unburned * last_snow) # apply decay rates to vis and infrared albedo - alb_v_d[burned] = alb_v[burned] * decay_burned[burned] - alb_ir_d[burned] = alb_ir[burned] * decay_burned[burned] + alb_v_d[burned] = alb_v[burned] * burned_exp[burned] + alb_ir_d[burned] = alb_ir[burned] * burned_exp[burned] - alb_v_d[unburned] = alb_v[unburned] * decay_unburned[unburned] - alb_ir_d[unburned] = alb_ir[unburned] * decay_unburned[unburned] + alb_v_d[unburned] = alb_v[unburned] * unburned_exp[unburned] + alb_ir_d[unburned] = alb_ir[unburned] * unburned_exp[unburned] return alb_v_d, alb_ir_d diff --git a/smrf/tests/envphys/test_albedo.py b/smrf/tests/envphys/test_albedo.py new file mode 100644 index 00000000..d258429b --- /dev/null +++ b/smrf/tests/envphys/test_albedo.py @@ -0,0 +1,55 @@ +import unittest + +import numpy as np + +from smrf.envphys.albedo import decay_burned + +ALBEDO_VIS = np.array([[0.8, 0.6], [0.7, 0.5]]) +ALBEDO_IR = np.array([[0.7, 0.5], [0.6, 0.4]]) +LAST_SNOW = np.array([[1.0, 1.5], [0.0, 2.0]]) +K_BURNED = 0.06 +K_UNBURNED = 0.02 + + +class TestDecayBurned(unittest.TestCase): + def test_decay_burned_with_burned_and_unburned_areas(self): + burn_mask = np.array([[1, 0], [0, 1]]) + + alb_vis, alb_ir = decay_burned( + ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, K_BURNED, K_UNBURNED + ) + + np.testing.assert_array_almost_equal( + np.array([[0.753412, 0.582267], [0.7, 0.44346]]), alb_vis, decimal=6 + ) + np.testing.assert_array_almost_equal( + np.array([[0.659235, 0.485223], [0.6, 0.354768]]), alb_ir, decimal=6 + ) + + def test_decay_burned_all_burned(self): + burn_mask = np.ones_like(ALBEDO_VIS) + + alb_vis, alb_ir = decay_burned( + ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, K_BURNED, K_UNBURNED + ) + + np.testing.assert_array_almost_equal( + np.array([[0.753412, 0.548359], [0.7, 0.44346]]), alb_vis, decimal=6 + ) + np.testing.assert_array_almost_equal( + np.array([[0.659235, 0.456966], [0.6, 0.354768]]), alb_ir, decimal=6 + ) + + def test_decay_burned_all_unburned(self): + burn_mask = np.zeros_like(ALBEDO_VIS) + + alb_vis, alb_ir = decay_burned( + ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, K_BURNED, K_UNBURNED + ) + + np.testing.assert_array_almost_equal( + np.array([[0.784159, 0.582267], [0.7, 0.480395]]), alb_vis, decimal=6 + ) + np.testing.assert_array_almost_equal( + np.array([[0.686139, 0.485223], [0.6, 0.384316]]), alb_ir, decimal=6 + ) From 17c1e32b78d2748d058d0f5af7b25f45b9ea3b52 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Fri, 10 Oct 2025 08:43:58 -0600 Subject: [PATCH 08/30] CoreConfig - Albedo - Change the name for dates one more time Change the names to configure dates to `decay_start` and `decay_end` to be consistent with the main parameter name `decay_method`. --- smrf/distribute/albedo.py | 4 ++-- smrf/framework/CoreConfig.ini | 4 ++-- smrf/framework/changelog.ini | 4 ++-- smrf/framework/model_framework.py | 10 +++++----- smrf/framework/recipes.ini | 4 ++-- smrf/tests/basins/Lakes/gold_hrrr/gold_config.ini | 4 ++-- smrf/tests/basins/RME/gold/gold_config.ini | 4 ++-- smrf/tests/basins/RME/gold_hrrr/gold_config.ini | 4 ++-- smrf/tests/distribute/test_albedo.py | 4 ++-- 9 files changed, 21 insertions(+), 21 deletions(-) diff --git a/smrf/distribute/albedo.py b/smrf/distribute/albedo.py index 8237d34d..1f9f277b 100755 --- a/smrf/distribute/albedo.py +++ b/smrf/distribute/albedo.py @@ -206,7 +206,7 @@ def distribute( def decay_window(self, current_timestep: datetime) -> Tuple[float, float]: # Calculate hour past start of decay - current_difference = current_timestep - self.config["start_decay"] + current_difference = current_timestep - self.config["decay_start"] current_hours = ( current_difference.days * 24.0 + current_difference.seconds / 3600.0 ) @@ -216,7 +216,7 @@ def decay_window(self, current_timestep: datetime) -> Tuple[float, float]: return -1, 0 # Calculate total time of decay - decay_difference = self.config["end_decay"] - self.config["start_decay"] + decay_difference = self.config["decay_end"] - self.config["decay_start"] decay_hours = ( decay_difference.days * 24.0 + decay_difference.seconds / 3600.0 ) diff --git a/smrf/framework/CoreConfig.ini b/smrf/framework/CoreConfig.ini index 664a8d85..009d45d5 100755 --- a/smrf/framework/CoreConfig.ini +++ b/smrf/framework/CoreConfig.ini @@ -794,12 +794,12 @@ default = None, options = [hardy2000 date_method post_fire None], description = Describe how the albedo decays in the late season -start_decay : +decay_start : default = None, type = DatetimeOrderedPair, description = Starting date for applying decay based methods -end_decay : +decay_end : default = None, type = DatetimeOrderedPair, description = Starting date for applying decay based methods diff --git a/smrf/framework/changelog.ini b/smrf/framework/changelog.ini index 482e1571..fd929554 100644 --- a/smrf/framework/changelog.ini +++ b/smrf/framework/changelog.ini @@ -107,8 +107,8 @@ precip/min_drift -> precip/winstral_min_drift precip/max_drift -> precip/winstral_max_drift # ALBEDO -albedo/date_method_start_decay -> albedo/start_decay -albedo/date_method_end_decay -> albedo/end_decay +albedo/date_method_start_decay -> albedo/decay_start +albedo/date_method_end_decay -> albedo/decay_end albedo/decay_power -> albedo/date_method_decay_power albedo/veg_default -> albedo/date_method_veg_default albedo/veg_41 -> albedo/date_method_veg_41 diff --git a/smrf/framework/model_framework.py b/smrf/framework/model_framework.py index b4338424..408c22b2 100755 --- a/smrf/framework/model_framework.py +++ b/smrf/framework/model_framework.py @@ -117,13 +117,13 @@ def __init__(self, config, external_logger=None): self._setup_date_and_time() # need to align date time - if self.config[Albedo.DISTRIBUTION_KEY].get("start_decay", None): - self.config[Albedo.DISTRIBUTION_KEY]["start_decay"] = self.config[ + if self.config[Albedo.DISTRIBUTION_KEY].get("decay_start", None): + self.config[Albedo.DISTRIBUTION_KEY]["decay_start"] = self.config[ Albedo.DISTRIBUTION_KEY - ]["start_decay"].replace(tzinfo=self.time_zone) - self.config[Albedo.DISTRIBUTION_KEY]["end_decay"] = self.config[ + ]["decay_start"].replace(tzinfo=self.time_zone) + self.config[Albedo.DISTRIBUTION_KEY]["decay_end"] = self.config[ Albedo.DISTRIBUTION_KEY - ]["end_decay"].replace(tzinfo=self.time_zone) + ]["decay_end"].replace(tzinfo=self.time_zone) # Gridded dataset self.forecast_flag = False diff --git a/smrf/framework/recipes.ini b/smrf/framework/recipes.ini index 55aa5628..7cdbcba6 100755 --- a/smrf/framework/recipes.ini +++ b/smrf/framework/recipes.ini @@ -107,7 +107,7 @@ trigger: has_value = [albedo decay_method None] albedo: - remove_item = [date_method_start_decay date_method_end_decay + remove_item = [decay_start decay_end date_method_decay_power date_method_veg_default date_method_veg_41 date_method_veg_42 date_method_veg_43 hardy2000_litter_albedo hardy2000_litter_default @@ -126,7 +126,7 @@ trigger: has_value = [albedo decay_method hardy2000] albedo: - remove_item = [date_method_start_decay date_method_end_decay + remove_item = [decay_start decay_end date_method_decay_power date_method_veg_default date_method_veg_41 date_method_veg_42 date_method_veg_43] diff --git a/smrf/tests/basins/Lakes/gold_hrrr/gold_config.ini b/smrf/tests/basins/Lakes/gold_hrrr/gold_config.ini index 3ceaa36d..b206f38d 100644 --- a/smrf/tests/basins/Lakes/gold_hrrr/gold_config.ini +++ b/smrf/tests/basins/Lakes/gold_hrrr/gold_config.ini @@ -135,8 +135,8 @@ max_grain: 700.0 dirt: 2.0 decay_method: None grid_mask: True -start_decay: None -end_decay: None +decay_start: None +decay_end: None date_method_decay_power: 0.714 date_method_veg_default: 0.25 date_method_veg_41: 0.36 diff --git a/smrf/tests/basins/RME/gold/gold_config.ini b/smrf/tests/basins/RME/gold/gold_config.ini index 9d4dda38..d96cece6 100644 --- a/smrf/tests/basins/RME/gold/gold_config.ini +++ b/smrf/tests/basins/RME/gold/gold_config.ini @@ -149,8 +149,8 @@ max_grain: 700.0 dirt: 2.0 decay_method: None grid_mask: True -start_decay: None -end_decay: None +decay_start: None +decay_end: None date_method_decay_power: 0.714 date_method_veg_default: 0.25 date_method_veg_41: 0.36 diff --git a/smrf/tests/basins/RME/gold_hrrr/gold_config.ini b/smrf/tests/basins/RME/gold_hrrr/gold_config.ini index b4c9db33..3396075d 100644 --- a/smrf/tests/basins/RME/gold_hrrr/gold_config.ini +++ b/smrf/tests/basins/RME/gold_hrrr/gold_config.ini @@ -145,8 +145,8 @@ max_grain: 2000.0 dirt: 2.0 decay_method: None grid_mask: True -start_decay: None -end_decay: None +decay_start: None +decay_end: None date_method_decay_power: 0.714 date_method_veg_default: 0.25 date_method_veg_41: 0.36 diff --git a/smrf/tests/distribute/test_albedo.py b/smrf/tests/distribute/test_albedo.py index 6fb6ab95..1971ced8 100644 --- a/smrf/tests/distribute/test_albedo.py +++ b/smrf/tests/distribute/test_albedo.py @@ -146,8 +146,8 @@ def test_distribute_calculated(self, mock_albedo_calc, mock_albedo_decay): CONFIG = { "decay_method": "date_method", - "start_decay": pd.to_datetime("2025-04-01"), - "end_decay": pd.to_datetime("2025-07-01"), + "decay_start": pd.to_datetime("2025-04-01"), + "decay_end": pd.to_datetime("2025-07-01"), "grain_size": 100.0, "max_grain": 800.0, "max": 1.0, From 14be402618057d4fe65d6008cbdb6d5bdbbb5d5e Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Fri, 10 Oct 2025 08:58:37 -0600 Subject: [PATCH 09/30] Tests - Distribute Albedo - Use proper data type for fake cos_z --- smrf/tests/distribute/test_albedo.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/smrf/tests/distribute/test_albedo.py b/smrf/tests/distribute/test_albedo.py index 1971ced8..19353558 100644 --- a/smrf/tests/distribute/test_albedo.py +++ b/smrf/tests/distribute/test_albedo.py @@ -158,6 +158,7 @@ def test_distribute_calculated(self, mock_albedo_calc, mock_albedo_decay): DECAY_TIME = pd.to_datetime("2025-04-30") NO_DECAY_TIME = pd.to_datetime("2025-03-01") STORM_DAYS = np.array([[0.0, 1.0], [1.0, 0.0]]) +COS_Z = np.array([[10.0, 10.0], [10.0, 10.0]]) TOPO = MagicMock( veg_type=MagicMock("veg_type"), ) @@ -184,11 +185,11 @@ def test_in_decay_window(self): @patch('smrf.distribute.albedo.albedo.albedo', return_value=(ALBEDO_VIS, ALBEDO_IR)) @patch('smrf.distribute.albedo.albedo.decay_alb_power') def test_distribute_date_method_not_in_window(self, decay_alb_power, envphys_albedo): - self.subject.distribute(NO_DECAY_TIME, 1, STORM_DAYS) + self.subject.distribute(NO_DECAY_TIME, COS_Z, STORM_DAYS) envphys_albedo.assert_called() envphys_albedo.assert_called_with( - STORM_DAYS, 1, CONFIG["grain_size"], CONFIG["max_grain"], CONFIG["dirt"] + STORM_DAYS, COS_Z, CONFIG["grain_size"], CONFIG["max_grain"], CONFIG["dirt"] ) decay_alb_power.assert_not_called() @@ -197,7 +198,7 @@ def test_distribute_date_method_not_in_window(self, decay_alb_power, envphys_alb def test_distribute_date_method_in_window(self, decay_alb_power, _envphys_albedo): self.subject.config["date_method_decay_power"] = 0.7 self.subject.config["date_method_veg_default"] = 0.2 - self.subject.distribute(DECAY_TIME, 1, STORM_DAYS) + self.subject.distribute(DECAY_TIME, COS_Z, STORM_DAYS) decay_alb_power.assert_called_with( {"default": 0.2}, TOPO.veg_type, 696.0, 2184.0, 0.7, ALBEDO_VIS, ALBEDO_IR @@ -207,11 +208,11 @@ def test_distribute_date_method_in_window(self, decay_alb_power, _envphys_albedo @patch('smrf.distribute.albedo.albedo.decay_burned') def test_distribute_post_fire_not_in_window(self, decay_burned, envphys_albedo): self.subject.config["decay_method"] = "post_fire" - self.subject.distribute(NO_DECAY_TIME, 1, STORM_DAYS) + self.subject.distribute(NO_DECAY_TIME, COS_Z, STORM_DAYS) envphys_albedo.assert_called() envphys_albedo.assert_called_with( - STORM_DAYS, 1, CONFIG["grain_size"], CONFIG["max_grain"], CONFIG["dirt"] + STORM_DAYS, COS_Z, CONFIG["grain_size"], CONFIG["max_grain"], CONFIG["dirt"] ) decay_burned.assert_not_called() @@ -221,7 +222,7 @@ def test_distribute_post_fire_in_window(self, decay_burned, _envphys_albedo): self.subject.config["decay_method"] = "post_fire" self.subject.config["post_fire_k_burned"] = 0.06 self.subject.config["post_fire_k_unburned"] = 0.02 - self.subject.distribute(DECAY_TIME, 1, STORM_DAYS) + self.subject.distribute(DECAY_TIME, COS_Z, STORM_DAYS) decay_burned.assert_called_with( ALBEDO_VIS, ALBEDO_IR, STORM_DAYS, TOPO.burn_mask, 0.06, 0.02 From 0922d6028b24773f06d513eeb92f7e8ea6cc24b9 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Fri, 10 Oct 2025 09:00:22 -0600 Subject: [PATCH 10/30] Distribute - Albedo - Rename variables to be PEP 8 conform and readable --- smrf/distribute/albedo.py | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/smrf/distribute/albedo.py b/smrf/distribute/albedo.py index 1f9f277b..1fdb391f 100755 --- a/smrf/distribute/albedo.py +++ b/smrf/distribute/albedo.py @@ -52,12 +52,6 @@ class Albedo(VariableBase): ) def __init__(self, config: dict, topo): - """ - Initialize albedo() - - Args: - config: configuration from [albedo] section - """ super().__init__(config, topo) # Variables for calculations of external files @@ -99,7 +93,7 @@ def initialize(self, metadata): self._logger.warning("No decay method is set!") def distribute( - self, current_time_step: datetime, cosz: npt.NDArray, storm_day: npt.NDArray + self, current_time_step: datetime, cos_z: npt.NDArray, storm_day: npt.NDArray ): """ Calculate or load snow albedo for given time step. When calculated, the configured @@ -119,13 +113,13 @@ def distribute( Args: current_time_step: Current time step in datetime object - cosz: Llumination angle for the current time step + cos_z: Illumination angle for the current time step storm_day: Decimal days since it last snowed at a grid cell """ self._logger.debug("%s Distributing albedo" % current_time_step) # Only calculated when the sun is up - if cosz is not None: + if cos_z is not None: if self.source_files is not None: # For files with VIS and IR values (typically a previous run) if ( @@ -162,7 +156,7 @@ def distribute( else: alb_v, alb_ir = albedo.albedo( storm_day, - cosz, + cos_z, self.config["grain_size"], self.config["max_grain"], self.config["dirt"], From ea0f0f6c5656548dc5707d7e3c21bc50e71af9cf Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 15 Oct 2025 18:02:15 -0600 Subject: [PATCH 11/30] Distribute - Albedo - Follow up fixes after merge with latest in main The main branch had an overhaul of the init method, changing a few pieces of the logic in this branch with was developed in parallel. This now also adds the accessor to veg_type from the topo with the ImageData base class. --- smrf/distribute/variable_base.py | 4 ++++ smrf/tests/distribute/test_albedo.py | 36 ++++++++++++++++++---------- 2 files changed, 27 insertions(+), 13 deletions(-) diff --git a/smrf/distribute/variable_base.py b/smrf/distribute/variable_base.py index ad8c9ebc..f526896e 100755 --- a/smrf/distribute/variable_base.py +++ b/smrf/distribute/variable_base.py @@ -142,6 +142,10 @@ def veg_tau(self): def veg_k(self): return self.topo.veg_k + @property + def veg_type(self): + return self.topo.veg_type + # END - Topo accessor methods # START - config accessors diff --git a/smrf/tests/distribute/test_albedo.py b/smrf/tests/distribute/test_albedo.py index 19353558..ccf80b1e 100644 --- a/smrf/tests/distribute/test_albedo.py +++ b/smrf/tests/distribute/test_albedo.py @@ -145,15 +145,17 @@ def test_distribute_calculated(self, mock_albedo_calc, mock_albedo_decay): from smrf.distribute.albedo import Albedo CONFIG = { - "decay_method": "date_method", - "decay_start": pd.to_datetime("2025-04-01"), - "decay_end": pd.to_datetime("2025-07-01"), - "grain_size": 100.0, - "max_grain": 800.0, - "max": 1.0, - "min": 0.0, - "dirt": 2, - "date_method_veg_default": 0.2, + "albedo": { + "decay_method": "date_method", + "decay_start": pd.to_datetime("2025-04-01"), + "decay_end": pd.to_datetime("2025-07-01"), + "grain_size": 100.0, + "max_grain": 800.0, + "max": 1.0, + "min": 0.0, + "dirt": 2, + "date_method_veg_default": 0.2, + } } DECAY_TIME = pd.to_datetime("2025-04-30") NO_DECAY_TIME = pd.to_datetime("2025-03-01") @@ -170,8 +172,8 @@ def test_distribute_calculated(self, mock_albedo_calc, mock_albedo_decay): class TestAlbedo(unittest.TestCase): def setUp(self): - self.subject = Albedo(CONFIG) - self.subject.initialize(TOPO, DATA) + self.subject = Albedo(CONFIG, TOPO) + self.subject.initialize(DATA) def test_before_decay_window(self): current, decay = self.subject.decay_window(NO_DECAY_TIME) @@ -189,7 +191,11 @@ def test_distribute_date_method_not_in_window(self, decay_alb_power, envphys_alb envphys_albedo.assert_called() envphys_albedo.assert_called_with( - STORM_DAYS, COS_Z, CONFIG["grain_size"], CONFIG["max_grain"], CONFIG["dirt"] + STORM_DAYS, + COS_Z, + CONFIG["albedo"]["grain_size"], + CONFIG["albedo"]["max_grain"], + CONFIG["albedo"]["dirt"] ) decay_alb_power.assert_not_called() @@ -212,7 +218,11 @@ def test_distribute_post_fire_not_in_window(self, decay_burned, envphys_albedo): envphys_albedo.assert_called() envphys_albedo.assert_called_with( - STORM_DAYS, COS_Z, CONFIG["grain_size"], CONFIG["max_grain"], CONFIG["dirt"] + STORM_DAYS, + COS_Z, + CONFIG["albedo"]["grain_size"], + CONFIG["albedo"]["max_grain"], + CONFIG["albedo"]["dirt"] ) decay_burned.assert_not_called() From 2dd4475f47d0dc64854a7e9b5876e19ae09fe7d8 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 18 Feb 2026 10:32:49 -0700 Subject: [PATCH 12/30] Copilot review - Remove duplicat code; fix description strings --- smrf/distribute/variable_base.py | 4 ---- smrf/framework/CoreConfig.ini | 2 +- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/smrf/distribute/variable_base.py b/smrf/distribute/variable_base.py index f526896e..ad8c9ebc 100755 --- a/smrf/distribute/variable_base.py +++ b/smrf/distribute/variable_base.py @@ -142,10 +142,6 @@ def veg_tau(self): def veg_k(self): return self.topo.veg_k - @property - def veg_type(self): - return self.topo.veg_type - # END - Topo accessor methods # START - config accessors diff --git a/smrf/framework/CoreConfig.ini b/smrf/framework/CoreConfig.ini index 009d45d5..f9d8e87c 100755 --- a/smrf/framework/CoreConfig.ini +++ b/smrf/framework/CoreConfig.ini @@ -802,7 +802,7 @@ description = Starting date for applying decay based methods decay_end : default = None, type = DatetimeOrderedPair, -description = Starting date for applying decay based methods +description = Ending date for applying decay based methods date_method_decay_power : default = 0.714, From 39c83f50e69c8f1902b7f6b73619ee66b1c6829a Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 18 Feb 2026 10:33:51 -0700 Subject: [PATCH 13/30] Tests - Fix after merge with latest in main. Setup changed when testing distribution variables. --- smrf/tests/distribute/test_albedo.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/smrf/tests/distribute/test_albedo.py b/smrf/tests/distribute/test_albedo.py index ccf80b1e..66322d6e 100644 --- a/smrf/tests/distribute/test_albedo.py +++ b/smrf/tests/distribute/test_albedo.py @@ -145,6 +145,10 @@ def test_distribute_calculated(self, mock_albedo_calc, mock_albedo_decay): from smrf.distribute.albedo import Albedo CONFIG = { + "time": { + "start_date": pd.to_datetime("2025-10-01"), + "time_zone": "utc", + }, "albedo": { "decay_method": "date_method", "decay_start": pd.to_datetime("2025-04-01"), From dc8bf6c953424e8cc48a26ce19c7e49321aa1ef0 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 18 Feb 2026 11:01:10 -0700 Subject: [PATCH 14/30] Envphys - Albedo - Fix decay rate for default values Bug introduced with the refactor to add burned albedo and used the current passed time versus the decay window hours. Thank you Copilot! Added tests to catch this in future. --- smrf/envphys/albedo.py | 6 +- smrf/tests/envphys/test_albedo.py | 145 ++++++++++++++++++++++++++++-- 2 files changed, 142 insertions(+), 9 deletions(-) diff --git a/smrf/envphys/albedo.py b/smrf/envphys/albedo.py index 02aabdd2..ad0ba094 100644 --- a/smrf/envphys/albedo.py +++ b/smrf/envphys/albedo.py @@ -143,8 +143,10 @@ def decay_alb_power( """ alb_dec = np.zeros_like(alb_v) + if current_hours <= 0: + return alb_v, alb_ir # Use max decay if after start - if current_hours > decay_hours: + elif current_hours > decay_hours: # Use default alb_dec = alb_dec + veg["default"] # Decay based on veg type @@ -156,7 +158,7 @@ def decay_alb_power( else: # Use defaults max_dec = veg["default"] - tao = current_hours / (max_dec ** (1.0 / pwr)) + tao = decay_hours / (max_dec ** (1.0 / pwr)) # Add default decay to array of zeros alb_dec = alb_dec + (current_hours / tao) ** pwr diff --git a/smrf/tests/envphys/test_albedo.py b/smrf/tests/envphys/test_albedo.py index d258429b..a05965bb 100644 --- a/smrf/tests/envphys/test_albedo.py +++ b/smrf/tests/envphys/test_albedo.py @@ -1,8 +1,9 @@ import unittest import numpy as np +import numpy.testing as npt -from smrf.envphys.albedo import decay_burned +from smrf.envphys.albedo import decay_alb_power, decay_burned ALBEDO_VIS = np.array([[0.8, 0.6], [0.7, 0.5]]) ALBEDO_IR = np.array([[0.7, 0.5], [0.6, 0.4]]) @@ -11,6 +12,136 @@ K_UNBURNED = 0.02 +class TestDecayAlbPower(unittest.TestCase): + def setUp(self): + self.veg = { + "default": 0.2, + "41": 0.25, + "42": 0.3, + } + self.veg_type = np.array( + [ + [0, 41], + [42, 0], + ] + ) + + def test_no_decay_before_start(self): + current_hours = -10 + decay_hours = 100 + pwr = 1.0 + + alb_v_d, alb_ir_d = decay_alb_power( + self.veg, + self.veg_type, + current_hours, + decay_hours, + pwr, + ALBEDO_VIS, + ALBEDO_IR + ) + + npt.assert_array_equal(ALBEDO_VIS, alb_v_d) + npt.assert_array_equal(ALBEDO_IR, alb_ir_d) + + def test_max_decay_after_end(self): + """Test with current_hours > decay_hours, ensuring maximum decay is applied.""" + current_hours = 150 + decay_hours = 100 + pwr = 1.0 + + expected_v = ALBEDO_VIS - np.array( + [ + [0.2, 0.25], + [0.3, 0.2], + ] + ) + expected_ir = ALBEDO_IR - np.array( + [ + [0.2, 0.25], + [0.3, 0.2], + ] + ) + + alb_v_d, alb_ir_d = decay_alb_power( + self.veg, + self.veg_type, + current_hours, + decay_hours, + pwr, + ALBEDO_VIS, + ALBEDO_IR + ) + + npt.assert_array_almost_equal(expected_v, alb_v_d) + npt.assert_array_almost_equal(expected_ir, alb_ir_d) + + def test_power_decay_during_window(self): + """Test power function decay during the decay period.""" + current_hours = 50 + decay_hours = 100 + pwr = 2.0 + + # Decay rates + tao_default = decay_hours / (self.veg["default"] ** (1.0 / pwr)) + tao_41 = decay_hours / (self.veg["41"] ** (1.0 / pwr)) + tao_42 = decay_hours / (self.veg["42"] ** (1.0 / pwr)) + + expected_decay = np.zeros_like(self.veg_type, dtype=float) + expected_decay[self.veg_type == 0] = (current_hours / tao_default) ** pwr + expected_decay[self.veg_type == 41] = (current_hours / tao_41) ** pwr + expected_decay[self.veg_type == 42] = (current_hours / tao_42) ** pwr + + expected_v = ALBEDO_VIS - expected_decay + expected_ir = ALBEDO_IR - expected_decay + + alb_v_d, alb_ir_d = decay_alb_power( + self.veg, + self.veg_type, + current_hours, + decay_hours, + pwr, + ALBEDO_VIS, + ALBEDO_IR + ) + + npt.assert_array_almost_equal(expected_v, alb_v_d, decimal=5) + npt.assert_array_almost_equal(expected_ir, alb_ir_d, decimal=5) + + def test_partial_coverage_during_decay(self): + """Test decay for mixed vegetation types.""" + current_hours = 75 + decay_hours = 150 + pwr = 1.5 + + # Calculate decay + tao_default = decay_hours / (self.veg["default"] ** (1.0 / pwr)) + decay_default = (current_hours / tao_default) ** pwr + + tao_41 = decay_hours / (self.veg["41"] ** (1.0 / pwr)) + tao_42 = decay_hours / (self.veg["42"] ** (1.0 / pwr)) + + expected_decay = np.zeros_like(self.veg_type, dtype=float) + expected_decay[self.veg_type == 0] = decay_default + expected_decay[self.veg_type == 41] = (current_hours / tao_41) ** pwr + expected_decay[self.veg_type == 42] = (current_hours / tao_42) ** pwr + + expected_v = ALBEDO_VIS - expected_decay + expected_ir = ALBEDO_IR - expected_decay + + alb_v_d, alb_ir_d = decay_alb_power( + self.veg, + self.veg_type, + current_hours, + decay_hours, + pwr, + ALBEDO_VIS, + ALBEDO_IR + ) + + npt.assert_array_almost_equal(expected_v, alb_v_d, decimal=5) + npt.assert_array_almost_equal(expected_ir, alb_ir_d, decimal=5) + class TestDecayBurned(unittest.TestCase): def test_decay_burned_with_burned_and_unburned_areas(self): burn_mask = np.array([[1, 0], [0, 1]]) @@ -19,10 +150,10 @@ def test_decay_burned_with_burned_and_unburned_areas(self): ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, K_BURNED, K_UNBURNED ) - np.testing.assert_array_almost_equal( + npt.assert_array_almost_equal( np.array([[0.753412, 0.582267], [0.7, 0.44346]]), alb_vis, decimal=6 ) - np.testing.assert_array_almost_equal( + npt.assert_array_almost_equal( np.array([[0.659235, 0.485223], [0.6, 0.354768]]), alb_ir, decimal=6 ) @@ -33,10 +164,10 @@ def test_decay_burned_all_burned(self): ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, K_BURNED, K_UNBURNED ) - np.testing.assert_array_almost_equal( + npt.assert_array_almost_equal( np.array([[0.753412, 0.548359], [0.7, 0.44346]]), alb_vis, decimal=6 ) - np.testing.assert_array_almost_equal( + npt.assert_array_almost_equal( np.array([[0.659235, 0.456966], [0.6, 0.354768]]), alb_ir, decimal=6 ) @@ -47,9 +178,9 @@ def test_decay_burned_all_unburned(self): ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, K_BURNED, K_UNBURNED ) - np.testing.assert_array_almost_equal( + npt.assert_array_almost_equal( np.array([[0.784159, 0.582267], [0.7, 0.480395]]), alb_vis, decimal=6 ) - np.testing.assert_array_almost_equal( + npt.assert_array_almost_equal( np.array([[0.686139, 0.485223], [0.6, 0.384316]]), alb_ir, decimal=6 ) From 1d231d414f78b68766425c834d0a89839f3f2ab4 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 18 Feb 2026 12:34:55 -0700 Subject: [PATCH 15/30] Enphys - Albedo - Use numexpr for decay_burned function --- smrf/envphys/albedo.py | 24 +++++++----------------- 1 file changed, 7 insertions(+), 17 deletions(-) diff --git a/smrf/envphys/albedo.py b/smrf/envphys/albedo.py index ad0ba094..acfe76bb 100644 --- a/smrf/envphys/albedo.py +++ b/smrf/envphys/albedo.py @@ -1,4 +1,5 @@ import math +import numexpr as ne from typing import Tuple import numpy as np @@ -263,23 +264,12 @@ def decay_burned( Returns: alb_v_d, alb_ir_d : numpy arrays of decayed albedo """ - # initialize output - alb_v_d = np.empty_like(alb_v) - alb_ir_d = np.empty_like(alb_ir) + decay_factor = ne.evaluate( + "exp(-where(burn_mask == 1, k_burned, k_unburned) * last_snow)" + ) - # masks - burned = burn_mask == 1 - unburned = burn_mask == 0 - - # exponential decay factors depending on burn condition - burned_exp = np.exp(-k_burned * last_snow) - unburned_exp = np.exp(-k_unburned * last_snow) - - # apply decay rates to vis and infrared albedo - alb_v_d[burned] = alb_v[burned] * burned_exp[burned] - alb_ir_d[burned] = alb_ir[burned] * burned_exp[burned] - - alb_v_d[unburned] = alb_v[unburned] * unburned_exp[unburned] - alb_ir_d[unburned] = alb_ir[unburned] * unburned_exp[unburned] + alb_v_d = ne.evaluate("alb_v * decay_factor") + alb_ir_d = ne.evaluate("alb_ir * decay_factor") return alb_v_d, alb_ir_d + From 99ea51f48e77b8edeba88a620833fe28c231b36b Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 18 Feb 2026 14:07:23 -0700 Subject: [PATCH 16/30] Albedo - Fix latest merge with main Couple of conflicts got commited unresolved. This mostly addresses the test class for the distribute/albedo. Also restores some logic improvement when deciding on the distribute method to use. --- smrf/distribute/albedo.py | 19 ++--- smrf/envphys/albedo.py | 6 +- smrf/framework/CoreConfig.ini | 4 +- smrf/tests/distribute/test_albedo.py | 112 ++++++++++++--------------- smrf/tests/envphys/test_albedo.py | 18 ----- 5 files changed, 61 insertions(+), 98 deletions(-) diff --git a/smrf/distribute/albedo.py b/smrf/distribute/albedo.py index 1fdb391f..ae8b9480 100755 --- a/smrf/distribute/albedo.py +++ b/smrf/distribute/albedo.py @@ -165,15 +165,16 @@ def distribute( # Perform litter decay if self.config["decay_method"] == "date_method": current_hours, decay_hours = self.decay_window(current_time_step) - alb_v, alb_ir = albedo.decay_alb_power( - self.veg, - self.veg_type, - current_hours, - decay_hours, - self.config["date_method_decay_power"], - alb_v, - alb_ir, - ) + if current_hours > 0: + alb_v, alb_ir = albedo.decay_alb_power( + self.veg, + self.veg_type, + current_hours, + decay_hours, + self.config["date_method_decay_power"], + alb_v, + alb_ir, + ) elif self.config["decay_method"] == "post_fire": current_hours, decay_hours = self.decay_window(current_time_step) if current_hours > 0: diff --git a/smrf/envphys/albedo.py b/smrf/envphys/albedo.py index acfe76bb..1ad0edb1 100644 --- a/smrf/envphys/albedo.py +++ b/smrf/envphys/albedo.py @@ -119,7 +119,7 @@ def decay_alb_power( """ Find a decrease in albedo due to litter accumulation. Decay is based on max decay, decay power, and start and end dates. No litter decay occurs - before start_date. Fore times between start and end of decay, + before start_date. For times between start and end of decay, .. math:: \\alpha = \\alpha - (dec_{max}^{\\frac{1.0}{pwr}} \\times @@ -144,10 +144,8 @@ def decay_alb_power( """ alb_dec = np.zeros_like(alb_v) - if current_hours <= 0: - return alb_v, alb_ir # Use max decay if after start - elif current_hours > decay_hours: + if current_hours > decay_hours: # Use default alb_dec = alb_dec + veg["default"] # Decay based on veg type diff --git a/smrf/framework/CoreConfig.ini b/smrf/framework/CoreConfig.ini index f9d8e87c..d483ab08 100755 --- a/smrf/framework/CoreConfig.ini +++ b/smrf/framework/CoreConfig.ini @@ -857,14 +857,13 @@ type = float, description = Litter rate for places where vegetation not specified for Hardy et al. 2000 decay method for vegetation classes NLCD 43 -<<<<<<< HEAD source_files : default = None, type = Directory, description = Full path to netcdf files that contain albedo values. Once set none of the other albedo variables will be effective and the raw values from the files will be used. -======= + post_fire_k_burned: default = None, type = float, @@ -876,7 +875,6 @@ default = None, type = float, description = Decay power for areas masked as 0 (unburned) in the burned mask of the topo file ->>>>>>> 2da0b90 (Albedo - Add new decay function based on binary burn mask) ################################################################################ # Cloud Factor - Fraction used to limit solar radiation Cloudy (0) - Sunny (1) diff --git a/smrf/tests/distribute/test_albedo.py b/smrf/tests/distribute/test_albedo.py index 66322d6e..31926746 100644 --- a/smrf/tests/distribute/test_albedo.py +++ b/smrf/tests/distribute/test_albedo.py @@ -1,5 +1,4 @@ import unittest -<<<<<<< HEAD from unittest.mock import MagicMock, patch import numpy as np @@ -20,14 +19,37 @@ }, "albedo": { "decay_method": "date_method", + "decay_start": pd.to_datetime("2025-04-01"), + "decay_end": pd.to_datetime("2025-07-01"), + "grain_size": 100.0, + "max_grain": 800.0, + "max": 1.0, + "min": 0.0, + "dirt": 2, + "date_method_veg_default": 0.2, + "date_method_decay_power": 1, }, } TIMEZONE = pytz.timezone(CONFIG["time"]["time_zone"]) -TIMESTEP = pd.to_datetime("2025-10-01 00:00:00") -STORM_DAY = np.array([0, 1, 0, 0, 5, 0]) +TIMESTEP = pd.to_datetime(CONFIG["time"]["start_date"]) + +DECAY_TIME = pd.to_datetime("2025-04-30") +NO_DECAY_TIME = pd.to_datetime("2025-03-01") + +STORM_DAYS = np.array([[0.0, 1.0], [1.0, 0.0]]) +COS_Z = np.array([[10.0, 10.0], [10.0, 10.0]]) + +DATA = MagicMock() + +ALBEDO_VIS = np.array([[0.0, 1.0], [1.0, 0.0]]) +ALBEDO_IR = np.array([[0.0, 1.0], [1.0, 0.0]]) class TestAlbedo(SMRFConfig, unittest.TestCase): + def setUp(self): + self.subject = Albedo(CONFIG, TOPO) + self.subject.initialize(DATA) + def test_init_default(self): subject = Albedo(CONFIG, TOPO) @@ -55,7 +77,7 @@ def test_distribute_file_broadband(self, mock_read_netcdf): subject = Albedo(config, TOPO) subject.initialize(pd.DataFrame()) - subject.distribute(TIMESTEP, np.ndarray([10]), STORM_DAY) + subject.distribute(TIMESTEP, np.ndarray([10]), STORM_DAYS) npt.assert_equal(subject.albedo, values) mock_source.load.assert_called_with("albedo", TIMESTEP) @@ -74,7 +96,7 @@ def test_distribute_file_vis_ir(self, mock_read_netcdf): subject = Albedo(config, TOPO) subject.initialize(pd.DataFrame()) - subject.distribute(TIMESTEP, np.ndarray([10]), STORM_DAY) + subject.distribute(TIMESTEP, np.ndarray([10]), STORM_DAYS) np.testing.assert_array_equal(subject.albedo_vis, values_vis) np.testing.assert_array_equal(subject.albedo_ir, values_ir) @@ -98,7 +120,7 @@ def test_distribute_file_direct_diffuse(self, mock_read_netcdf): subject = Albedo(config, TOPO) subject.initialize(pd.DataFrame()) - subject.distribute(TIMESTEP, np.ndarray([10]), STORM_DAY) + subject.distribute(TIMESTEP, np.ndarray([10]), STORM_DAYS) np.testing.assert_array_equal(subject.albedo_direct, values_direct) np.testing.assert_array_equal(subject.albedo_diffuse, values_diffuse) @@ -122,62 +144,17 @@ def test_distribute_calculated(self, mock_albedo_calc, mock_albedo_decay): config = self._copy_config(CONFIG) config["albedo"]["grain_size"] = "100" config["albedo"]["max_grain"] = "700" - config["albedo"]["dirt"] = "1.2" - config["albedo"]["date_method_start_decay"] = "2026-04-01" - config["albedo"]["date_method_end_decay"] = "2026-06-30" config["albedo"]["date_method_decay_power"] = "1.2" subject = Albedo(config, TOPO) subject.initialize(pd.DataFrame()) - subject.distribute(TIMESTEP, np.ndarray([10]), STORM_DAY) + subject.distribute(TIMESTEP, np.ndarray([10]), STORM_DAYS) mock_albedo_calc.assert_called_once() mock_albedo_decay.assert_called_once() np.testing.assert_array_equal(subject.albedo_vis, values_vis_decay) np.testing.assert_array_equal(subject.albedo_ir, values_ir_decay) -======= -from unittest.mock import patch, MagicMock - -import numpy as np -import pandas as pd - -from smrf.distribute.albedo import Albedo - -CONFIG = { - "time": { - "start_date": pd.to_datetime("2025-10-01"), - "time_zone": "utc", - }, - "albedo": { - "decay_method": "date_method", - "decay_start": pd.to_datetime("2025-04-01"), - "decay_end": pd.to_datetime("2025-07-01"), - "grain_size": 100.0, - "max_grain": 800.0, - "max": 1.0, - "min": 0.0, - "dirt": 2, - "date_method_veg_default": 0.2, - } -} -DECAY_TIME = pd.to_datetime("2025-04-30") -NO_DECAY_TIME = pd.to_datetime("2025-03-01") -STORM_DAYS = np.array([[0.0, 1.0], [1.0, 0.0]]) -COS_Z = np.array([[10.0, 10.0], [10.0, 10.0]]) -TOPO = MagicMock( - veg_type=MagicMock("veg_type"), -) -DATA = MagicMock() - -ALBEDO_VIS = np.array([[0.0, 1.0], [1.0, 0.0]]) -ALBEDO_IR = np.array([[0.0, 1.0], [1.0, 0.0]]) - - -class TestAlbedo(unittest.TestCase): - def setUp(self): - self.subject = Albedo(CONFIG, TOPO) - self.subject.initialize(DATA) def test_before_decay_window(self): current, decay = self.subject.decay_window(NO_DECAY_TIME) @@ -188,9 +165,11 @@ def test_in_decay_window(self): self.assertEqual(696.0, current) self.assertEqual(2184.0, decay) - @patch('smrf.distribute.albedo.albedo.albedo', return_value=(ALBEDO_VIS, ALBEDO_IR)) - @patch('smrf.distribute.albedo.albedo.decay_alb_power') - def test_distribute_date_method_not_in_window(self, decay_alb_power, envphys_albedo): + @patch("smrf.distribute.albedo.albedo.albedo", return_value=(ALBEDO_VIS, ALBEDO_IR)) + @patch("smrf.distribute.albedo.albedo.decay_alb_power") + def test_distribute_date_method_not_in_window( + self, decay_alb_power, envphys_albedo + ): self.subject.distribute(NO_DECAY_TIME, COS_Z, STORM_DAYS) envphys_albedo.assert_called() @@ -199,12 +178,15 @@ def test_distribute_date_method_not_in_window(self, decay_alb_power, envphys_alb COS_Z, CONFIG["albedo"]["grain_size"], CONFIG["albedo"]["max_grain"], - CONFIG["albedo"]["dirt"] + CONFIG["albedo"]["dirt"], ) decay_alb_power.assert_not_called() - @patch('smrf.distribute.albedo.albedo.albedo', return_value=(ALBEDO_VIS, ALBEDO_IR)) - @patch('smrf.distribute.albedo.albedo.decay_alb_power', return_value=(ALBEDO_VIS, ALBEDO_IR)) + @patch("smrf.distribute.albedo.albedo.albedo", return_value=(ALBEDO_VIS, ALBEDO_IR)) + @patch( + "smrf.distribute.albedo.albedo.decay_alb_power", + return_value=(ALBEDO_VIS, ALBEDO_IR), + ) def test_distribute_date_method_in_window(self, decay_alb_power, _envphys_albedo): self.subject.config["date_method_decay_power"] = 0.7 self.subject.config["date_method_veg_default"] = 0.2 @@ -214,8 +196,8 @@ def test_distribute_date_method_in_window(self, decay_alb_power, _envphys_albedo {"default": 0.2}, TOPO.veg_type, 696.0, 2184.0, 0.7, ALBEDO_VIS, ALBEDO_IR ) - @patch('smrf.distribute.albedo.albedo.albedo', return_value=(ALBEDO_VIS, ALBEDO_IR)) - @patch('smrf.distribute.albedo.albedo.decay_burned') + @patch("smrf.distribute.albedo.albedo.albedo", return_value=(ALBEDO_VIS, ALBEDO_IR)) + @patch("smrf.distribute.albedo.albedo.decay_burned") def test_distribute_post_fire_not_in_window(self, decay_burned, envphys_albedo): self.subject.config["decay_method"] = "post_fire" self.subject.distribute(NO_DECAY_TIME, COS_Z, STORM_DAYS) @@ -226,12 +208,15 @@ def test_distribute_post_fire_not_in_window(self, decay_burned, envphys_albedo): COS_Z, CONFIG["albedo"]["grain_size"], CONFIG["albedo"]["max_grain"], - CONFIG["albedo"]["dirt"] + CONFIG["albedo"]["dirt"], ) decay_burned.assert_not_called() - @patch('smrf.distribute.albedo.albedo.albedo', return_value=(ALBEDO_VIS, ALBEDO_IR)) - @patch('smrf.distribute.albedo.albedo.decay_burned', return_value=(ALBEDO_VIS, ALBEDO_IR)) + @patch("smrf.distribute.albedo.albedo.albedo", return_value=(ALBEDO_VIS, ALBEDO_IR)) + @patch( + "smrf.distribute.albedo.albedo.decay_burned", + return_value=(ALBEDO_VIS, ALBEDO_IR), + ) def test_distribute_post_fire_in_window(self, decay_burned, _envphys_albedo): self.subject.config["decay_method"] = "post_fire" self.subject.config["post_fire_k_burned"] = 0.06 @@ -241,4 +226,3 @@ def test_distribute_post_fire_in_window(self, decay_burned, _envphys_albedo): decay_burned.assert_called_with( ALBEDO_VIS, ALBEDO_IR, STORM_DAYS, TOPO.burn_mask, 0.06, 0.02 ) ->>>>>>> 2da0b90 (Albedo - Add new decay function based on binary burn mask) diff --git a/smrf/tests/envphys/test_albedo.py b/smrf/tests/envphys/test_albedo.py index a05965bb..9b64b090 100644 --- a/smrf/tests/envphys/test_albedo.py +++ b/smrf/tests/envphys/test_albedo.py @@ -26,24 +26,6 @@ def setUp(self): ] ) - def test_no_decay_before_start(self): - current_hours = -10 - decay_hours = 100 - pwr = 1.0 - - alb_v_d, alb_ir_d = decay_alb_power( - self.veg, - self.veg_type, - current_hours, - decay_hours, - pwr, - ALBEDO_VIS, - ALBEDO_IR - ) - - npt.assert_array_equal(ALBEDO_VIS, alb_v_d) - npt.assert_array_equal(ALBEDO_IR, alb_ir_d) - def test_max_decay_after_end(self): """Test with current_hours > decay_hours, ensuring maximum decay is applied.""" current_hours = 150 From 846f83c67df75d39ccb6e490137d8cc348333a7b Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 18 Feb 2026 19:37:33 -0700 Subject: [PATCH 17/30] Envphys - Albedo - Simplify the power decay albedo logic Reduce the complexity of the method by combining common steps between conditional steps. This eliminates all but one if case. --- smrf/envphys/albedo.py | 49 ++++++++++++++---------------------------- 1 file changed, 16 insertions(+), 33 deletions(-) diff --git a/smrf/envphys/albedo.py b/smrf/envphys/albedo.py index 1ad0edb1..1d34411e 100644 --- a/smrf/envphys/albedo.py +++ b/smrf/envphys/albedo.py @@ -1,7 +1,7 @@ import math -import numexpr as ne from typing import Tuple +import numexpr as ne import numpy as np # define some constants @@ -142,37 +142,21 @@ def decay_alb_power( alb_v_d, alb_ir_d : numpy arrays of decayed albedo """ - alb_dec = np.zeros_like(alb_v) - - # Use max decay if after start - if current_hours > decay_hours: - # Use default - alb_dec = alb_dec + veg["default"] - # Decay based on veg type - for k, v in veg.items(): - if isint(k): - alb_dec[veg_type == int(k)] = v - - # Power function decay if during decay period - else: - # Use defaults - max_dec = veg["default"] - tao = decay_hours / (max_dec ** (1.0 / pwr)) - - # Add default decay to array of zeros - alb_dec = alb_dec + (current_hours / tao) ** pwr - - # Decay based on veg type - for k, v in veg.items(): - max_dec = v - tao = decay_hours / (max_dec ** (1.0 / pwr)) - - # Set albedo decay at correct veg types - if isint(k): - alb_dec[veg_type == int(k)] = (current_hours / tao) ** pwr - - alb_v_d = alb_v - alb_dec - alb_ir_d = alb_ir - alb_dec + decay_rates = np.full(veg_type.shape, veg["default"]) + + # Map vegetation-specific decay values to the topo grid + for k, v in veg.items(): + if isint(k): + decay_rates[veg_type == int(k)] = v + + if current_hours < decay_hours: + inv_pwr = 1.0 / pwr + decay_rates = ne.evaluate( + "((current_hours * (decay_rates ** inv_pwr)) / decay_hours) ** pwr" + ) + + alb_v_d = alb_v - decay_rates + alb_ir_d = alb_ir - decay_rates return alb_v_d, alb_ir_d @@ -270,4 +254,3 @@ def decay_burned( alb_ir_d = ne.evaluate("alb_ir * decay_factor") return alb_v_d, alb_ir_d - From 485a0d2709067eca4e508a0e4e5311d7b15c3dd6 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Thu, 19 Feb 2026 08:34:16 -0700 Subject: [PATCH 18/30] Albedo - Add debug logging for chosen method and add error handling for post_fire Added a check that k burned and unburned are set when chossing this method. --- smrf/distribute/albedo.py | 4 +++- smrf/envphys/albedo.py | 6 ++++++ smrf/tests/envphys/test_albedo.py | 30 ++++++++++++++++++++++++++++++ 3 files changed, 39 insertions(+), 1 deletion(-) diff --git a/smrf/distribute/albedo.py b/smrf/distribute/albedo.py index ae8b9480..ff6f698f 100755 --- a/smrf/distribute/albedo.py +++ b/smrf/distribute/albedo.py @@ -162,8 +162,8 @@ def distribute( self.config["dirt"], ) - # Perform litter decay if self.config["decay_method"] == "date_method": + self._logger.debug(" Using date method decay") current_hours, decay_hours = self.decay_window(current_time_step) if current_hours > 0: alb_v, alb_ir = albedo.decay_alb_power( @@ -176,6 +176,7 @@ def distribute( alb_ir, ) elif self.config["decay_method"] == "post_fire": + self._logger.debug(" Using post fire decay") current_hours, decay_hours = self.decay_window(current_time_step) if current_hours > 0: alb_v, alb_ir = albedo.decay_burned( @@ -188,6 +189,7 @@ def distribute( ) elif self.config["decay_method"] == "hardy2000": + self._logger.debug(" Using hardy2000 decay") alb_v, alb_ir = albedo.decay_alb_hardy( self.litter, self.veg_type, storm_day, alb_v, alb_ir ) diff --git a/smrf/envphys/albedo.py b/smrf/envphys/albedo.py index 1d34411e..936444ea 100644 --- a/smrf/envphys/albedo.py +++ b/smrf/envphys/albedo.py @@ -246,6 +246,12 @@ def decay_burned( Returns: alb_v_d, alb_ir_d : numpy arrays of decayed albedo """ + if k_burned is None or k_unburned is None: + raise ValueError( + "Post fire albedo decay is configured, but k_burned and k_unburned are not " + "set in the config file." + ) + decay_factor = ne.evaluate( "exp(-where(burn_mask == 1, k_burned, k_unburned) * last_snow)" ) diff --git a/smrf/tests/envphys/test_albedo.py b/smrf/tests/envphys/test_albedo.py index 9b64b090..8cd4fea6 100644 --- a/smrf/tests/envphys/test_albedo.py +++ b/smrf/tests/envphys/test_albedo.py @@ -166,3 +166,33 @@ def test_decay_burned_all_unburned(self): npt.assert_array_almost_equal( np.array([[0.686139, 0.485223], [0.6, 0.384316]]), alb_ir, decimal=6 ) + + def test_decay_burned_k_burned_none(self): + burn_mask = np.array([[1, 0], [0, 1]]) + + with self.assertRaises(ValueError) as context: + decay_burned( + ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, None, K_UNBURNED + ) + + self.assertIn("k_burned and k_unburned", str(context.exception)) + + def test_decay_burned_k_unburned_none(self): + burn_mask = np.array([[1, 0], [0, 1]]) + + with self.assertRaises(ValueError) as context: + decay_burned( + ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, K_BURNED, None + ) + + self.assertIn("k_burned and k_unburned", str(context.exception)) + + def test_decay_burned_both_k_none(self): + burn_mask = np.array([[1, 0], [0, 1]]) + + with self.assertRaises(ValueError) as context: + decay_burned( + ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, None, None + ) + + self.assertIn("k_burned and k_unburned", str(context.exception)) From 98825849ccc7612a32ccbc0f7c4e0930da185651 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Thu, 19 Feb 2026 11:21:51 -0700 Subject: [PATCH 19/30] Albedo - Post Fire - Only use a different decay for burned pixels. Change of strategy and update the logic to apply standard time-decay albedo for unburned pixels and only a different decay for burned ones. This also uses the initial "clean" albedo after a storm for burned pixels for the burned pixels and then applies the burn decay rate when there is no more new snow. --- smrf/distribute/albedo.py | 114 +++++++++++++++----- smrf/envphys/albedo.py | 115 +++++++++++++------- smrf/framework/CoreConfig.ini | 15 +-- smrf/tests/distribute/test_albedo.py | 151 +++++++++++++++------------ smrf/tests/envphys/test_albedo.py | 99 +++++++++++------- 5 files changed, 315 insertions(+), 179 deletions(-) diff --git a/smrf/distribute/albedo.py b/smrf/distribute/albedo.py index ff6f698f..0967edef 100755 --- a/smrf/distribute/albedo.py +++ b/smrf/distribute/albedo.py @@ -1,6 +1,7 @@ from datetime import datetime from typing import Tuple +import numexpr as ne import numpy as np import numpy.typing as npt @@ -60,6 +61,7 @@ def __init__(self, config: dict, topo): self.albedo_ir = None self.albedo_direct = None self.albedo_diffuse = None + self.burn_mask = None # Get the veg values for the decay methods. Date method uses self.veg # Hardy2000 uses self.litter @@ -84,7 +86,7 @@ def initialize(self, metadata): """ super().initialize(metadata) - self.burn_mask = self.topo.burn_mask + self.load_burn_mask() if ( self.config.get("decay_method", None) is None @@ -92,6 +94,19 @@ def initialize(self, metadata): ): self._logger.warning("No decay method is set!") + def load_burn_mask(self): + """ + Load a post fire burn mask from the topo file. If none is set all values + will be marked as 0. + + Sets: :py:attr:`burn_mask` + """ + burn_mask = getattr(self.topo, "burn_mask", None) + if burn_mask is None: + self.burn_mask = np.zeros_like(self.topo.dem) + else: + self.burn_mask = self.topo.burn_mask + def distribute( self, current_time_step: datetime, cos_z: npt.NDArray, storm_day: npt.NDArray ): @@ -165,27 +180,9 @@ def distribute( if self.config["decay_method"] == "date_method": self._logger.debug(" Using date method decay") current_hours, decay_hours = self.decay_window(current_time_step) - if current_hours > 0: - alb_v, alb_ir = albedo.decay_alb_power( - self.veg, - self.veg_type, - current_hours, - decay_hours, - self.config["date_method_decay_power"], - alb_v, - alb_ir, - ) - elif self.config["decay_method"] == "post_fire": - self._logger.debug(" Using post fire decay") - current_hours, decay_hours = self.decay_window(current_time_step) - if current_hours > 0: - alb_v, alb_ir = albedo.decay_burned( - alb_v, - alb_ir, - storm_day, - self.burn_mask, - self.config["post_fire_k_burned"], - self.config["post_fire_k_unburned"], + if current_hours >= 0: + alb_v, alb_ir = self.date_method( + alb_v, alb_ir, current_hours, decay_hours, storm_day ) elif self.config["decay_method"] == "hardy2000": @@ -201,8 +198,73 @@ def distribute( self.albedo_vis = np.zeros(storm_day.shape) self.albedo_ir = np.zeros(storm_day.shape) + def date_method( + self, + alb_v: npt.NDArray, + alb_ir: npt.NDArray, + current_hours: float, + decay_hours: float, + storm_day: npt.NDArray, + ) -> Tuple[np.ndarray, np.ndarray]: + """ + Apply a power law decay to the albedo and an optional post fire decay if configured. + The post fire decay uses the initial albedo value of the time decay method after + a fresh snowfall. + + :param alb_v: Visibility albedo + :param alb_ir: Infrared albedo + :param current_hours: Time in hours since start of decay window + :param decay_hours: Total hours of the decay window + :param storm_day: Time (days) since last snowfall + + :return: + alb_v, alb_ir : Decayed albedo for visible and infrared spectrum + """ + # Create a mask of burned pixels with no new snowfall + if self.config.get("post_fire", False): + burned_no_snowfall = ne.evaluate( + "(burn_mask == 1) & (storm_day > 0)", + local_dict={"burn_mask": self.burn_mask, "storm_day": storm_day}, + ) + else: + burned_no_snowfall = np.zeros_like(storm_day) + + alb_v, alb_ir = albedo.decay_alb_power( + self.veg, + self.veg_type, + current_hours, + decay_hours, + self.config["date_method_decay_power"], + alb_v, + alb_ir, + burned_no_snowfall, + ) + + if self.config.get("post_fire", False): + self._logger.debug(" Applying post fire decay") + alb_v, alb_ir = albedo.decay_burned( + alb_v, + alb_ir, + storm_day, + self.burn_mask, + self.config.get("post_fire_k_burned", None), + ) + + return alb_v, alb_ir + def decay_window(self, current_timestep: datetime) -> Tuple[float, float]: - # Calculate hour past start of decay + """ + Determine if the current time is within the decay window. + When the window has not been reached, a return value of -1 is returned. + When the window has been reached, the current and total decay hours are returned. + + :param current_timestep: Timestep to check if it is within the decay window + + :return: Tuple(current_hours, decay_hours) + current_hours: Time in hours since start of decay window or -1 if outside the window + decay_hours: Total hours of the decay window + """ + # Calculate current hours past the start of the decay window current_difference = current_timestep - self.config["decay_start"] current_hours = ( current_difference.days * 24.0 + current_difference.seconds / 3600.0 @@ -212,10 +274,8 @@ def decay_window(self, current_timestep: datetime) -> Tuple[float, float]: if current_hours < 0: return -1, 0 - # Calculate total time of decay + # Calculate total time of decay window decay_difference = self.config["decay_end"] - self.config["decay_start"] - decay_hours = ( - decay_difference.days * 24.0 + decay_difference.seconds / 3600.0 - ) + decay_hours = decay_difference.days * 24.0 + decay_difference.seconds / 3600.0 return current_hours, decay_hours diff --git a/smrf/envphys/albedo.py b/smrf/envphys/albedo.py index 936444ea..37ad79fe 100644 --- a/smrf/envphys/albedo.py +++ b/smrf/envphys/albedo.py @@ -3,17 +3,17 @@ import numexpr as ne import numpy as np +import numpy.typing as npt -# define some constants -MAXV = 1.0 # vis albedo when gsize = 0 -MAXIR = 0.85447 # IR albedo when gsize = 0 -IRFAC = -0.02123 # IR decay factor -VFAC = 500.0 # visible decay factor -VZRG = 1.375e-3 # vis zenith increase range factor -IRZRG = 2.0e-3 # ir zenith increase range factor -IRZ0 = 0.1 # ir zenith increase range, gsize=0 -BOIL = 373.15 # boiling temperature K -GRAVITY = 9.80665 # gravity (m/s^2) +MAXV = 1.0 # vis albedo when gsize = 0 +MAXIR = 0.85447 # IR albedo when gsize = 0 +IRFAC = -0.02123 # IR decay factor +VFAC = 500.0 # visible decay factor +VZRG = 1.375e-3 # vis zenith increase range factor +IRZRG = 2.0e-3 # ir zenith increase range factor +IRZ0 = 0.1 # ir zenith increase range, gsize=0 +BOIL = 373.15 # boiling temperature K +GRAVITY = 9.80665 # gravity (m/s^2) # TODO - Need to raise here instead of silently fail to make the user aware @@ -109,13 +109,14 @@ def albedo(telapsed, cosz, gsize, maxgsz, dirt=2) -> Tuple[np.ndarray, np.ndarra def decay_alb_power( veg: dict, - veg_type: np.ndarray, + veg_type: npt.NDArray, current_hours: float, decay_hours: float, pwr: float, - alb_v: np.ndarray, - alb_ir: np.ndarray, -) -> Tuple[np.ndarray, np.ndarray]: + alb_v: npt.NDArray, + alb_ir: npt.NDArray, + burned_no_snowfall: npt.NDArray, +) -> Tuple[npt.NDArray, npt.NDArray]: """ Find a decrease in albedo due to litter accumulation. Decay is based on max decay, decay power, and start and end dates. No litter decay occurs @@ -137,6 +138,7 @@ def decay_alb_power( pwr: power for power law decay alb_v: numpy array of albedo for visible spectrum alb_ir: numpy array of albedo for IR spectrum + burned_no_snowfall: numpy array of burned pixels with no new snowfall in a burn mask Returns: Tuple alb_v_d, alb_ir_d : numpy arrays of decayed albedo @@ -151,14 +153,41 @@ def decay_alb_power( if current_hours < decay_hours: inv_pwr = 1.0 / pwr - decay_rates = ne.evaluate( - "((current_hours * (decay_rates ** inv_pwr)) / decay_hours) ** pwr" + ne.evaluate( + "((current_hours * (decay_rates ** inv_pwr)) / decay_hours) ** pwr", + out=decay_rates, + local_dict={ + "current_hours": current_hours, + "decay_rates": decay_rates, + "decay_hours": decay_hours, + "inv_pwr": inv_pwr, + "pwr": pwr, + }, ) - alb_v_d = alb_v - decay_rates - alb_ir_d = alb_ir - decay_rates + # Only apply the standard decay rate if it is not in a pixel that was burned and + # no new snow was deposited. Burned pixels still get the initial albedo after a storm + # in the decay window from the power law + ne.evaluate( + "where(burned_no_snowfall == 0, alb_v - decay_rates, alb_v)", + out=alb_v, + local_dict={ + "alb_v": alb_v, + "decay_rates": decay_rates, + "burned_no_snowfall": burned_no_snowfall, + }, + ) + ne.evaluate( + "where(burned_no_snowfall == 0, alb_ir - decay_rates, alb_ir)", + out=alb_ir, + local_dict={ + "alb_ir": alb_ir, + "decay_rates": decay_rates, + "burned_no_snowfall": burned_no_snowfall, + }, + ) - return alb_v_d, alb_ir_d + return alb_v, alb_ir def decay_alb_hardy(litter, veg_type, storm_day, alb_v, alb_ir): @@ -201,39 +230,37 @@ def decay_alb_hardy(litter, veg_type, storm_day, alb_v, alb_ir): # array for decimal percent snow coverage sc = np.zeros_like(alb_v) # calculate snow coverage default veg type - l_rate = litter['default'] - alb_litter = litter['albedo'] + l_rate = litter["default"] + alb_litter = litter["albedo"] - sc = sc + (1.0-l_rate)**(storm_day) + sc = sc + (1.0 - l_rate) ** (storm_day) # calculate snow coverage based on veg type for k, v in litter.items(): - l_rate = litter[k] if isint(k): - sc[veg_type == int(k)] = ( - 1.0 - l_rate)**(storm_day[veg_type == int(k)]) + sc[veg_type == int(k)] = (1.0 - l_rate) ** (storm_day[veg_type == int(k)]) # calculate litter coverage lc = np.ones_like(alb_v) - sc # weighted average to find decayed albedo - alb_v_d = alb_v*sc + alb_litter*lc - alb_ir_d = alb_ir*sc + alb_litter*lc + alb_v_d = alb_v * sc + alb_litter * lc + alb_ir_d = alb_ir * sc + alb_litter * lc return alb_v_d, alb_ir_d + def decay_burned( alb_v: np.ndarray, alb_ir: np.ndarray, last_snow: np.ndarray, burn_mask: np.ndarray, k_burned: float, - k_unburned: float, ) -> Tuple[np.ndarray, np.ndarray]: """ Apply exponential albedo decay as a function of time since last snowfall. - Different decay rates are applied for burned and unburned pixels. + This changes all pixels which are marked via the burn mask in the topo file. Args: alb_v: Visible albedo @@ -241,22 +268,34 @@ def decay_burned( last_snow: Time since last snow storm (decimal days) burn_mask: Mask of burned area k_burned: Decay rate for burned area - k_unburned: Decay rate unburned area - Returns: + Returns: Tuple alb_v_d, alb_ir_d : numpy arrays of decayed albedo """ - if k_burned is None or k_unburned is None: + if k_burned is None: raise ValueError( - "Post fire albedo decay is configured, but k_burned and k_unburned are not " - "set in the config file." + "Post fire albedo decay is configured, but post_fire_k_burned is not configured " + "in the config file." ) decay_factor = ne.evaluate( - "exp(-where(burn_mask == 1, k_burned, k_unburned) * last_snow)" + "exp(-where(burn_mask == 1, k_burned, 0) * last_snow)", + local_dict={ + "burn_mask": burn_mask, + "k_burned": k_burned, + "last_snow": last_snow, + }, ) - alb_v_d = ne.evaluate("alb_v * decay_factor") - alb_ir_d = ne.evaluate("alb_ir * decay_factor") + ne.evaluate( + "alb_v * decay_factor", + out=alb_v, + local_dict={"alb_v": alb_v, "decay_factor": decay_factor}, + ) + ne.evaluate( + "alb_ir * decay_factor", + out=alb_ir, + local_dict={"alb_ir": alb_ir, "decay_factor": decay_factor}, + ) - return alb_v_d, alb_ir_d + return alb_v, alb_ir diff --git a/smrf/framework/CoreConfig.ini b/smrf/framework/CoreConfig.ini index d483ab08..75d120cc 100755 --- a/smrf/framework/CoreConfig.ini +++ b/smrf/framework/CoreConfig.ini @@ -791,7 +791,7 @@ description = Effective contamination for adjustment to visible albedo (usually decay_method : default = None, -options = [hardy2000 date_method post_fire None], +options = [hardy2000 date_method None], description = Describe how the albedo decays in the late season decay_start : @@ -864,18 +864,19 @@ description = Full path to netcdf files that contain albedo values. Once set none of the other albedo variables will be effective and the raw values from the files will be used. +post_fire : +default = False, +type = bool, +description = Apply a different post fire decay method for areas masked in the topo + file via a 'burn_mask' layer. This also requires to set the + 'post_fire_k_burned' value in this configuration section + post_fire_k_burned: default = None, type = float, description = Decay power for areas masked as 1 (burned) in the burned mask of the topo file -post_fire_k_unburned: -default = None, -type = float, -description = Decay power for areas masked as 0 (unburned) in the burned mask of the - topo file - ################################################################################ # Cloud Factor - Fraction used to limit solar radiation Cloudy (0) - Sunny (1) ################################################################################ diff --git a/smrf/tests/distribute/test_albedo.py b/smrf/tests/distribute/test_albedo.py index 31926746..d9f36728 100644 --- a/smrf/tests/distribute/test_albedo.py +++ b/smrf/tests/distribute/test_albedo.py @@ -41,8 +41,8 @@ DATA = MagicMock() -ALBEDO_VIS = np.array([[0.0, 1.0], [1.0, 0.0]]) -ALBEDO_IR = np.array([[0.0, 1.0], [1.0, 0.0]]) +ALBEDO_VIS = np.array([[0.9, 1.0], [1.0, 0.8]]) +ALBEDO_IR = np.array([[0.8, 1.0], [1.0, 0.9]]) class TestAlbedo(SMRFConfig, unittest.TestCase): @@ -64,6 +64,22 @@ def test_init_default(self): self.assertEqual(TIMESTEP.strftime("%Y%m%d"), subject.start_date) self.assertEqual(TIMEZONE, subject.time_zone) + def test_load_burn_mask_with_defined_mask(self): + self.subject.topo.burn_mask = np.array([[1, 0], [0, 1]]) + + self.subject.load_burn_mask() + + np.testing.assert_array_equal(self.subject.burn_mask, self.subject.topo.burn_mask) + + def test_load_burn_mask_without_defined_mask(self): + self.subject.topo.burn_mask = None + self.subject.topo.dem = np.array([[1, 0], [0, 1]]) + + self.subject.load_burn_mask() + + expected_burn_mask = np.zeros_like(self.subject.topo.dem) + np.testing.assert_array_equal(self.subject.burn_mask, expected_burn_mask) + @patch("smrf.distribute.variable_base.ReadNetCDF") def test_distribute_file_broadband(self, mock_read_netcdf): values = np.array([0.9, 0.5]) @@ -77,7 +93,7 @@ def test_distribute_file_broadband(self, mock_read_netcdf): subject = Albedo(config, TOPO) subject.initialize(pd.DataFrame()) - subject.distribute(TIMESTEP, np.ndarray([10]), STORM_DAYS) + subject.distribute(TIMESTEP, COS_Z, STORM_DAYS) npt.assert_equal(subject.albedo, values) mock_source.load.assert_called_with("albedo", TIMESTEP) @@ -96,7 +112,7 @@ def test_distribute_file_vis_ir(self, mock_read_netcdf): subject = Albedo(config, TOPO) subject.initialize(pd.DataFrame()) - subject.distribute(TIMESTEP, np.ndarray([10]), STORM_DAYS) + subject.distribute(TIMESTEP, COS_Z, STORM_DAYS) np.testing.assert_array_equal(subject.albedo_vis, values_vis) np.testing.assert_array_equal(subject.albedo_ir, values_ir) @@ -120,7 +136,7 @@ def test_distribute_file_direct_diffuse(self, mock_read_netcdf): subject = Albedo(config, TOPO) subject.initialize(pd.DataFrame()) - subject.distribute(TIMESTEP, np.ndarray([10]), STORM_DAYS) + subject.distribute(TIMESTEP, COS_Z, STORM_DAYS) np.testing.assert_array_equal(subject.albedo_direct, values_direct) np.testing.assert_array_equal(subject.albedo_diffuse, values_diffuse) @@ -130,31 +146,21 @@ def test_distribute_file_direct_diffuse(self, mock_read_netcdf): mock_source.load.assert_any_call("albedo_direct", TIMESTEP) mock_source.load.assert_any_call("albedo_diffuse", TIMESTEP) - @patch("smrf.envphys.albedo.decay_alb_power") + @patch("smrf.distribute.albedo.Albedo.date_method") @patch("smrf.envphys.albedo.albedo") - def test_distribute_calculated(self, mock_albedo_calc, mock_albedo_decay): - values_vis = np.array([0.8, 0.7]) - values_ir = np.array([0.9, 0.5]) - mock_albedo_calc.return_value = (values_vis, values_ir) - - values_vis_decay = np.array([0.75, 0.65]) - values_ir_decay = np.array([0.85, 0.45]) - mock_albedo_decay.return_value = (values_vis_decay, values_ir_decay) + def test_distribute_calculated(self, mock_albedo_calc, mock_date_method): + mock_albedo_calc.return_value = (ALBEDO_VIS, ALBEDO_IR) + date_albedo_vis = ALBEDO_VIS - 0.1 + date_albedo_ir = ALBEDO_IR - 0.2 + mock_date_method.return_value = (date_albedo_vis, date_albedo_ir) - config = self._copy_config(CONFIG) - config["albedo"]["grain_size"] = "100" - config["albedo"]["max_grain"] = "700" - config["albedo"]["date_method_decay_power"] = "1.2" - - subject = Albedo(config, TOPO) - subject.initialize(pd.DataFrame()) - - subject.distribute(TIMESTEP, np.ndarray([10]), STORM_DAYS) + self.subject.distribute(TIMESTEP, COS_Z, STORM_DAYS) mock_albedo_calc.assert_called_once() - mock_albedo_decay.assert_called_once() - np.testing.assert_array_equal(subject.albedo_vis, values_vis_decay) - np.testing.assert_array_equal(subject.albedo_ir, values_ir_decay) + mock_date_method.assert_called_once() + + npt.assert_array_equal(self.subject.albedo_vis, date_albedo_vis) + npt.assert_array_equal(self.subject.albedo_ir, date_albedo_ir) def test_before_decay_window(self): current, decay = self.subject.decay_window(NO_DECAY_TIME) @@ -166,9 +172,9 @@ def test_in_decay_window(self): self.assertEqual(2184.0, decay) @patch("smrf.distribute.albedo.albedo.albedo", return_value=(ALBEDO_VIS, ALBEDO_IR)) - @patch("smrf.distribute.albedo.albedo.decay_alb_power") + @patch("smrf.distribute.albedo.Albedo.date_method") def test_distribute_date_method_not_in_window( - self, decay_alb_power, envphys_albedo + self, date_decay_method, envphys_albedo ): self.subject.distribute(NO_DECAY_TIME, COS_Z, STORM_DAYS) @@ -180,49 +186,60 @@ def test_distribute_date_method_not_in_window( CONFIG["albedo"]["max_grain"], CONFIG["albedo"]["dirt"], ) - decay_alb_power.assert_not_called() + date_decay_method.assert_not_called() - @patch("smrf.distribute.albedo.albedo.albedo", return_value=(ALBEDO_VIS, ALBEDO_IR)) - @patch( - "smrf.distribute.albedo.albedo.decay_alb_power", - return_value=(ALBEDO_VIS, ALBEDO_IR), - ) - def test_distribute_date_method_in_window(self, decay_alb_power, _envphys_albedo): - self.subject.config["date_method_decay_power"] = 0.7 - self.subject.config["date_method_veg_default"] = 0.2 - self.subject.distribute(DECAY_TIME, COS_Z, STORM_DAYS) - - decay_alb_power.assert_called_with( - {"default": 0.2}, TOPO.veg_type, 696.0, 2184.0, 0.7, ALBEDO_VIS, ALBEDO_IR + @patch("smrf.envphys.albedo.decay_alb_power") + def test_date_method_no_post_fire(self, mock_decay_power): + expected_v = ALBEDO_VIS - 0.05 + expected_ir = ALBEDO_IR - 0.05 + mock_decay_power.return_value = (expected_v, expected_ir) + + self.subject.config["post_fire"] = False + current_hours, decay_hours = self.subject.decay_window(DECAY_TIME) + + res_v, res_ir = self.subject.date_method( + ALBEDO_VIS, ALBEDO_IR, current_hours, decay_hours, STORM_DAYS ) - @patch("smrf.distribute.albedo.albedo.albedo", return_value=(ALBEDO_VIS, ALBEDO_IR)) - @patch("smrf.distribute.albedo.albedo.decay_burned") - def test_distribute_post_fire_not_in_window(self, decay_burned, envphys_albedo): - self.subject.config["decay_method"] = "post_fire" - self.subject.distribute(NO_DECAY_TIME, COS_Z, STORM_DAYS) + # Burned_no_snowfall mask should be all zeros (last argument) + burn_mask = mock_decay_power.call_args[0][-1] + npt.assert_array_equal(burn_mask, np.zeros_like(STORM_DAYS)) - envphys_albedo.assert_called() - envphys_albedo.assert_called_with( - STORM_DAYS, - COS_Z, - CONFIG["albedo"]["grain_size"], - CONFIG["albedo"]["max_grain"], - CONFIG["albedo"]["dirt"], + npt.assert_array_equal(res_v, expected_v) + npt.assert_array_equal(res_ir, expected_ir) + + @patch("smrf.envphys.albedo.decay_burned") + @patch("smrf.envphys.albedo.decay_alb_power") + def test_date_method_with_post_fire(self, mock_decay_power, mock_decay_burned): + self.subject.config["post_fire"] = True + self.subject.config["post_fire_k_burned"] = 1.2 + burn_mask = np.array([[1.0, 1.0], [0.0, 0.0]]) + self.subject.burn_mask = burn_mask + + power_v = ALBEDO_VIS - 0.02 + power_ir = ALBEDO_IR - 0.02 + mock_decay_power.return_value = (power_v, power_ir) + + final_v = power_v - 0.1 + final_ir = power_ir - 0.1 + mock_decay_burned.return_value = (final_v, final_ir) + + current_hours, decay_hours = self.subject.decay_window(DECAY_TIME) + + res_v, res_ir = self.subject.date_method( + ALBEDO_VIS, ALBEDO_IR, current_hours, decay_hours, STORM_DAYS ) - decay_burned.assert_not_called() - @patch("smrf.distribute.albedo.albedo.albedo", return_value=(ALBEDO_VIS, ALBEDO_IR)) - @patch( - "smrf.distribute.albedo.albedo.decay_burned", - return_value=(ALBEDO_VIS, ALBEDO_IR), - ) - def test_distribute_post_fire_in_window(self, decay_burned, _envphys_albedo): - self.subject.config["decay_method"] = "post_fire" - self.subject.config["post_fire_k_burned"] = 0.06 - self.subject.config["post_fire_k_unburned"] = 0.02 - self.subject.distribute(DECAY_TIME, COS_Z, STORM_DAYS) - - decay_burned.assert_called_with( - ALBEDO_VIS, ALBEDO_IR, STORM_DAYS, TOPO.burn_mask, 0.06, 0.02 + expected_mask = np.array([[False, True], [False, False]]) + received_mask = mock_decay_power.call_args[0][-1] + npt.assert_array_equal(received_mask, expected_mask) + + mock_decay_burned.assert_called_once_with( + power_v, + power_ir, + STORM_DAYS, + self.subject.burn_mask, + self.subject.config["post_fire_k_burned"], ) + + npt.assert_array_equal(res_v, final_v) diff --git a/smrf/tests/envphys/test_albedo.py b/smrf/tests/envphys/test_albedo.py index 8cd4fea6..c8c547a9 100644 --- a/smrf/tests/envphys/test_albedo.py +++ b/smrf/tests/envphys/test_albedo.py @@ -8,8 +8,8 @@ ALBEDO_VIS = np.array([[0.8, 0.6], [0.7, 0.5]]) ALBEDO_IR = np.array([[0.7, 0.5], [0.6, 0.4]]) LAST_SNOW = np.array([[1.0, 1.5], [0.0, 2.0]]) +BURNED_NO_SNOWFALL = np.array([[0.0, 0.0], [0.0, 0.0]]) K_BURNED = 0.06 -K_UNBURNED = 0.02 class TestDecayAlbPower(unittest.TestCase): @@ -51,8 +51,9 @@ def test_max_decay_after_end(self): current_hours, decay_hours, pwr, - ALBEDO_VIS, - ALBEDO_IR + ALBEDO_VIS.copy(), + ALBEDO_IR.copy(), + BURNED_NO_SNOWFALL, ) npt.assert_array_almost_equal(expected_v, alb_v_d) @@ -83,8 +84,46 @@ def test_power_decay_during_window(self): current_hours, decay_hours, pwr, - ALBEDO_VIS, - ALBEDO_IR + ALBEDO_VIS.copy(), + ALBEDO_IR.copy(), + BURNED_NO_SNOWFALL, + ) + + npt.assert_array_almost_equal(expected_v, alb_v_d, decimal=5) + npt.assert_array_almost_equal(expected_ir, alb_ir_d, decimal=5) + + def test_power_decay_during_window_with_burn(self): + """Burned pixels will not be changed by power decay after snowfall""" + current_hours = 50 + decay_hours = 100 + pwr = 2.0 + burned_no_snowfall = np.array([[0.0, 1.0], [0.0, 0.0]]) + + # Decay rates + tao_default = decay_hours / (self.veg["default"] ** (1.0 / pwr)) + tao_41 = decay_hours / (self.veg["41"] ** (1.0 / pwr)) + tao_42 = decay_hours / (self.veg["42"] ** (1.0 / pwr)) + + expected_decay = np.zeros_like(self.veg_type, dtype=float) + expected_decay[self.veg_type == 0] = (current_hours / tao_default) ** pwr + expected_decay[self.veg_type == 41] = (current_hours / tao_41) ** pwr + expected_decay[self.veg_type == 42] = (current_hours / tao_42) ** pwr + + # Burned pixels should not be changed by power decay + expected_v = ALBEDO_VIS - expected_decay + expected_v[0][1] = ALBEDO_VIS[0][1] + expected_ir = ALBEDO_IR - expected_decay + expected_ir[0][1] = ALBEDO_IR[0][1] + + alb_v_d, alb_ir_d = decay_alb_power( + self.veg, + self.veg_type, + current_hours, + decay_hours, + pwr, + ALBEDO_VIS.copy(), + ALBEDO_IR.copy(), + burned_no_snowfall, ) npt.assert_array_almost_equal(expected_v, alb_v_d, decimal=5) @@ -117,8 +156,9 @@ def test_partial_coverage_during_decay(self): current_hours, decay_hours, pwr, - ALBEDO_VIS, - ALBEDO_IR + ALBEDO_VIS.copy(), + ALBEDO_IR.copy(), + BURNED_NO_SNOWFALL, ) npt.assert_array_almost_equal(expected_v, alb_v_d, decimal=5) @@ -128,22 +168,24 @@ class TestDecayBurned(unittest.TestCase): def test_decay_burned_with_burned_and_unburned_areas(self): burn_mask = np.array([[1, 0], [0, 1]]) + # Arrays are updated in place, so pass a copy alb_vis, alb_ir = decay_burned( - ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, K_BURNED, K_UNBURNED + ALBEDO_VIS.copy(), ALBEDO_IR.copy(), LAST_SNOW, burn_mask, K_BURNED ) npt.assert_array_almost_equal( - np.array([[0.753412, 0.582267], [0.7, 0.44346]]), alb_vis, decimal=6 + np.array([[0.753412, ALBEDO_VIS[0][1]], [ALBEDO_VIS[1][0], 0.44346]]), alb_vis, decimal=6 ) npt.assert_array_almost_equal( - np.array([[0.659235, 0.485223], [0.6, 0.354768]]), alb_ir, decimal=6 + np.array([[0.659235, ALBEDO_IR[0][1]], [ALBEDO_IR[1][0], 0.354768]]), alb_ir, decimal=6 ) def test_decay_burned_all_burned(self): burn_mask = np.ones_like(ALBEDO_VIS) + # Arrays are updated in place, so pass a copy alb_vis, alb_ir = decay_burned( - ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, K_BURNED, K_UNBURNED + ALBEDO_VIS.copy(), ALBEDO_IR.copy(), LAST_SNOW, burn_mask, K_BURNED ) npt.assert_array_almost_equal( @@ -156,43 +198,20 @@ def test_decay_burned_all_burned(self): def test_decay_burned_all_unburned(self): burn_mask = np.zeros_like(ALBEDO_VIS) + # Arrays are updated in place, so pass a copy alb_vis, alb_ir = decay_burned( - ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, K_BURNED, K_UNBURNED + ALBEDO_VIS.copy(), ALBEDO_IR.copy(), LAST_SNOW, burn_mask, K_BURNED ) - npt.assert_array_almost_equal( - np.array([[0.784159, 0.582267], [0.7, 0.480395]]), alb_vis, decimal=6 - ) - npt.assert_array_almost_equal( - np.array([[0.686139, 0.485223], [0.6, 0.384316]]), alb_ir, decimal=6 - ) + npt.assert_array_equal(ALBEDO_VIS, alb_vis) + npt.assert_array_equal(ALBEDO_IR, alb_ir) def test_decay_burned_k_burned_none(self): burn_mask = np.array([[1, 0], [0, 1]]) with self.assertRaises(ValueError) as context: decay_burned( - ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, None, K_UNBURNED - ) - - self.assertIn("k_burned and k_unburned", str(context.exception)) - - def test_decay_burned_k_unburned_none(self): - burn_mask = np.array([[1, 0], [0, 1]]) - - with self.assertRaises(ValueError) as context: - decay_burned( - ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, K_BURNED, None - ) - - self.assertIn("k_burned and k_unburned", str(context.exception)) - - def test_decay_burned_both_k_none(self): - burn_mask = np.array([[1, 0], [0, 1]]) - - with self.assertRaises(ValueError) as context: - decay_burned( - ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, None, None + ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, None ) - self.assertIn("k_burned and k_unburned", str(context.exception)) + self.assertIn("k_burned", str(context.exception)) From 52070d06c4af87fcbb8b7fad7cdbf05d786d8b02 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Mon, 23 Feb 2026 08:32:51 -0700 Subject: [PATCH 20/30] Envphys - Albedo - Ensure that all constants are float32 Numexpr will throw an error when trying to write float64 results back into a float32. Also partially addresses issue in #59 --- smrf/envphys/albedo.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/smrf/envphys/albedo.py b/smrf/envphys/albedo.py index 37ad79fe..a1f374ff 100644 --- a/smrf/envphys/albedo.py +++ b/smrf/envphys/albedo.py @@ -144,7 +144,7 @@ def decay_alb_power( alb_v_d, alb_ir_d : numpy arrays of decayed albedo """ - decay_rates = np.full(veg_type.shape, veg["default"]) + decay_rates = np.full(veg_type.shape, veg["default"]).astype(np.float32, order="C") # Map vegetation-specific decay values to the topo grid for k, v in veg.items(): @@ -157,11 +157,11 @@ def decay_alb_power( "((current_hours * (decay_rates ** inv_pwr)) / decay_hours) ** pwr", out=decay_rates, local_dict={ - "current_hours": current_hours, + "current_hours": np.float32(current_hours), "decay_rates": decay_rates, - "decay_hours": decay_hours, - "inv_pwr": inv_pwr, - "pwr": pwr, + "decay_hours": np.float32(decay_hours), + "inv_pwr": np.float32(inv_pwr), + "pwr": np.float32(pwr), }, ) @@ -282,7 +282,7 @@ def decay_burned( "exp(-where(burn_mask == 1, k_burned, 0) * last_snow)", local_dict={ "burn_mask": burn_mask, - "k_burned": k_burned, + "k_burned": np.float32(k_burned), "last_snow": last_snow, }, ) From 460b21aecb20425d8ef046525f09c9c50f3fb42c Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Mon, 23 Feb 2026 08:33:31 -0700 Subject: [PATCH 21/30] Albedo - Update argument and return types using numpy's typing module --- smrf/distribute/albedo.py | 2 +- smrf/envphys/albedo.py | 14 ++++++++------ 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/smrf/distribute/albedo.py b/smrf/distribute/albedo.py index 0967edef..84724070 100755 --- a/smrf/distribute/albedo.py +++ b/smrf/distribute/albedo.py @@ -205,7 +205,7 @@ def date_method( current_hours: float, decay_hours: float, storm_day: npt.NDArray, - ) -> Tuple[np.ndarray, np.ndarray]: + ) -> Tuple[npt.NDArray, npt.NDArray]: """ Apply a power law decay to the albedo and an optional post fire decay if configured. The post fire decay uses the initial albedo value of the time decay method after diff --git a/smrf/envphys/albedo.py b/smrf/envphys/albedo.py index a1f374ff..47a6ca2b 100644 --- a/smrf/envphys/albedo.py +++ b/smrf/envphys/albedo.py @@ -40,7 +40,9 @@ def growth(t): return 1.0 - factor -def albedo(telapsed, cosz, gsize, maxgsz, dirt=2) -> Tuple[np.ndarray, np.ndarray]: +def albedo( + telapsed: npt.NDArray, cosz: npt.NDArray, gsize: float, maxgsz: float, dirt=2 +) -> Tuple[npt.NDArray, npt.NDArray]: """ Calculate the albedo, adapted from IPW function albedo @@ -252,12 +254,12 @@ def decay_alb_hardy(litter, veg_type, storm_day, alb_v, alb_ir): def decay_burned( - alb_v: np.ndarray, - alb_ir: np.ndarray, - last_snow: np.ndarray, - burn_mask: np.ndarray, + alb_v: npt.NDArray, + alb_ir: npt.NDArray, + last_snow: npt.NDArray, + burn_mask: npt.NDArray, k_burned: float, -) -> Tuple[np.ndarray, np.ndarray]: +) -> Tuple[npt.NDArray, npt.NDArray]: """ Apply exponential albedo decay as a function of time since last snowfall. This changes all pixels which are marked via the burn mask in the topo file. From d5dfbfc0cf818124e96b3ea142202937fc93cadf Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Mon, 16 Mar 2026 11:16:45 -0600 Subject: [PATCH 22/30] Tests - Topo - Fix faulty merge with main --- smrf/tests/data/test_topo.py | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) diff --git a/smrf/tests/data/test_topo.py b/smrf/tests/data/test_topo.py index 8eed6f11..ff6f8380 100644 --- a/smrf/tests/data/test_topo.py +++ b/smrf/tests/data/test_topo.py @@ -29,13 +29,6 @@ def setUp(cls): def tearDown(cls): cls.ds.close() -<<<<<<< HEAD - @mock.patch.object(Topo, "gradient") - @mock.patch.object(Topo, "readNetCDF") - def test_init(self, read_nc, gradient): - topo = Topo(TOPO_CONFIG) - self.assertEqual(TOPO_CONFIG, topo.topoConfig) -======= def test_read_topo_images(self): self.assertEqual( [ @@ -50,7 +43,6 @@ def test_read_topo_images(self): Topo.IMAGES ) - def test_topo_image_as_variables(self): for variable in Topo.IMAGES: self.assertTrue(hasattr(self.topo, variable)) @@ -60,10 +52,11 @@ def test_no_burn_mask_required(self): @mock.patch.object(Topo, 'readNetCDF') @mock.patch.object(Topo, 'gradient') def test_init(self, gradient, read_nc): - topo = Topo(self.TOPO_CONFIG) - self.assertEqual(self.TOPO_CONFIG, topo.topoConfig) + # Need to create a local topo instance to test the mock patch calls + topo = Topo(TOPO_CONFIG) + self.assertEqual(TOPO_CONFIG, topo.topoConfig) gradient.assert_called_once() ->>>>>>> 190ae1b (Topo - Read in a burn mask if present.) + read_nc.assert_called_once() gradient.assert_called_once() From 440d12dc99ac1492d7d11f2c9dd8dddad63b1b12 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Mon, 16 Mar 2026 11:38:52 -0600 Subject: [PATCH 23/30] Envphys - albedo - burn decay - Start decay once one day has passed With the `exp` function using the last snow as a multiplier, any period less than a day caused an increase an albedo before. To prevent these short-term spikes, the burn decay is applied after one full day sine a storm. --- smrf/distribute/albedo.py | 8 ++- smrf/envphys/albedo.py | 22 ++++--- smrf/tests/distribute/test_albedo.py | 99 ++++++++++++++++++++++++++-- smrf/tests/envphys/test_albedo.py | 13 ++-- 4 files changed, 121 insertions(+), 21 deletions(-) diff --git a/smrf/distribute/albedo.py b/smrf/distribute/albedo.py index 84724070..1361bd7e 100755 --- a/smrf/distribute/albedo.py +++ b/smrf/distribute/albedo.py @@ -220,10 +220,12 @@ def date_method( :return: alb_v, alb_ir : Decayed albedo for visible and infrared spectrum """ - # Create a mask of burned pixels with no new snowfall + # Create a mask of burned pixels with snow more than one day ago. + # Only burned pixels that have at least one day since the last snow will + # have a different decay rate applied. See :py:meth:`albedo.decay_burned` if self.config.get("post_fire", False): burned_no_snowfall = ne.evaluate( - "(burn_mask == 1) & (storm_day > 0)", + "(burn_mask == 1) & (storm_day >= 1)", local_dict={"burn_mask": self.burn_mask, "storm_day": storm_day}, ) else: @@ -246,7 +248,7 @@ def date_method( alb_v, alb_ir, storm_day, - self.burn_mask, + burned_no_snowfall, self.config.get("post_fire_k_burned", None), ) diff --git a/smrf/envphys/albedo.py b/smrf/envphys/albedo.py index 47a6ca2b..6fb864cb 100644 --- a/smrf/envphys/albedo.py +++ b/smrf/envphys/albedo.py @@ -140,7 +140,8 @@ def decay_alb_power( pwr: power for power law decay alb_v: numpy array of albedo for visible spectrum alb_ir: numpy array of albedo for IR spectrum - burned_no_snowfall: numpy array of burned pixels with no new snowfall in a burn mask + burned_no_snowfall: numpy array of burned pixels with at least one day since + last snowfall Returns: Tuple alb_v_d, alb_ir_d : numpy arrays of decayed albedo @@ -257,19 +258,22 @@ def decay_burned( alb_v: npt.NDArray, alb_ir: npt.NDArray, last_snow: npt.NDArray, - burn_mask: npt.NDArray, + burned_no_snowfall: npt.NDArray, k_burned: float, ) -> Tuple[npt.NDArray, npt.NDArray]: """ Apply exponential albedo decay as a function of time since last snowfall. - This changes all pixels which are marked via the burn mask in the topo file. + This changes all pixels which are marked via the burn mask in the topo file and + had at least one day since it snowed. The one day minimum is applied to account + for the `exp()` function, where less than 1 values would cause a growth instead + of a decline. Args: - alb_v: Visible albedo - alb_ir: Infrared albedo - last_snow: Time since last snow storm (decimal days) - burn_mask: Mask of burned area - k_burned: Decay rate for burned area + alb_v: Visible albedo + alb_ir: Infrared albedo + last_snow: Time since last snow storm (decimal days) + burned_no_snowfall: Mask of burned area with at least one day since snowfall + k_burned: Decay rate for burned area Returns: Tuple alb_v_d, alb_ir_d : numpy arrays of decayed albedo @@ -283,7 +287,7 @@ def decay_burned( decay_factor = ne.evaluate( "exp(-where(burn_mask == 1, k_burned, 0) * last_snow)", local_dict={ - "burn_mask": burn_mask, + "burn_mask": burned_no_snowfall, "k_burned": np.float32(k_burned), "last_snow": last_snow, }, diff --git a/smrf/tests/distribute/test_albedo.py b/smrf/tests/distribute/test_albedo.py index d9f36728..6c2b03ec 100644 --- a/smrf/tests/distribute/test_albedo.py +++ b/smrf/tests/distribute/test_albedo.py @@ -1,5 +1,5 @@ import unittest -from unittest.mock import MagicMock, patch +from unittest.mock import ANY, MagicMock, patch import numpy as np import numpy.testing as npt @@ -7,10 +7,11 @@ import pytz from smrf.distribute.albedo import Albedo +from smrf.envphys.albedo import decay_alb_power, decay_burned from smrf.tests.smrf_config import SMRFConfig TOPO = MagicMock( - veg_type=MagicMock("veg_type"), + veg_type=np.array([[1.0, 1.0], [1.0, 1.0]]), ) CONFIG = { "time": { @@ -69,7 +70,9 @@ def test_load_burn_mask_with_defined_mask(self): self.subject.load_burn_mask() - np.testing.assert_array_equal(self.subject.burn_mask, self.subject.topo.burn_mask) + np.testing.assert_array_equal( + self.subject.burn_mask, self.subject.topo.burn_mask + ) def test_load_burn_mask_without_defined_mask(self): self.subject.topo.burn_mask = None @@ -238,8 +241,96 @@ def test_date_method_with_post_fire(self, mock_decay_power, mock_decay_burned): power_v, power_ir, STORM_DAYS, - self.subject.burn_mask, + ANY, self.subject.config["post_fire_k_burned"], ) + npt.assert_array_equal(mock_decay_burned.call_args[0][3], expected_mask) npt.assert_array_equal(res_v, final_v) + + @patch("smrf.envphys.albedo.decay_burned") + @patch("smrf.envphys.albedo.decay_alb_power") + def test_date_method_mask_with_post_fire_post_storm( + self, mock_decay_power, mock_decay_burned + ): + """ + Test the window where a snow storm was within the last 24 hrs. There, the burn + decay is noe applied. + """ + self.subject.config["post_fire"] = True + self.subject.config["post_fire_k_burned"] = 1.2 + burn_mask = np.array([[1.0, 1.0], [0.0, 0.0]]) + self.subject.burn_mask = burn_mask + storm_days = np.array([[0.3, 2.0], [1.0, 0.0]]) + + # Preserve the method return signature of mocked methods + power_v = 1 + power_ir = 1 + mock_decay_power.return_value = (power_v, power_ir) + mock_decay_burned.return_value = (1, 1) + + current_hours, decay_hours = self.subject.decay_window(DECAY_TIME) + + self.subject.date_method( + ALBEDO_VIS, ALBEDO_IR, current_hours, decay_hours, storm_days + ) + + expected_mask = np.array([[False, True], [False, False]]) + npt.assert_array_equal(mock_decay_power.call_args[0][-1], expected_mask) + + mock_decay_burned.assert_called_once_with( + power_v, + power_ir, + storm_days, + ANY, + self.subject.config["post_fire_k_burned"], + ) + npt.assert_array_equal(mock_decay_burned.call_args[0][3], expected_mask) + + def test_date_method_with_post_fire_post_storm(self): + """ + Integration test that uses both decay methods + """ + self.subject.config["post_fire"] = True + self.subject.config["post_fire_k_burned"] = 1.2 + self.subject.burn_mask = np.array([[1.0, 1.0], [1.0, 0.0]]) + storm_days = np.array([[0.3, 2.0], [1.0, 0.0]]) + + current_hours, decay_hours = self.subject.decay_window(DECAY_TIME) + expected_mask = np.array([[False, True], [True, False]]) + power_vis, power_ir = decay_alb_power( + self.subject.veg, + self.subject.veg_type, + current_hours, + decay_hours, + self.subject.config["date_method_decay_power"], + ALBEDO_VIS.copy(), + ALBEDO_IR.copy(), + expected_mask + ) + burn_vis, burn_ir = decay_burned( + power_vis.copy(), + power_ir.copy(), + storm_days, + expected_mask, + self.subject.config["post_fire_k_burned"], + ) + + albedo_v, albedo_ir = self.subject.date_method( + ALBEDO_VIS.copy(), ALBEDO_IR.copy(), current_hours, decay_hours, storm_days + ) + + npt.assert_array_equal( + albedo_v, + np.array([ + [power_vis[0][0], burn_vis[0][1]], + [burn_vis[1][0], power_vis[1][1]] + ]), + ) + npt.assert_array_equal( + albedo_ir, + np.array([ + [power_ir[0][0], burn_ir[0][1]], + [burn_ir[1][0], power_ir[1][1]] + ]), + ) diff --git a/smrf/tests/envphys/test_albedo.py b/smrf/tests/envphys/test_albedo.py index c8c547a9..3c764856 100644 --- a/smrf/tests/envphys/test_albedo.py +++ b/smrf/tests/envphys/test_albedo.py @@ -164,6 +164,7 @@ def test_partial_coverage_during_decay(self): npt.assert_array_almost_equal(expected_v, alb_v_d, decimal=5) npt.assert_array_almost_equal(expected_ir, alb_ir_d, decimal=5) + class TestDecayBurned(unittest.TestCase): def test_decay_burned_with_burned_and_unburned_areas(self): burn_mask = np.array([[1, 0], [0, 1]]) @@ -174,10 +175,14 @@ def test_decay_burned_with_burned_and_unburned_areas(self): ) npt.assert_array_almost_equal( - np.array([[0.753412, ALBEDO_VIS[0][1]], [ALBEDO_VIS[1][0], 0.44346]]), alb_vis, decimal=6 + np.array([[0.753412, ALBEDO_VIS[0][1]], [ALBEDO_VIS[1][0], 0.44346]]), + alb_vis, + decimal=6, ) npt.assert_array_almost_equal( - np.array([[0.659235, ALBEDO_IR[0][1]], [ALBEDO_IR[1][0], 0.354768]]), alb_ir, decimal=6 + np.array([[0.659235, ALBEDO_IR[0][1]], [ALBEDO_IR[1][0], 0.354768]]), + alb_ir, + decimal=6, ) def test_decay_burned_all_burned(self): @@ -210,8 +215,6 @@ def test_decay_burned_k_burned_none(self): burn_mask = np.array([[1, 0], [0, 1]]) with self.assertRaises(ValueError) as context: - decay_burned( - ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, None - ) + decay_burned(ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, None) self.assertIn("k_burned", str(context.exception)) From 43a35e843bffaecfc5e91d8aeb143cd5da825c1c Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Thu, 16 Apr 2026 11:13:51 -0600 Subject: [PATCH 24/30] Distribute - Albedo - Change threshold for the burned and no snowfall mask Changing the threshold for the burn mask to be active to be larger than one full day. At a storm day value of 1.0, the time decay method a decay rate of 0.09817055, versus the exponential burned rate of 0.94176453. This causes a jump and increase of the previous albedo. --- smrf/distribute/albedo.py | 5 +++-- smrf/tests/distribute/test_albedo.py | 10 +++++----- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/smrf/distribute/albedo.py b/smrf/distribute/albedo.py index 1361bd7e..7029df1d 100755 --- a/smrf/distribute/albedo.py +++ b/smrf/distribute/albedo.py @@ -220,12 +220,13 @@ def date_method( :return: alb_v, alb_ir : Decayed albedo for visible and infrared spectrum """ - # Create a mask of burned pixels with snow more than one day ago. + # Create a mask of burned pixels with snow more than one full day ago. + # The comparison uses 1.03 since one hour is ~0.42 in decimal days (1/24) # Only burned pixels that have at least one day since the last snow will # have a different decay rate applied. See :py:meth:`albedo.decay_burned` if self.config.get("post_fire", False): burned_no_snowfall = ne.evaluate( - "(burn_mask == 1) & (storm_day >= 1)", + "(burn_mask == 1) & (storm_day > 1.03)", local_dict={"burn_mask": self.burn_mask, "storm_day": storm_day}, ) else: diff --git a/smrf/tests/distribute/test_albedo.py b/smrf/tests/distribute/test_albedo.py index 6c2b03ec..9cad597c 100644 --- a/smrf/tests/distribute/test_albedo.py +++ b/smrf/tests/distribute/test_albedo.py @@ -37,7 +37,7 @@ DECAY_TIME = pd.to_datetime("2025-04-30") NO_DECAY_TIME = pd.to_datetime("2025-03-01") -STORM_DAYS = np.array([[0.0, 1.0], [1.0, 0.0]]) +STORM_DAYS = np.array([[0.0, 1.042], [1.0, 0.0]]) COS_Z = np.array([[10.0, 10.0], [10.0, 10.0]]) DATA = MagicMock() @@ -254,12 +254,12 @@ def test_date_method_mask_with_post_fire_post_storm( self, mock_decay_power, mock_decay_burned ): """ - Test the window where a snow storm was within the last 24 hrs. There, the burn - decay is noe applied. + Test the window where a snow storm was within the last 24 hrs. + Only after one full day, the burn decay is applied. """ self.subject.config["post_fire"] = True self.subject.config["post_fire_k_burned"] = 1.2 - burn_mask = np.array([[1.0, 1.0], [0.0, 0.0]]) + burn_mask = np.array([[1.0, 1.0], [1.0, 0.0]]) self.subject.burn_mask = burn_mask storm_days = np.array([[0.3, 2.0], [1.0, 0.0]]) @@ -294,7 +294,7 @@ def test_date_method_with_post_fire_post_storm(self): self.subject.config["post_fire"] = True self.subject.config["post_fire_k_burned"] = 1.2 self.subject.burn_mask = np.array([[1.0, 1.0], [1.0, 0.0]]) - storm_days = np.array([[0.3, 2.0], [1.0, 0.0]]) + storm_days = np.array([[0.3, 2.0], [1.042, 0.0]]) current_hours, decay_hours = self.subject.decay_window(DECAY_TIME) expected_mask = np.array([[False, True], [True, False]]) From f7ad81aae7b3c4c5d8a92019286950092fce1a56 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 22 Apr 2026 12:02:26 -0600 Subject: [PATCH 25/30] Tests - Distribute - Add more attributes to the Topo mock That is DEM and vegetation type --- smrf/tests/distribute/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/smrf/tests/distribute/__init__.py b/smrf/tests/distribute/__init__.py index a459de78..d9c5ad30 100644 --- a/smrf/tests/distribute/__init__.py +++ b/smrf/tests/distribute/__init__.py @@ -7,8 +7,10 @@ SKY_VIEW_FACTOR_MOCK = np.array([[1.0, 1.0]]) TOPO_MOCK = MagicMock( spec=Topo, + dem=np.array([[200.0, 150.0], [300.0, 225.0]]), sky_view_factor=SKY_VIEW_FACTOR_MOCK, veg_height=np.array([[5.0, 10.0]]), + veg_type=np.array([[1.0, 1.0], [1.0, 1.0]]), veg_k=np.array([[0.8, 0.1]]), veg_tau=np.array([[0.6, 0.7]]), instance=True, From 0d40ee2f0487e0f2c82865bd0313fd2b2c38eeab Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 22 Apr 2026 12:10:28 -0600 Subject: [PATCH 26/30] Albedo - Post Fire - Change fresh snowfall logic one more time. Third times the charm and this changes the logic of a post fire configuration right after fresh snow. Previously, to catch the lower albedo values late in the melt season after a "reset" event, the post fire decay was only applied after one day since snowfall. This, however, created edge cases where the albedo jumped up right after the one-day mark when the post fire decay rate still was slower than the time decay. To alleviate the need to find a minimum time since last snowfall, this logic now always applies time decay and then the post fire decay for pixels in the burn mask and that are lower than the time decay value. --- smrf/distribute/albedo.py | 27 +++-- smrf/envphys/albedo.py | 69 +++++------- smrf/tests/distribute/test_albedo.py | 163 ++++++++------------------- smrf/tests/envphys/test_albedo.py | 100 ++++++++-------- 4 files changed, 139 insertions(+), 220 deletions(-) diff --git a/smrf/distribute/albedo.py b/smrf/distribute/albedo.py index 7029df1d..7e9fcd70 100755 --- a/smrf/distribute/albedo.py +++ b/smrf/distribute/albedo.py @@ -1,7 +1,6 @@ from datetime import datetime from typing import Tuple -import numexpr as ne import numpy as np import numpy.typing as npt @@ -209,7 +208,7 @@ def date_method( """ Apply a power law decay to the albedo and an optional post fire decay if configured. The post fire decay uses the initial albedo value of the time decay method after - a fresh snowfall. + a fresh snowfall until the power decay reduces the albedo faster. :param alb_v: Visibility albedo :param alb_ir: Infrared albedo @@ -220,17 +219,16 @@ def date_method( :return: alb_v, alb_ir : Decayed albedo for visible and infrared spectrum """ - # Create a mask of burned pixels with snow more than one full day ago. - # The comparison uses 1.03 since one hour is ~0.42 in decimal days (1/24) - # Only burned pixels that have at least one day since the last snow will - # have a different decay rate applied. See :py:meth:`albedo.decay_burned` + if self.config.get("post_fire", False): - burned_no_snowfall = ne.evaluate( - "(burn_mask == 1) & (storm_day > 1.03)", - local_dict={"burn_mask": self.burn_mask, "storm_day": storm_day}, - ) + burned_mask = self.burn_mask == 1 + # Keep originals before power decay + alb_v_initial = alb_v.copy() + alb_ir_initial = alb_ir.copy() else: - burned_no_snowfall = np.zeros_like(storm_day) + burned_mask = None + alb_v_initial = None + alb_ir_initial = None alb_v, alb_ir = albedo.decay_alb_power( self.veg, @@ -240,21 +238,22 @@ def date_method( self.config["date_method_decay_power"], alb_v, alb_ir, - burned_no_snowfall, ) if self.config.get("post_fire", False): - self._logger.debug(" Applying post fire decay") alb_v, alb_ir = albedo.decay_burned( alb_v, alb_ir, + alb_v_initial, + alb_ir_initial, storm_day, - burned_no_snowfall, + burned_mask, self.config.get("post_fire_k_burned", None), ) return alb_v, alb_ir + def decay_window(self, current_timestep: datetime) -> Tuple[float, float]: """ Determine if the current time is within the decay window. diff --git a/smrf/envphys/albedo.py b/smrf/envphys/albedo.py index 6fb864cb..5d8f5e3d 100644 --- a/smrf/envphys/albedo.py +++ b/smrf/envphys/albedo.py @@ -117,7 +117,6 @@ def decay_alb_power( pwr: float, alb_v: npt.NDArray, alb_ir: npt.NDArray, - burned_no_snowfall: npt.NDArray, ) -> Tuple[npt.NDArray, npt.NDArray]: """ Find a decrease in albedo due to litter accumulation. Decay is based on @@ -140,12 +139,9 @@ def decay_alb_power( pwr: power for power law decay alb_v: numpy array of albedo for visible spectrum alb_ir: numpy array of albedo for IR spectrum - burned_no_snowfall: numpy array of burned pixels with at least one day since - last snowfall Returns: Tuple alb_v_d, alb_ir_d : numpy arrays of decayed albedo - """ decay_rates = np.full(veg_type.shape, veg["default"]).astype(np.float32, order="C") @@ -168,27 +164,8 @@ def decay_alb_power( }, ) - # Only apply the standard decay rate if it is not in a pixel that was burned and - # no new snow was deposited. Burned pixels still get the initial albedo after a storm - # in the decay window from the power law - ne.evaluate( - "where(burned_no_snowfall == 0, alb_v - decay_rates, alb_v)", - out=alb_v, - local_dict={ - "alb_v": alb_v, - "decay_rates": decay_rates, - "burned_no_snowfall": burned_no_snowfall, - }, - ) - ne.evaluate( - "where(burned_no_snowfall == 0, alb_ir - decay_rates, alb_ir)", - out=alb_ir, - local_dict={ - "alb_ir": alb_ir, - "decay_rates": decay_rates, - "burned_no_snowfall": burned_no_snowfall, - }, - ) + alb_v -= decay_rates + alb_ir -= decay_rates return alb_v, alb_ir @@ -253,26 +230,27 @@ def decay_alb_hardy(litter, veg_type, storm_day, alb_v, alb_ir): return alb_v_d, alb_ir_d - def decay_burned( alb_v: npt.NDArray, alb_ir: npt.NDArray, + alb_v_initial: npt.NDArray, + alb_ir_initial: npt.NDArray, last_snow: npt.NDArray, - burned_no_snowfall: npt.NDArray, + burn_mask: npt.NDArray, k_burned: float, ) -> Tuple[npt.NDArray, npt.NDArray]: """ Apply exponential albedo decay as a function of time since last snowfall. This changes all pixels which are marked via the burn mask in the topo file and - had at least one day since it snowed. The one day minimum is applied to account - for the `exp()` function, where less than 1 values would cause a growth instead - of a decline. + result in faster decay through this exponential method over the time decay. Args: alb_v: Visible albedo alb_ir: Infrared albedo + alb_v_initial: Visible albedo before the time decay + alb_ir_initial: Infrared albedo after the time decay last_snow: Time since last snow storm (decimal days) - burned_no_snowfall: Mask of burned area with at least one day since snowfall + burn_mask: Mask of burned area from the topo k_burned: Decay rate for burned area Returns: Tuple @@ -284,24 +262,29 @@ def decay_burned( "in the config file." ) - decay_factor = ne.evaluate( - "exp(-where(burn_mask == 1, k_burned, 0) * last_snow)", - local_dict={ - "burn_mask": burned_no_snowfall, - "k_burned": np.float32(k_burned), - "last_snow": last_snow, - }, - ) + # Pre-calculate the exponential decay factor once + decay_factor = np.exp(-k_burned * last_snow) ne.evaluate( - "alb_v * decay_factor", + "where(burn_mask == 1, where(alb_v < (alb_v_initial * decay_factor), alb_v, (alb_v_initial * decay_factor)), alb_v)", out=alb_v, - local_dict={"alb_v": alb_v, "decay_factor": decay_factor}, + local_dict={ + "burn_mask": burn_mask, + "alb_v": alb_v, + "alb_v_initial": alb_v_initial, + "decay_factor": decay_factor, + }, ) ne.evaluate( - "alb_ir * decay_factor", + "where(burn_mask == 1, where(alb_ir < (alb_ir_initial * decay_factor), alb_ir, (alb_ir_initial * decay_factor)), alb_ir)", out=alb_ir, - local_dict={"alb_ir": alb_ir, "decay_factor": decay_factor}, + local_dict={ + "burn_mask": burn_mask, + "alb_ir": alb_ir, + "alb_ir_initial": alb_ir_initial, + "decay_factor": decay_factor, + }, ) + return alb_v, alb_ir diff --git a/smrf/tests/distribute/test_albedo.py b/smrf/tests/distribute/test_albedo.py index 9cad597c..fc3ae3b7 100644 --- a/smrf/tests/distribute/test_albedo.py +++ b/smrf/tests/distribute/test_albedo.py @@ -1,5 +1,5 @@ import unittest -from unittest.mock import ANY, MagicMock, patch +from unittest.mock import MagicMock, patch import numpy as np import numpy.testing as npt @@ -7,12 +7,9 @@ import pytz from smrf.distribute.albedo import Albedo -from smrf.envphys.albedo import decay_alb_power, decay_burned from smrf.tests.smrf_config import SMRFConfig +from smrf.tests.distribute import TOPO_MOCK -TOPO = MagicMock( - veg_type=np.array([[1.0, 1.0], [1.0, 1.0]]), -) CONFIG = { "time": { "time_zone": "utc", @@ -48,11 +45,11 @@ class TestAlbedo(SMRFConfig, unittest.TestCase): def setUp(self): - self.subject = Albedo(CONFIG, TOPO) + self.subject = Albedo(CONFIG, TOPO_MOCK) self.subject.initialize(DATA) def test_init_default(self): - subject = Albedo(CONFIG, TOPO) + subject = Albedo(CONFIG, TOPO_MOCK) self.assertIsNone(subject.albedo) self.assertIsNone(subject.albedo_vis) @@ -93,7 +90,7 @@ def test_distribute_file_broadband(self, mock_read_netcdf): config = self._copy_config(CONFIG) config["albedo"]["source_files"] = "path/to/files" - subject = Albedo(config, TOPO) + subject = Albedo(config, TOPO_MOCK) subject.initialize(pd.DataFrame()) subject.distribute(TIMESTEP, COS_Z, STORM_DAYS) @@ -112,7 +109,7 @@ def test_distribute_file_vis_ir(self, mock_read_netcdf): config = self._copy_config(CONFIG) config["albedo"]["source_files"] = "path/to/files" - subject = Albedo(config, TOPO) + subject = Albedo(config, TOPO_MOCK) subject.initialize(pd.DataFrame()) subject.distribute(TIMESTEP, COS_Z, STORM_DAYS) @@ -136,7 +133,7 @@ def test_distribute_file_direct_diffuse(self, mock_read_netcdf): config = self._copy_config(CONFIG) config["albedo"]["source_files"] = "path/to/files" - subject = Albedo(config, TOPO) + subject = Albedo(config, TOPO_MOCK) subject.initialize(pd.DataFrame()) subject.distribute(TIMESTEP, COS_Z, STORM_DAYS) @@ -191,8 +188,9 @@ def test_distribute_date_method_not_in_window( ) date_decay_method.assert_not_called() + @patch("smrf.envphys.albedo.decay_burned") @patch("smrf.envphys.albedo.decay_alb_power") - def test_date_method_no_post_fire(self, mock_decay_power): + def test_date_method_no_post_fire(self, mock_decay_power, mock_decay_burned): expected_v = ALBEDO_VIS - 0.05 expected_ir = ALBEDO_IR - 0.05 mock_decay_power.return_value = (expected_v, expected_ir) @@ -204,13 +202,34 @@ def test_date_method_no_post_fire(self, mock_decay_power): ALBEDO_VIS, ALBEDO_IR, current_hours, decay_hours, STORM_DAYS ) - # Burned_no_snowfall mask should be all zeros (last argument) - burn_mask = mock_decay_power.call_args[0][-1] - npt.assert_array_equal(burn_mask, np.zeros_like(STORM_DAYS)) + self.assertEqual(mock_decay_power.call_count, 1) + args, _ = mock_decay_power.call_args + + npt.assert_array_equal(args[0], self.subject.veg) + npt.assert_array_equal(args[1], TOPO_MOCK.veg_type) + self.assertEqual(args[2], current_hours) + self.assertEqual(args[3], decay_hours) + npt.assert_array_equal(args[4], self.subject.config["date_method_decay_power"]) + npt.assert_array_equal(args[5], ALBEDO_VIS) + npt.assert_array_equal(args[6], ALBEDO_IR) npt.assert_array_equal(res_v, expected_v) npt.assert_array_equal(res_ir, expected_ir) + mock_decay_burned.asssert_not_called() + + @patch("smrf.envphys.albedo.decay_alb_power") + def test_date_method_no_post_fire_config(self, mock_decay_power): + if "post_fire" in self.subject.config: + del self.subject.config["post_fire"] + + mock_decay_power.return_value = MagicMock(), MagicMock() + + self.subject.date_method( + ALBEDO_VIS, ALBEDO_IR, 1, 2, STORM_DAYS + ) + mock_decay_power.assert_called_once() + @patch("smrf.envphys.albedo.decay_burned") @patch("smrf.envphys.albedo.decay_alb_power") def test_date_method_with_post_fire(self, mock_decay_power, mock_decay_burned): @@ -223,9 +242,9 @@ def test_date_method_with_post_fire(self, mock_decay_power, mock_decay_burned): power_ir = ALBEDO_IR - 0.02 mock_decay_power.return_value = (power_v, power_ir) - final_v = power_v - 0.1 - final_ir = power_ir - 0.1 - mock_decay_burned.return_value = (final_v, final_ir) + final_vis = MagicMock("AlBEDO_VIS") + final_ir = MagicMock("ALBEDO_IR") + mock_decay_burned.return_value = final_vis, final_ir current_hours, decay_hours = self.subject.decay_window(DECAY_TIME) @@ -233,104 +252,16 @@ def test_date_method_with_post_fire(self, mock_decay_power, mock_decay_burned): ALBEDO_VIS, ALBEDO_IR, current_hours, decay_hours, STORM_DAYS ) - expected_mask = np.array([[False, True], [False, False]]) - received_mask = mock_decay_power.call_args[0][-1] - npt.assert_array_equal(received_mask, expected_mask) - - mock_decay_burned.assert_called_once_with( - power_v, - power_ir, - STORM_DAYS, - ANY, - self.subject.config["post_fire_k_burned"], - ) - npt.assert_array_equal(mock_decay_burned.call_args[0][3], expected_mask) - - npt.assert_array_equal(res_v, final_v) - - @patch("smrf.envphys.albedo.decay_burned") - @patch("smrf.envphys.albedo.decay_alb_power") - def test_date_method_mask_with_post_fire_post_storm( - self, mock_decay_power, mock_decay_burned - ): - """ - Test the window where a snow storm was within the last 24 hrs. - Only after one full day, the burn decay is applied. - """ - self.subject.config["post_fire"] = True - self.subject.config["post_fire_k_burned"] = 1.2 - burn_mask = np.array([[1.0, 1.0], [1.0, 0.0]]) - self.subject.burn_mask = burn_mask - storm_days = np.array([[0.3, 2.0], [1.0, 0.0]]) - - # Preserve the method return signature of mocked methods - power_v = 1 - power_ir = 1 - mock_decay_power.return_value = (power_v, power_ir) - mock_decay_burned.return_value = (1, 1) - - current_hours, decay_hours = self.subject.decay_window(DECAY_TIME) - - self.subject.date_method( - ALBEDO_VIS, ALBEDO_IR, current_hours, decay_hours, storm_days - ) - - expected_mask = np.array([[False, True], [False, False]]) - npt.assert_array_equal(mock_decay_power.call_args[0][-1], expected_mask) - - mock_decay_burned.assert_called_once_with( - power_v, - power_ir, - storm_days, - ANY, - self.subject.config["post_fire_k_burned"], - ) - npt.assert_array_equal(mock_decay_burned.call_args[0][3], expected_mask) - - def test_date_method_with_post_fire_post_storm(self): - """ - Integration test that uses both decay methods - """ - self.subject.config["post_fire"] = True - self.subject.config["post_fire_k_burned"] = 1.2 - self.subject.burn_mask = np.array([[1.0, 1.0], [1.0, 0.0]]) - storm_days = np.array([[0.3, 2.0], [1.042, 0.0]]) - - current_hours, decay_hours = self.subject.decay_window(DECAY_TIME) - expected_mask = np.array([[False, True], [True, False]]) - power_vis, power_ir = decay_alb_power( - self.subject.veg, - self.subject.veg_type, - current_hours, - decay_hours, - self.subject.config["date_method_decay_power"], - ALBEDO_VIS.copy(), - ALBEDO_IR.copy(), - expected_mask - ) - burn_vis, burn_ir = decay_burned( - power_vis.copy(), - power_ir.copy(), - storm_days, - expected_mask, - self.subject.config["post_fire_k_burned"], - ) + self.assertEqual(mock_decay_burned.call_count, 1) + args, _ = mock_decay_burned.call_args - albedo_v, albedo_ir = self.subject.date_method( - ALBEDO_VIS.copy(), ALBEDO_IR.copy(), current_hours, decay_hours, storm_days - ) + npt.assert_array_equal(args[0], power_v) + npt.assert_array_equal(args[1], power_ir) + npt.assert_array_equal(args[2], ALBEDO_VIS) + npt.assert_array_equal(args[3], ALBEDO_IR) + npt.assert_array_equal(args[4], STORM_DAYS) + npt.assert_array_equal(args[5], self.subject.burn_mask) + self.assertEqual(args[6], 1.2) - npt.assert_array_equal( - albedo_v, - np.array([ - [power_vis[0][0], burn_vis[0][1]], - [burn_vis[1][0], power_vis[1][1]] - ]), - ) - npt.assert_array_equal( - albedo_ir, - np.array([ - [power_ir[0][0], burn_ir[0][1]], - [burn_ir[1][0], power_ir[1][1]] - ]), - ) + self.assertIs(res_v, final_vis) + self.assertIs(res_ir, final_ir) diff --git a/smrf/tests/envphys/test_albedo.py b/smrf/tests/envphys/test_albedo.py index 3c764856..f00dfd83 100644 --- a/smrf/tests/envphys/test_albedo.py +++ b/smrf/tests/envphys/test_albedo.py @@ -53,7 +53,6 @@ def test_max_decay_after_end(self): pwr, ALBEDO_VIS.copy(), ALBEDO_IR.copy(), - BURNED_NO_SNOWFALL, ) npt.assert_array_almost_equal(expected_v, alb_v_d) @@ -86,44 +85,6 @@ def test_power_decay_during_window(self): pwr, ALBEDO_VIS.copy(), ALBEDO_IR.copy(), - BURNED_NO_SNOWFALL, - ) - - npt.assert_array_almost_equal(expected_v, alb_v_d, decimal=5) - npt.assert_array_almost_equal(expected_ir, alb_ir_d, decimal=5) - - def test_power_decay_during_window_with_burn(self): - """Burned pixels will not be changed by power decay after snowfall""" - current_hours = 50 - decay_hours = 100 - pwr = 2.0 - burned_no_snowfall = np.array([[0.0, 1.0], [0.0, 0.0]]) - - # Decay rates - tao_default = decay_hours / (self.veg["default"] ** (1.0 / pwr)) - tao_41 = decay_hours / (self.veg["41"] ** (1.0 / pwr)) - tao_42 = decay_hours / (self.veg["42"] ** (1.0 / pwr)) - - expected_decay = np.zeros_like(self.veg_type, dtype=float) - expected_decay[self.veg_type == 0] = (current_hours / tao_default) ** pwr - expected_decay[self.veg_type == 41] = (current_hours / tao_41) ** pwr - expected_decay[self.veg_type == 42] = (current_hours / tao_42) ** pwr - - # Burned pixels should not be changed by power decay - expected_v = ALBEDO_VIS - expected_decay - expected_v[0][1] = ALBEDO_VIS[0][1] - expected_ir = ALBEDO_IR - expected_decay - expected_ir[0][1] = ALBEDO_IR[0][1] - - alb_v_d, alb_ir_d = decay_alb_power( - self.veg, - self.veg_type, - current_hours, - decay_hours, - pwr, - ALBEDO_VIS.copy(), - ALBEDO_IR.copy(), - burned_no_snowfall, ) npt.assert_array_almost_equal(expected_v, alb_v_d, decimal=5) @@ -158,7 +119,6 @@ def test_partial_coverage_during_decay(self): pwr, ALBEDO_VIS.copy(), ALBEDO_IR.copy(), - BURNED_NO_SNOWFALL, ) npt.assert_array_almost_equal(expected_v, alb_v_d, decimal=5) @@ -166,12 +126,16 @@ def test_partial_coverage_during_decay(self): class TestDecayBurned(unittest.TestCase): + """ + NOTE: All arrays passed in as arguments are updated in place, so pass a `copy()` + """ def test_decay_burned_with_burned_and_unburned_areas(self): burn_mask = np.array([[1, 0], [0, 1]]) - # Arrays are updated in place, so pass a copy alb_vis, alb_ir = decay_burned( - ALBEDO_VIS.copy(), ALBEDO_IR.copy(), LAST_SNOW, burn_mask, K_BURNED + ALBEDO_VIS.copy(), ALBEDO_IR.copy(), + ALBEDO_VIS.copy(), ALBEDO_IR.copy(), + LAST_SNOW, burn_mask, K_BURNED ) npt.assert_array_almost_equal( @@ -188,9 +152,10 @@ def test_decay_burned_with_burned_and_unburned_areas(self): def test_decay_burned_all_burned(self): burn_mask = np.ones_like(ALBEDO_VIS) - # Arrays are updated in place, so pass a copy alb_vis, alb_ir = decay_burned( - ALBEDO_VIS.copy(), ALBEDO_IR.copy(), LAST_SNOW, burn_mask, K_BURNED + ALBEDO_VIS.copy(), ALBEDO_IR.copy(), + ALBEDO_VIS.copy(), ALBEDO_IR.copy(), + LAST_SNOW, burn_mask, K_BURNED ) npt.assert_array_almost_equal( @@ -203,18 +168,59 @@ def test_decay_burned_all_burned(self): def test_decay_burned_all_unburned(self): burn_mask = np.zeros_like(ALBEDO_VIS) - # Arrays are updated in place, so pass a copy alb_vis, alb_ir = decay_burned( - ALBEDO_VIS.copy(), ALBEDO_IR.copy(), LAST_SNOW, burn_mask, K_BURNED + ALBEDO_VIS.copy(), ALBEDO_IR.copy(), + ALBEDO_VIS.copy(), ALBEDO_IR.copy(), + LAST_SNOW, burn_mask, K_BURNED ) npt.assert_array_equal(ALBEDO_VIS, alb_vis) npt.assert_array_equal(ALBEDO_IR, alb_ir) + def test_decay_burned_initial_lower_than_calculated(self): + """ + Test that the minimum of current albedo and calculated burn decay is taken. + Tests the nested where condition: where(alb_v < calculated, alb_v, calculated) + """ + burn_mask = np.ones_like(ALBEDO_VIS) + k_burned = 0.1 + last_snow = np.array([[10.0, 10.0], [10.0, 10.0]]) + + alb_v_current = np.array([[0.2, 0.6], [0.7, 0.1]]) + # [0,0]: Time decay + # [0,1]: Burn decay + # [1,0]: Burn decay + # [1,1]: Time decay + + initial_v = ALBEDO_VIS.copy() + initial_ir = ALBEDO_IR.copy() + + decay_factor = np.exp(-k_burned * last_snow) + expected_v = np.minimum(alb_v_current, initial_v * decay_factor) + # ALBEDO_IR values are all higher than decay here + expected_ir = np.minimum(ALBEDO_IR, initial_ir * decay_factor) + + alb_v_out, alb_ir_out = decay_burned( + alb_v_current.copy(), ALBEDO_IR.copy(), + initial_v, initial_ir, + last_snow, burn_mask, k_burned + ) + + npt.assert_array_almost_equal(expected_v, alb_v_out, decimal=6) + npt.assert_array_almost_equal(expected_ir, alb_ir_out, decimal=6) + def test_decay_burned_k_burned_none(self): burn_mask = np.array([[1, 0], [0, 1]]) with self.assertRaises(ValueError) as context: - decay_burned(ALBEDO_VIS, ALBEDO_IR, LAST_SNOW, burn_mask, None) + decay_burned( + ALBEDO_VIS.copy(), + ALBEDO_IR.copy(), + ALBEDO_VIS.copy(), + ALBEDO_IR.copy(), + LAST_SNOW, + burn_mask, + None, + ) self.assertIn("k_burned", str(context.exception)) From 4bd315c5818b29c00e8c05eec974d5986a3a55f7 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 20 May 2026 08:33:12 -0600 Subject: [PATCH 27/30] Envphys - Constants - Re-use already defined values for albedo Great catch by @brentwilder that the albedo module had the same constants defined as the global constant file. Used this chance to rename them to (hopefully) be a bit more readable. Also organized the constants file and group values with sorting by alphabetically order. --- smrf/envphys/albedo.py | 18 ++++-------- smrf/envphys/constants.py | 52 ++++++++++++++++++----------------- smrf/envphys/solar/toporad.py | 6 ++-- 3 files changed, 35 insertions(+), 41 deletions(-) diff --git a/smrf/envphys/albedo.py b/smrf/envphys/albedo.py index 5d8f5e3d..0d2f189b 100644 --- a/smrf/envphys/albedo.py +++ b/smrf/envphys/albedo.py @@ -5,15 +5,7 @@ import numpy as np import numpy.typing as npt -MAXV = 1.0 # vis albedo when gsize = 0 -MAXIR = 0.85447 # IR albedo when gsize = 0 -IRFAC = -0.02123 # IR decay factor -VFAC = 500.0 # visible decay factor -VZRG = 1.375e-3 # vis zenith increase range factor -IRZRG = 2.0e-3 # ir zenith increase range factor -IRZ0 = 0.1 # ir zenith increase range, gsize=0 -BOIL = 373.15 # boiling temperature K -GRAVITY = 9.80665 # gravity (m/s^2) +from .constants import IR_FACTOR, IR_Z_0, IR_Z_RF, IR_MAX_0, VIS_MAX_0, VIS_FACTOR, VIS_Z_RF # TODO - Need to raise here instead of silently fail to make the user aware @@ -88,14 +80,14 @@ def albedo( gir = radius_ir + (range_ir * growth_factor) # calc albedo for cos(z)=1 - alb_v_1 = MAXV - (gv / VFAC) - alb_ir_1 = MAXIR * np.exp(IRFAC * gir) + alb_v_1 = VIS_MAX_0 - (gv / VIS_FACTOR) + alb_ir_1 = IR_MAX_0 * np.exp(IR_FACTOR * gir) # calculate effect of cos(z)<1 # adjust diurnal increase range - dzv = gv * VZRG - dzir = (gir * IRZRG) + IRZ0 + dzv = gv * VIS_Z_RF + dzir = (gir * IR_Z_RF) + IR_Z_0 # calculate albedo alb_v = alb_v_1 diff --git a/smrf/envphys/constants.py b/smrf/envphys/constants.py index b21b9cd2..98ea496b 100644 --- a/smrf/envphys/constants.py +++ b/smrf/envphys/constants.py @@ -1,29 +1,31 @@ -VISIBLE_MIN = .28 -VISIBLE_MAX = .7 -IR_MIN = .7 +# Albedo related +IR_FACTOR = -0.02123 # IR decay factor IR_MAX = 2.8 -# visible solar irradiance wavelengths -VISIBLE_WAVELENGTHS = [VISIBLE_MIN, VISIBLE_MAX] -# infrared solar irradiance wavelengths +IR_MAX_0 = 0.85447 # IR albedo when grain_size = 0 +IR_MIN = .7 +IR_Z_RF = 2.0e-3 # IR zenith increase range factor +IR_Z_0 = 0.1 # IR zenith increase range, grain_size = 0 +VIS_FACTOR = 500.0 # Visible decay factor +VIS_MAX = .7 +VIS_MAX_0 = 1.0 # Visible albedo when grain_size = 0 +VIS_MIN = .28 +VIS_Z_RF = 1.375e-3 # Visible zenith increase range factor + +# Visible solar irradiance wavelengths +VISIBLE_WAVELENGTHS = [VIS_MIN, VIS_MAX] +# Infrared solar irradiance wavelengths IR_WAVELENGTHS = [IR_MIN, IR_MAX] -SOLAR_CONSTANT = 1368.0 # solar constant in W/m**2 -MAXV = 1.0 # vis albedo when gsize = 0 -MAXIR = 0.85447 # IR albedo when gsize = 0 -IRFAC = -0.02123 # IR decay factor -VFAC = 500.0 # visible decay factor -VZRG = 1.375e-3 # vis zenith increase range factor -IRZRG = 2.0e-3 # ir zenith increase range factor -IRZ0 = 0.1 # ir zenith increase range, gsize=0 -STEF_BOLTZ = 5.6697e-8 # Stefan Boltzmann constant -EMISS_TERRAIN = 0.98 # emissivity of the terrain -EMISS_VEG = 0.96 # emissivity of the vegetation -FREEZE = 273.16 # freezing temp K -BOIL = 373.15 # boiling temperature K -STD_LAPSE_M = -0.0065 # lapse rate (K/m) -STD_LAPSE = -6.5 # lapse rate (K/km) -SEA_LEVEL = 1.013246e5 # sea level pressure -RGAS = 8.31432e3 # gas constant (J / kmole / deg) -GRAVITY = 9.80665 # gravity (m/s^2) -MOL_AIR = 28.9644 # molecular weight of air (kg / kmole) +BOIL = 373.15 # Boiling temperature K +EMISS_TERRAIN = 0.98 # Emissivity of the terrain +EMISS_VEG = 0.96 # Emissivity of the vegetation +FREEZE = 273.16 # Freezing temperature K +GRAVITY = 9.80665 # Gravity (m/s^2) +MOL_AIR = 28.9644 # Molecular weight of air (kg / kmole) +RGAS = 8.31432e3 # Gas constant (J / kmole / deg) +SEA_LEVEL = 1.013246e5 # Sea level pressure +SOLAR_CONSTANT = 1368.0 # Solar constant in W/m**2 STD_AIRTMP = 2.88e2 +STD_LAPSE = -6.5 # Lapse rate (K/km) +STD_LAPSE_M = -0.0065 # Lapse rate (K/m) +STEF_BOLTZ = 5.6697e-8 # Stefan Boltzmann constant diff --git a/smrf/envphys/solar/toporad.py b/smrf/envphys/solar/toporad.py index 92eceea6..5476a0e2 100644 --- a/smrf/envphys/solar/toporad.py +++ b/smrf/envphys/solar/toporad.py @@ -11,8 +11,8 @@ SEA_LEVEL, STD_AIRTMP, STD_LAPSE, - VISIBLE_MAX, - VISIBLE_MIN, + VIS_MAX, + VIS_MIN, ) from smrf.envphys.solar.irradiance import direct_solar_irradiance from smrf.envphys.solar.twostream import twostream @@ -21,7 +21,7 @@ def check_wavelengths(wavelength_range): - if wavelength_range[0] >= VISIBLE_MIN and wavelength_range[1] <= VISIBLE_MAX: + if wavelength_range[0] >= VIS_MIN and wavelength_range[1] <= VIS_MAX: return elif wavelength_range[0] >= IR_MIN and wavelength_range[1] <= IR_MAX: return From 465daa6dd71d23139f85a924987d21240ca1ec2d Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 20 May 2026 08:47:07 -0600 Subject: [PATCH 28/30] Albedo - Move post_fire_k_burned check to distribute modules Great feedback by @brentwilder that math should be done in the envphys module while distribute checks for presence of all relevant parameters. --- smrf/distribute/albedo.py | 10 +++++++++- smrf/envphys/albedo.py | 6 ------ smrf/tests/distribute/test_albedo.py | 11 +++++++++++ smrf/tests/envphys/test_albedo.py | 16 ---------------- 4 files changed, 20 insertions(+), 23 deletions(-) diff --git a/smrf/distribute/albedo.py b/smrf/distribute/albedo.py index 7e9fcd70..df8e4e5e 100755 --- a/smrf/distribute/albedo.py +++ b/smrf/distribute/albedo.py @@ -225,8 +225,16 @@ def date_method( # Keep originals before power decay alb_v_initial = alb_v.copy() alb_ir_initial = alb_ir.copy() + + k_burned = self.config.get("post_fire_k_burned", None) + if k_burned is None: + raise ValueError( + "Post fire albedo decay is configured, but post_fire_k_burned is " + "not set in the config file." + ) else: burned_mask = None + k_burned = None alb_v_initial = None alb_ir_initial = None @@ -248,7 +256,7 @@ def date_method( alb_ir_initial, storm_day, burned_mask, - self.config.get("post_fire_k_burned", None), + k_burned, ) return alb_v, alb_ir diff --git a/smrf/envphys/albedo.py b/smrf/envphys/albedo.py index 0d2f189b..463dcaca 100644 --- a/smrf/envphys/albedo.py +++ b/smrf/envphys/albedo.py @@ -248,12 +248,6 @@ def decay_burned( Returns: Tuple alb_v_d, alb_ir_d : numpy arrays of decayed albedo """ - if k_burned is None: - raise ValueError( - "Post fire albedo decay is configured, but post_fire_k_burned is not configured " - "in the config file." - ) - # Pre-calculate the exponential decay factor once decay_factor = np.exp(-k_burned * last_snow) diff --git a/smrf/tests/distribute/test_albedo.py b/smrf/tests/distribute/test_albedo.py index fc3ae3b7..aa29400f 100644 --- a/smrf/tests/distribute/test_albedo.py +++ b/smrf/tests/distribute/test_albedo.py @@ -265,3 +265,14 @@ def test_date_method_with_post_fire(self, mock_decay_power, mock_decay_burned): self.assertIs(res_v, final_vis) self.assertIs(res_ir, final_ir) + + def test_date_method_with_post_fire_missing_k(self): + self.subject.config["post_fire"] = True + self.subject.config["post_fire_k_burned"] = None + + with self.assertRaises(ValueError) as context: + self.subject.date_method( + ALBEDO_VIS, ALBEDO_IR, 1, 1, STORM_DAYS + ) + + self.assertIn("post_fire_k_burned is not set", str(context.exception)) diff --git a/smrf/tests/envphys/test_albedo.py b/smrf/tests/envphys/test_albedo.py index f00dfd83..ab31eb34 100644 --- a/smrf/tests/envphys/test_albedo.py +++ b/smrf/tests/envphys/test_albedo.py @@ -208,19 +208,3 @@ def test_decay_burned_initial_lower_than_calculated(self): npt.assert_array_almost_equal(expected_v, alb_v_out, decimal=6) npt.assert_array_almost_equal(expected_ir, alb_ir_out, decimal=6) - - def test_decay_burned_k_burned_none(self): - burn_mask = np.array([[1, 0], [0, 1]]) - - with self.assertRaises(ValueError) as context: - decay_burned( - ALBEDO_VIS.copy(), - ALBEDO_IR.copy(), - ALBEDO_VIS.copy(), - ALBEDO_IR.copy(), - LAST_SNOW, - burn_mask, - None, - ) - - self.assertIn("k_burned", str(context.exception)) From 1b9d28fe7b0def67e3b84a662d45d2add714a6b0 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 20 May 2026 08:59:57 -0600 Subject: [PATCH 29/30] CoreConfig - Add TODO to post_fire_k_burned parameter --- smrf/framework/CoreConfig.ini | 1 + 1 file changed, 1 insertion(+) diff --git a/smrf/framework/CoreConfig.ini b/smrf/framework/CoreConfig.ini index 75d120cc..ae7b194d 100755 --- a/smrf/framework/CoreConfig.ini +++ b/smrf/framework/CoreConfig.ini @@ -871,6 +871,7 @@ description = Apply a different post fire decay method for areas masked in the t file via a 'burn_mask' layer. This also requires to set the 'post_fire_k_burned' value in this configuration section +# Issue #74: Add paper citation describing this parameter post_fire_k_burned: default = None, type = float, From 224ace08674f6c079b5a9fd92b30261fdfa5d7f2 Mon Sep 17 00:00:00 2001 From: Joachim Meyer Date: Wed, 20 May 2026 09:10:02 -0600 Subject: [PATCH 30/30] Distribute - Albedo - Improve docstring describing the date_method --- smrf/distribute/albedo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/smrf/distribute/albedo.py b/smrf/distribute/albedo.py index df8e4e5e..2f81ef7b 100755 --- a/smrf/distribute/albedo.py +++ b/smrf/distribute/albedo.py @@ -208,7 +208,7 @@ def date_method( """ Apply a power law decay to the albedo and an optional post fire decay if configured. The post fire decay uses the initial albedo value of the time decay method after - a fresh snowfall until the power decay reduces the albedo faster. + a fresh snowfall until the burn decay reduces the albedo faster. :param alb_v: Visibility albedo :param alb_ir: Infrared albedo