Merge pull request #244 from UW-Hydro/develop

Update for 2.4

docs/configuration.rst

@@ -38,6 +38,11 @@ the ``forcing`` entry.  Can be one of the following: ``ascii``, ``binary``,
 
 **Optional Variables**
 
+``method ::str``: The method to use for estimation of meteorological quantities.
+This can be either ``mtclim`` to estimate missing variables or ``passthrough`` if
+some of the meteorological variables have already been estimated (for example, by
+DayMet, PRISM, or GridMET). Defaults to ``mtclim``.
+
 ``out_prefix :: str``: The output file base name. Defaults to ``forcing``.
 
 ``out_precision :: str``: Precision to use when writing output.  Defaults to

docs/whats-new.rst

@@ -3,6 +3,18 @@
 What's New
 ==========
 
+.. _whats-new.2.4.0:
+
+v2.4.0
+------
+Enchancements
+~~~~~~~~~~~~~
+- Allow for passing already estimated met variables
+  (such as shortwave and/or longwave radiation)
+  through to the disaggregation routines. This
+  functionality can be accessed by setting the
+  ``method`` to ``passthrough`` in the configuration
+
 .. _whats-new.2.3.0:
 
 v2.3.3

metsim/disaggregate.py

@@ -105,9 +105,13 @@ def disaggregate(df_daily: pd.DataFrame, params: dict,
 
     df_disagg['tskc'] = tskc(df_daily['tskc'].values, ts, params)
 
+    if 'longwave' in df_daily:
+        daily_lw = df_daily['longwave']
+    else:
+        daily_lw = None
     df_disagg['longwave'] = longwave(
         df_disagg['temp'].values, df_disagg['vapor_pressure'].values,
-        df_disagg['tskc'].values, params)
+        df_disagg['tskc'].values, params, daily_lw)
     df_disagg['prec'] = prec(df_daily['prec'], df_daily['t_min'], ts, params,
                              df_daily.get('t_pk'), df_daily.get('dur'))
 
@@ -478,7 +482,7 @@ def vapor_pressure(vp_daily: np.array, temp: np.array, t_t_min: np.array,
 
 
 def longwave(air_temp: np.array, vapor_pressure: np.array,
-             tskc: np.array, params: dict) -> np.array:
+             tskc: np.array, params: dict, daily_lw=None) -> np.array:
     """
     Calculate longwave. This calculation can be performed
     using a variety of parameterizations for both the
@@ -569,6 +573,13 @@ def longwave(air_temp: np.array, vapor_pressure: np.array,
     emiss_func = cloud_calc[params['lw_cloud'].upper()]
     emissivity = emiss_func(emissivity_clear, tskc)
     lwrad = emissivity * cnst.STEFAN_B * np.power(air_temp, 4)
+    if daily_lw is not  None:
+        ts = int(params['time_step'])
+        ts_per_day = int(cnst.HOURS_PER_DAY * cnst.MIN_PER_HOUR / ts)
+        factor = np.mean(lwrad.reshape(-1, ts_per_day), axis=1).flatten()
+        factor = daily_lw.values / factor
+        factor = np.repeat(factor, ts_per_day)
+        lwrad *= factor
     return lwrad
 
 
@@ -623,14 +634,17 @@ def shortwave(sw_rad: np.array, daylength: np.array, day_of_year: np.array,
     """
     ts = int(params['time_step'])
     ts_hourly = float(ts) / cnst.MIN_PER_HOUR
-    tmp_rad = (sw_rad * daylength) / (cnst.SEC_PER_HOUR * ts_hourly)
+    if params['method'] == 'mtclim' or params.get('sw_averaging', '') == 'daylight':
+        tmp_rad = (sw_rad * daylength) / (cnst.SEC_PER_HOUR * ts_hourly)
+    else:
+        tmp_rad = sw_rad * 24
     n_days = len(tmp_rad)
     ts_per_day = int(cnst.HOURS_PER_DAY * cnst.MIN_PER_HOUR / ts)
     disaggrad = np.zeros(int(n_days * ts_per_day))
     rad_fract_per_day = int(cnst.SEC_PER_DAY / cnst.SW_RAD_DT)
     tmp_rad = np.repeat(tmp_rad, rad_fract_per_day)
     if params['utc_offset']:
