diff --git a/docs/configuration.rst b/docs/configuration.rst index ccd46b2..7adc720 100644 --- a/docs/configuration.rst +++ b/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 diff --git a/docs/whats-new.rst b/docs/whats-new.rst index 7e62d7b..b79e9c9 100644 --- a/docs/whats-new.rst +++ b/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 diff --git a/metsim/data/passthrough.nc b/metsim/data/passthrough.nc new file mode 100644 index 0000000..5405bae Binary files /dev/null and b/metsim/data/passthrough.nc differ diff --git a/metsim/disaggregate.py b/metsim/disaggregate.py index 3d453f3..5ecb2ca 100644 --- a/metsim/disaggregate.py +++ b/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() diff --git a/metsim/methods/passthrough.py b/metsim/methods/passthrough.py new file mode 100644 index 0000000..9e99a57 --- /dev/null +++ b/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 . + +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 diff --git a/metsim/metsim.py b/metsim/metsim.py index c3a9f85..8224ae3 100644 --- a/metsim/metsim.py +++ b/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 diff --git a/metsim/tests/test_metsim.py b/metsim/tests/test_metsim.py index 2dab1fa..0868bfd 100755 --- a/metsim/tests/test_metsim.py +++ b/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'] diff --git a/metsim/units.py b/metsim/units.py index d182ce9..c784f5d 100644 --- a/metsim/units.py +++ b/metsim/units.py @@ -52,5 +52,8 @@ }, 'wind': { 'm s-1': lambda x, ts: x, + }, + 'daylength': { + 's': lambda x, ts: x, } }