Merge pull request #236 from UW-Hydro/develop

arbennett · web-flow · commit c7fd5d53b658 · 2020-08-01T13:54:30.000-07:00
Develop
diff --git a/README.md b/README.md
@@ -145,6 +145,11 @@ Or in BibTeX:
 }
 ```
 
+Acknowledgements
+================
+MetSim has greatly benefited from the user community, who have contributed code, tested features, provided feedback, and helped with documentation.
+We would like to thank and acknowledge the work of Ted Bohn, Andy Wood, Kristen Whitney, Yifan Cheng, Liz Clark, Oriana Chegwidden, Ethan Gutmann, Kostas Andreadis, Thomas Remke, Ed Maurer, and Philipp Sommer for their help.
+
 
 References
 ==========
diff --git a/docs/whats-new.rst b/docs/whats-new.rst
@@ -5,6 +5,17 @@ What's New
 
 .. _whats-new.2.3.0:
 
+v2.3.2
+------
+
+Bug fixes
+~~~~~~~~~
+- Fix a bug in ascii input reading due to pandas change
+
+Enhancements
+~~~~~~~~~~~~
+- Drastically speed up PITRI precipitation disaggregation
+
 v2.3.1
 ------
 
@@ -15,7 +26,7 @@ Bug fixes
 
 v2.3.0
 ------
-Enchancements
+Enhancements
 ~~~~~~~~~~~~~
 - Allow for variable renaming in INI configuration files.
 - Added capability for new YAML configuration file format
diff --git a/metsim/disaggregate.py b/metsim/disaggregate.py
@@ -107,9 +107,8 @@ def disaggregate(df_daily: pd.DataFrame, params: dict,
     df_disagg['longwave'] = longwave(
         df_disagg['temp'].values, df_disagg['vapor_pressure'].values,
         df_disagg['tskc'].values, params)
-
-    df_disagg['prec'] = prec(df_daily['prec'].values, df_daily['t_min'].values,
-                             ts, params, df_daily.index.month)
+    df_disagg['prec'] = prec(df_daily['prec'], df_daily['t_min'], ts, params,
+                             df_daily.get('t_pk'), df_daily.get('dur'))
 
     if 'wind' in df_daily:
         df_disagg['wind'] = wind(df_daily['wind'].values, ts, params)
@@ -252,7 +251,7 @@ def temp(t_min: np.array, t_max: np.array, out_len: int,
 
 
 def prec(prec: pd.Series, t_min: pd.Series, ts: float, params: dict,
-         month_of_year: int):
+         t_pk=None, dur=None):
     """
     Distributes sub-daily precipitation either evenly (uniform) or with a
     triangular (triangle) distribution, depending upon the chosen method.
@@ -281,127 +280,59 @@ def prec(prec: pd.Series, t_min: pd.Series, ts: float, params: dict,
     prec:
         A sub-daily timeseries of precipitation. [mm]
     """
-    def prec_TRIANGLE(prec: pd.Series, t_min: pd.Series,
-                      month_of_year: int, do_mix: bool, params: dict):
-
-        def P_kernel(t_corners, m, t):
-            # calculating precipitation intensity of current time t
-            t_start = t_corners[0]
-            t_pk = t_corners[1]
-            t_end = t_corners[2]
-            if (t < t_start) | (t > t_end):
-                P = 0
-            elif t <= t_pk:
-                P = m * (t - t_start)
-            elif t >= t_pk:
-                P = m * (t_end - t)
-            return P
-
-        dur = params['dur']
-        t_pk = params['t_pk']
-        n_days = len(prec)
-        steps_per_day = cnst.MIN_PER_HOUR * cnst.HOURS_PER_DAY
-        offset = np.ceil(steps_per_day / 2)
-        output_length = int(steps_per_day * n_days)
-        index = np.arange(output_length)
-        steps_per_two_days = int(((np.ceil(steps_per_day / 2)) * 2) +
-                                 steps_per_day)
-        # time of step on next day [minutes]
-        t_next = 2 * cnst.MIN_PER_DAY + 2
-        P_return = pd.Series(np.zeros(output_length, dtype='float'),
-                             index=index)
-
-        # create kernel of unit hyetographs, one for each month
-        kernels = np.zeros(shape=(12, steps_per_two_days))
-
-        for month in np.arange(12, dtype=int):
-            # Calculating the unit hyetograph (kernels) with area of 1 mm
-
-            # Computing monthly constants
-            P_pk = 2. / dur[month]  # peak precipitation intensity
-            m = P_pk / (dur[month] / 2.)  # slope of precipitation rising limb
-            # time of storm start [minutes]
-            t_start = (t_pk[month] - (dur[month] / 2.))
-            # time of storm end [minutes]
-            t_end = (t_pk[month] + (dur[month] / 2.))
-            # time step of storm start
-            i_start = int(np.floor(t_start) + offset)
-            # time step of storm end
-            i_end = int(np.floor(t_end) + offset)
-
-            # Initializing timestep variables
-            i = i_start  # current timestep
-            t = i_start  # current time [minutes]
-            t_0 = t     # start time of current timestep [minutes]
-            # times of key corners of unit hyetograph
-            t_corners = [t_start, t_pk[month], t_end, t_next] + offset
-            c = 0   # current corner index
-            # end time of current timestep [minutes]
-            t_1 = min(t_0, t_corners[c])
-            area = 0  # area under curve for current timestep
-
-            # Looping through kernel timesteps
-            while i <= i_end:
-                P_0 = P_kernel(t_corners, m, t_0)
-                P_1 = P_kernel(t_corners, m, t_1)
-                area += 0.5 * (P_0 + P_1) * (t_1 - t_0)
-                if t_1 == t:  # end of timestep check
-                    kernels[month, i] = area
-                    area = 0
-                    t += 1
-                    i += 1
-                if t_1 == t_corners[c]:
-                    c += 1
-                t_0 = t_1
-                t_1 = min(t, t_corners[c])
-
-        # Loop through each rain day of the timeseries and apply the kernel for
-        # the appropriate month of year
-        rain_days = np.asarray(np.where(prec > 0))[0]
-        if len(rain_days) > 0:
-            for d in rain_days:
-                mon = month_of_year[d] - 1
-                if do_mix and t_min[d] < 0:
-                    i0 = d * steps_per_day
-                    i1 = i0 + steps_per_day
-                    P_return[i0:i1] += prec[d] * 1 / steps_per_day
-                elif d == 0:
-                    # beginning of kernel is clipped;
-                    # rescale so that clipped kernel sums to original total
-                    k0 = int(np.ceil(steps_per_day / 2))
-                    i1 = int(np.ceil(steps_per_day * 1.5))
-                    P_return[:i1] += (prec[d] / sum(kernels[mon, k0:])
-                                      * kernels[mon, k0:])
-                elif d == (n_days - 1):
-                    # end of kernel is clipped;
-                    # rescale so that clipped kernel sums to original total
-                    k1 = int(np.ceil(steps_per_day * 1.5))
-                    i0 = int(np.floor((d - 0.5) * steps_per_day))
-                    P_return[i0:] += (prec[d] / sum(kernels[mon, :k1])
-                                      * kernels[mon, :k1])
-                else:
-                    i0 = int(np.floor((d - 0.5) * steps_per_day))
-                    i1 = int(i0 + (2 * steps_per_day))
-                    P_return[i0:i1] += prec[d] * kernels[mon]
-        P_return = np.around(P_return, decimals=5)
-        return P_return.values
+    #def prec_TRIANGLE(prec: pd.Series, t_min: pd.Series,
+    #                  month_of_year: int, do_mix: bool, params: dict):
+    def prec_TRIANGLE(daily_prec, t_min, prec_peak, prec_dur, ts, do_mix):
+        '''Triangular disaggregation'''
+        ts_per_day = int(cnst.HOURS_PER_DAY * cnst.MIN_PER_HOUR)
+        n_days = len(daily_prec)
+        # Scale to desired timestep
+        #prec_peak /= ts
+        #prec_dur /= ts
+        disagg_prec = np.zeros(int(n_days*ts_per_day))
+
+        # Loop over days
+        for i, (t, P) in enumerate(daily_prec.iteritems()):
+            times_day = np.arange(ts_per_day)
+            prec_day = np.zeros(ts_per_day)
+
+            # Look up climatology
+            t_pk = prec_peak[t]
+            t_dur = prec_dur[t]
+
+            # Rising and falling time
+            t_start = t_pk - 0.5 * t_dur
+            t_stop = t_pk + 0.5 * t_dur
+            t_plus = times_day[
+                    np.logical_and(times_day <= t_pk, times_day >= t_start)]
+            t_minus = times_day[
+                    np.logical_and(times_day >= t_pk, times_day <= t_stop)]
+
+            # Begin with relative intensity
+            prec_day[t_plus] = np.linspace(0, 1.0, len(t_plus))
+            prec_day[t_minus] = np.linspace(1.0, 0, len(t_minus))
+
+            # Scale to input precipitation
+            prec_day = (P / np.sum(prec_day)) * prec_day
+            disagg_prec[i*ts_per_day:(i+1)*ts_per_day] = prec_day
+        return disagg_prec
 
     def prec_UNIFORM(prec: pd.Series, *args):
-        return np.repeat(prec / cnst.MIN_PER_DAY, cnst.MIN_PER_DAY)
+        return np.repeat(prec.values / cnst.MIN_PER_DAY, cnst.MIN_PER_DAY)
 
     prec_function = {
         'UNIFORM': prec_UNIFORM,
         'TRIANGLE': prec_TRIANGLE,
         'MIX': prec_TRIANGLE,
     }
-    if params['prec_type'].upper() == 'MIX':
-        do_mix = True
+    if params['prec_type'].upper() in ['TRIANGLE', 'MIX']:
+        if params['prec_type'].upper() == 'MIX':
+            do_mix = True
+        else:
+            do_mix = False
+        P_return = prec_TRIANGLE(prec, t_min, t_pk, dur, ts, do_mix)
     else:
-        do_mix = False
-
-    P_return = prec_function[params['prec_type'].upper()](prec, t_min,
-                                                          month_of_year,
-                                                          do_mix, params)
+        P_return = prec_UNIFORM(prec)
 
     if params['utc_offset']:
         offset = int(((params["lon"] / cnst.DEG_PER_REV) * cnst.MIN_PER_DAY))
diff --git a/metsim/io.py b/metsim/io.py
@@ -254,7 +254,7 @@ def read_ascii(data_handle, domain=None, is_worker=False,
                var_dict=None) -> xr.Dataset:
     """Read in an ascii forcing file"""
     dates = date_range(start, stop, calendar=calendar)
-    names = var_dict.keys()
+    names = list(var_dict.keys())
     ds = pd.read_csv(data_handle, header=None, delim_whitespace=True,
                      names=names).head(len(dates))
     ds.index = dates
diff --git a/metsim/metsim.py b/metsim/metsim.py
@@ -152,6 +152,7 @@ class MetSim(object):
         "utc_offset": False,
         "lw_cloud": 'cloud_deardorff',
         "lw_type": 'prata',
+        "prec_type": 'uniform',
         "tdew_tol": 1e-6,
         "tmax_daylength_fraction": 0.67,
         "rain_scalar": 0.75,
@@ -235,7 +236,6 @@ def _update_unit_attrs(self, out_vars):
                 v['units'] = available_outputs[k]['units']
 
     def _validate_force_times(self, force_times):
-
         for p, i in [('start', 0), ('stop', -1)]:
             # infer times from force_times
             if isinstance(self.params[p], str):
@@ -270,6 +270,14 @@ def _validate_force_times(self, force_times):
                              'long_name': 'local time at grid location',
                              'standard_name': 'local_time'}
 
+    def convert_monthly_param(self, name):
+        self.met_data[name] = self.met_data['prec'].copy()
+        months = self.met_data['time'].dt.month
+        for m in range(12):
+            param = self.domain[name].sel(month=m)
+            locations = {'time': self.met_data['time'].isel(time=months == m)}
+            self.met_data[name].loc[locations] = param
+
     @property
     def domain(self):
         if self._domain is None:
@@ -458,14 +466,12 @@ def run_slice(self):
         params['tday_coef'] = float(params['tday_coef'])
         params['tmax_daylength_fraction'] = float(params['tmax_daylength_fraction'])
         params['lapse_rate'] = float(params['lapse_rate'])
+        if self.params['prec_type'].upper() in ['TRIANGLE', 'MIX']:
+            self.convert_monthly_param('dur')
+            self.convert_monthly_param('t_pk')
         for index, mask_val in np.ndenumerate(self.domain['mask'].values):
             if mask_val > 0:
                 locs = {d: i for d, i in zip(self.domain['mask'].dims, index)}
-                if self.params['prec_type'].upper() in ['TRIANGLE', 'MIX']:
-                    # add variables for triangle precipitation disgregation
-                    # method to parameters
-                    params['dur'], params['t_pk'] = (
-                        add_prec_tri_vars(self.domain.isel(**locs)))
             else:
                 continue
 
@@ -826,52 +832,3 @@ def __exit__(self, *args):
     @property
     def locked(self):
         return False
-
-
-def add_prec_tri_vars(domain):
-    """
-    Check that variables for triangle precipitation method exist and have
-    values that are within allowable ranges. Return these variables.
-
-    Parameters
-    ----------
-    domain:
-        Dataset of domain variables for given location.
-
-    Returns
-    -------
-    dur
-        Array of climatological monthly storm durations. [minutes]
-    t_pk
-        Array of climatological monthly times to storm peaks. [minutes]
-    """
-    # Check that variables exist
-    try:
-        dur = domain['dur']
-    except Exception as e:
-        raise e("Storm duration and time to peak values are "
-                "required in the domain file for the triangle "
-                "preciptation disagregation method.")
-    try:
-        t_pk = domain['t_pk']
-    except Exception as e:
-        raise e("Storm duration and time to peak values are "
-                "required in the domain file for the triangle "
-                "preciptation disagregation method.")
-
-    # Check that variable values are within allowable ranges.
-    day_length = cnst.MIN_PER_HOUR * cnst.HOURS_PER_DAY
-    dur_zero_test = dur <= 0
-    dur_day_test = dur > day_length
-    if dur_zero_test.any() or dur_day_test.any():
-        raise ValueError('Storm duration must be greater than 0 and less than',
-                         day_length, '(i.e. the day length in minutes)')
-
-    t_pk_zero_test = t_pk < 0
-    t_pk_day_test = t_pk > day_length
-    if t_pk_zero_test.any() or t_pk_day_test.any():
-        raise ValueError('Storm time to peak must be greater than or equal to '
-                         '0, and less than', day_length,
-                         '(i.e. the end of a day)')
-
-    return dur, t_pk