-        utc_offset = int((params['lon']  / cnst.DEG_PER_REV) * rad_fract_per_day)
+        utc_offset = int((params['lon'] / cnst.DEG_PER_REV) * rad_fract_per_day)
         tiny_rad_fract = np.roll(tiny_rad_fract.flatten(), -utc_offset)
         tmp_rad = np.roll(tmp_rad.flatten(), -utc_offset)
         tiny_rad_fract = tiny_rad_fract.flatten()

metsim/methods/passthrough.py

@@ -0,0 +1,46 @@
+"""
+Passthrough
+"""
+# Meteorology Simulator
+# Copyright (C) 2017  The Computational Hydrology Group, Department of Civil
+# and Environmental Engineering, University of Washington.
+
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+import numpy as np
+import pandas as pd
+
+import metsim.constants as cnst
+from metsim.physics import atm_pres, calc_pet, svp
+from metsim.methods.mtclim import t_day, tfmax, tskc, pet, tdew, vapor_pressure
+
+def run(df, params):
+    assert 'shortwave' in df
+
+    if 't_day' not in df:
+        df['t_day'] = t_day(df['t_min'].values, df['t_max'].values, params)
+    if 'tfmax' not in df:
+        df['tfmax'] = tfmax(df['dtr'].values, df['smoothed_dtr'].values,
+                            df['prec'].values, params)
+    if 'tskc' not in df:
+        df['tskc'] = tskc(df['tfmax'].values, params)
+    if 'pet' not in df:
+        df['pet'] = pet(df['shortwave'].values, df['t_day'].values,
+                        df['daylength'].values, params)
+    if 'tdew' not in df:
+        df['tdew'] = tdew(df['pet'].values, df['t_min'].values,
+                          df['seasonal_prec'].values, df['dtr'].values)
+    if 'vapor_pressure' not in df:
+        df['vapor_pressure'] = vapor_pressure(df['tdew'].values)
+    return df

metsim/metsim.py

@@ -53,7 +53,7 @@
 from metsim import io
 from metsim.datetime import date_range
 from metsim.disaggregate import disaggregate
-from metsim.methods import mtclim
+from metsim.methods import mtclim, passthrough
 from metsim.physics import solar_geom
 from metsim.units import converters
 
@@ -104,6 +104,9 @@
          'rel_humid': {'units': '%', 'long_name': 'relative humidity',
                        'standard_name': 'relative_humidity',
                        'missing_value': np.nan, 'fill_value': np.nan},
+         'daylength': {'units': 's', 'long_name': 'daylength',
+                       'standard_name': 'length of day',
+                       'missing_value': np.nan, 'fill_value': np.nan},
          'spec_humid': {'units': 'g g-1', 'long_name': 'specific humidity',
                         'standard_name': 'specific_humidity',
                         'missing_value': np.nan, 'fill_value': np.nan},
