Merge pull request #1 from UofU-Cryosphere/storm_day_comparison

Storm day fix

Author: jomey

smrf/distribute/precipitation.py

@@ -1,7 +1,8 @@
+from datetime import timedelta
 
 import netCDF4 as nc
 import numpy as np
-
+from dateutil.parser import parse
 from smrf.distribute import image_data
 from smrf.envphys import precip, Snow, storms
 from smrf.utils import utils
@@ -126,35 +127,42 @@ def initialize(self, topo, data, date_time=None):
         self.snow_density = np.zeros((topo.ny, topo.nx))
         self.storm_days = np.zeros((topo.ny, topo.nx))
         self.storm_total = np.zeros((topo.ny, topo.nx))
-        self.last_storm_day = np.zeros((topo.ny, topo.nx))
         self.dem = topo.dem
 
         # Assign storm_days array if given
         if self.config["storm_days_restart"] is not None:
             self._logger.debug('Reading {} from {}'.format(
-                'storm_days', self.config['storm_days_restart']))
-            f = nc.Dataset(self.config['storm_days_restart'], 'r')
+                'storm_days', self.config["storm_days_restart"])
+            )
+            f = nc.Dataset(self.config["storm_days_restart"])
             f.set_always_mask(False)
 
             if 'storm_days' in f.variables:
                 t = f.variables['time']
+                t_max = t[:].max()
                 time = nc.num2date(
-                    t[:],
+                    t_max,
                     t.getncattr('units'),
                     t.getncattr('calendar'),
-                    only_use_cftime_datetimes=False)
-                time = np.array(
-                    [ti.replace(tzinfo=self.start_date.tzinfo) for ti in time])
-                time_ind = np.where(time == self.start_date.to_pydatetime())[0]
+                    only_use_cftime_datetimes=False,
+                    only_use_python_datetimes=True
+                )
+                # Check whether the last storm day entry and the start
+                # of this run is an hour apart (3600 seconds)
+                max_time = time.replace(tzinfo=self.start_date.tzinfo)
+                delta_seconds = (
+                    self.start_date.to_pydatetime() - parse(str(max_time))
+                )
 
-                if time_ind.size == 0:
+                # Python timedelta are handled in seconds and days
+                if delta_seconds > timedelta(seconds=(self.time_step*60)):
                     self._logger.warning(
                         'Invalid storm_days input! Setting to 0.0')
                     self.storm_days = np.zeros((topo.ny, topo.nx))
 
                 else:
                     # start at index of storm_days - 1
-                    self.storm_days = f.variables['storm_days'][time_ind - 1, :, :][0]  # noqa
+                    self.storm_days = f.variables['storm_days'][t_max]
             else:
                 self._logger.error(
                     'Variable storm_days not in {}'.format(
@@ -401,10 +409,6 @@ def distribute_for_marks2017(self, data, precip_temp, ta, time):
         self.percent_snow = perc_snow
         self.snow_density = snow_den
 
-        # day of last storm, this will be used in albedo
-        self.last_storm_day = utils.water_day(data.name)[0] - \
-            self.storm_days - 0.001
-
     def distribute_for_susong1999(self, data, ppt_temp, time):
         """Susong 1999 estimates percent snow and snow density based on
         Susong et al, (1999) :cite:`Susong&al:1999`.
@@ -430,11 +434,11 @@ def distribute_for_susong1999(self, data, ppt_temp, time):
             stormDays, stormPrecip = storms.time_since_storm(
                 self.precip,
                 perc_snow,
+                storm_days=self.storm_days,
+                storm_precip=self.storm_total,
                 time_step=self.time_step/60/24,
-                mass=self.ppt_threshold,
-                time=self.time_to_end_storm,
-                stormDays=self.storm_days,
-                stormPrecip=self.storm_total)
+                mass_threshold=self.ppt_threshold,
+            )
 
             # save the model state
             self.percent_snow = perc_snow
@@ -451,10 +455,6 @@ def distribute_for_susong1999(self, data, ppt_temp, time):
             self.percent_snow = np.zeros(self.storm_days.shape)
             self.snow_density = np.zeros(self.storm_days.shape)
 
-        # day of last storm, this will be used in albedo
-        self.last_storm_day = utils.water_day(data.name)[0] - \
-            self.storm_days - 0.001
-
     def distribute_thread(self, smrf_queue, data_queue):
         """
         Distribute the data using threading and queue. All data is provided and

smrf/envphys/storms.py

@@ -2,69 +2,68 @@
 import pandas as pd
 
 
-def time_since_storm(precipitation, perc_snow, time_step=1/24, mass=1.0,
-                     time=4, stormDays=None, stormPrecip=None, ps_thresh=0.5):
+def time_since_storm(
+    precipitation,
+    percent_snow_precipitation,
+    storm_days,
+    storm_precip,
+    time_step,
+    mass_threshold=1.0,
+    percent_snow_threshold=0.5
+):
     """
-    Calculate the decimal days since the last storm given a precip time series,
-    percent snow, mass threshold, and time threshold
+    Increase or reset the given storm days and storm precipitation based of
+    given percent snow and mass threshold.
 
-     - Will look for pixels where perc_snow > 50% as storm locations
-     - A new storm will start if the mass at the pixel has exceeded the mass
-         limit, this ensures that the enough has accumulated
+    Steps:
+     - Look for pixels where percent snow precipitation is above the threshold.
+     - Reset storm days if the precipitation at the pixel is above the mass
+         threshold or increase the counter if it is below.
+     - Add the precipitation amount to the storm precipitation.
+
+    *NOTE*: Storm day precipitation is initialized and set to 0 with each
+        processed day. Hence, no cross-day storm tracking is possible if the
+        period between the days falls below the threshold.
 
     Args:
         precipitation: Precipitation values
-        perc_snow: Percent of precipitation that was snow
+        percent_snow_precipitation: Percent of precipitation that was snow
+        storm_days: Storm days to keep track of
+        storm_precip: Keeps track of the total storm precip
         time_step: Step in days of the model run
-        mass: Threshold for the mass to start a new storm
-        time: Threshold for the time to start a new storm
-        stormDays: If specified, this is the output from a previous run of
-            storms else it will be set to the date_time value
-        stormPrecip: Keeps track of the total storm precip
+        (Optional)
+        mass_threshold: Minimum amount of precipitation required to be a storm
+                        (snow mass). Default: 0.5
+        percent_snow_threshold: Minimum fraction for values in 
+                                `percent_snow_precipitation` to be considered 
+                                a snow event. Default: 0.5 (50%)
 
     Returns:
         tuple:
-        - **stormDays** - Array representing the days since the last storm at
-            a pixel
-        - **stormPrecip** - Array representing the precip accumulated during
-            the most recent storm
+        - **stormDays** - Updated storm days
+        - **stormPrecip** - Updated storm precipitation
 
-    Created Janurary 5, 2016
+    Created January 5, 2016
     @author: Scott Havens
-    """
-    # either preallocate or use the input
-    if stormDays is None:
-        stormDays = np.zeros(precipitation.shape)
-
-    if stormPrecip is None:
-        stormPrecip = np.zeros(precipitation.shape)
-
-    # if there is no snow, don't reset the counter
-    # This ensures that the albedo won't be reset
-    stormDays += 1
-    if np.sum(perc_snow) == 0:
-        #         stormDays = np.add(stormDays, 1)
-        stormPrecip = np.zeros(precipitation.shape)
-        return stormDays, stormPrecip
 
-    # determine locations where it has snowed
-    idx = perc_snow >= ps_thresh
-
-    # determine locations where the time threshold has passed
-    # these areas, the stormPrecip will be set back to zero
-    idx_time = stormDays >= time
-    stormPrecip[idx_time] = 0
+    Updated: February 07, 2022
+    @author: Joachim Meyer, Dillon Ragar
+    """
 
-    # add the values to the stormPrecip
-    stormPrecip[idx] = + precipitation[idx]
+    # Step 1: Pixels above snow percent threshold
+    location_index = (percent_snow_precipitation >= percent_snow_threshold)
+    storm_precip[location_index] += precipitation[location_index]
 
-    # see if the mass threshold has been passed
-    idx_mass = stormPrecip >= mass
+    # Step 2: Reset locations above mass threshold or increase counter when
+    # below
+    location_index = (storm_precip >= mass_threshold)
+    storm_days[location_index] = 0
+    storm_days[~location_index] += time_step
 
-    # reset the stormDays to zero where the storm is present
-    stormDays[idx_mass] = 0
+    # Step 3: Increase the storm precipitation total for the day
+    storm_precip[~location_index] = 0
 
-    return stormDays, stormPrecip
+    return storm_days, storm_precip
 
 
 def time_since_storm_pixel(precipitation, dpt, perc_snow, storming,