@@ -131,7 +134,7 @@ class MetSim(object):
     """
 
     # Class variables
-    methods = {'mtclim': mtclim}
+    methods = {'mtclim': mtclim, 'passthrough': passthrough}
     params = {
         "period_ending": False,
         "is_worker": False,
@@ -474,7 +477,6 @@ def run_slice(self):
                 locs = {d: i for d, i in zip(self.domain['mask'].dims, index)}
             else:
                 continue
-
             df, state = wrap_run_cell(self.method.run, params,
                                       self.met_data.isel(**locs),
                                       self.state.isel(**locs),
@@ -654,7 +656,7 @@ def _validate_setup(self):
 
         # Check output variables are valid
         daily_out_vars = ['t_min', 't_max', 't_day', 'prec', 'vapor_pressure',
-                          'shortwave', 'tskc', 'pet', 'wind']
+                          'shortwave', 'tskc', 'pet', 'wind', 'daylength']
         out_var_check = ['temp', 'prec', 'shortwave', 'vapor_pressure',
                          'air_pressure', 'rel_humid', 'spec_humid',
                          'longwave', 'tskc', 'wind']
@@ -778,7 +780,13 @@ def wrap_run_cell(func: callable, params: dict,
         # If we're outputting daily values, we dont' need to
         # change the output dates - see inside of `if` condition
         # above for more explanation
-        new_times = out_times
+        start = out_times.values[0]
+        stop = (out_times.values[-1] + pd.Timedelta('1 days') -
+                pd.Timedelta("{} minutes".format(params['time_step'])))
+        new_times = date_range(
+            start, stop, freq='{}T'.format(params['time_step']),
+            calendar=params['calendar'])
+        df_base.index = new_times
         df_complete = df_base
 
     # Cut the returned data down to the correct time index

metsim/tests/test_metsim.py

@@ -277,6 +277,68 @@ def test_variable_rename():
     assert 'SWRadAtm' in ds.variables
 
 
+def test_passthrough():
+    """Test to make sure passing through previously estimated
+       variables doesn't alter the values from disaggregation"""
+    loc = data_locations['binary']
+    data_files = [os.path.join(loc, f) for f in os.listdir(loc)]
+    out_dir = '.'
+    params = {'start': dates['binary'][0],
+              'stop': dates['binary'][1],
+              'forcing_fmt': 'binary',
+              'domain_fmt': 'netcdf',
+              'state_fmt': 'netcdf',
+              'domain': './metsim/data/stehekin.nc',
+              'state': './metsim/data/state_vic.nc',
+              'forcing': data_files,
+              'method': 'mtclim',
+              'scheduler': 'threading',
+              'time_step': "60",
+              'out_dir': out_dir,
+              'out_state': os.path.join(out_dir, 'state.nc'),
+              'out_vars': {'prec': {'out_name': 'prec'},
+                           'temp': {'out_name': 'temp'},
+                           'longwave': {'out_name': 'longwave'},
+                           'vapor_pressure': {'out_name': 'vapor_pressure'},
+                           'shortwave': {'out_name': 'shortwave'}},
+              'forcing_vars': in_vars_section['binary'],
+              'domain_vars': domain_section['binary']}
+
+    # Second run will be to use mtclim and hourly disagg
+    params1 = dict()
+    params1.update(params)
+    params1['out_prefix'] = 'mtclim'
+    ms1 = MetSim(params1)
+    ms1.run()
+    with ms1.open_output() as ds:
+        mtclim_ds = ds.load()
+
+    # Third run will be to use passthrough and hourly disagg
+    # with input data from teh first run
+    params2 = dict()
+    params2.update(params)
+    params2['method'] = 'passthrough'
+    params2['out_prefix'] = 'passthrough'
+    params2['forcing_vars'] = OrderedDict(
+        prec='prec', t_max='t_max', t_min='t_min',
+        wind='wind', shortwave='shortwave', vapor_pressure='vapor_pressure')
+    params2['forcing_fmt'] = 'netcdf'
+    params2['forcing'] = './metsim/data/passthrough.nc'
+    ms2 = MetSim(params2)
+    ms2.run()
+    with ms2.open_output() as ds:
+        passthrough_ds = ds.load()
+
+    tol = 1e-4
+    assert np.allclose(passthrough_ds['shortwave'].mean(),
+                       mtclim_ds['shortwave'].mean(), atol=tol)
+    assert np.allclose(passthrough_ds['vapor_pressure'].mean(),
+                       mtclim_ds['vapor_pressure'].mean(), atol=tol)
+    assert np.allclose(passthrough_ds['longwave'].mean(),
+                       mtclim_ds['longwave'].mean(), atol=tol)
+
+
+
 def test_unit_conversion():
     """Tests to make sure that variable renaming works"""
     loc = data_locations['binary']

metsim/units.py

@@ -52,5 +52,8 @@
             },
         'wind': {
             'm s-1': lambda x, ts: x,
+            },
+        'daylength': {
+            's': lambda x, ts: x,
             }
         }