Commit 941f876e authored by Cynthia Crowley's avatar Cynthia Crowley

Integrate forcing variables into files separate from LSM results

parent 93c0986f
......@@ -141,6 +141,11 @@ class GLDAS20_NoahConfig(ConfigBase):
assert False
@classmethod
def forcing_integrated_vars(cls, basis=None):
return {}
def result_postprocess_steps(self, yearmon=None, target=None, member=None):
input_file = os.path.join(self._source,
'GLDAS_NOAH025_M.A{}.020.nc4'.format(yearmon))
......
......@@ -60,7 +60,7 @@ def create_forcing_file(workspace: DefaultWorkspace,
'Pr:units',
'Pr:standard_name'
])),
output=workspace.forcing(yearmon=yearmon, target=target, member=member)
output=workspace.forcing(yearmon=yearmon, target=target, member=member, window=1)
)
]
......@@ -204,7 +204,7 @@ def compute_return_periods(workspace: DefaultWorkspace, *,
obs=[
read_vars(workspace.results(**args), *result_vars)
if result_vars else None,
read_vars(workspace.forcing(yearmon=yearmon, target=target, member=member), *forcing_vars)
read_vars(workspace.forcing(yearmon=yearmon, target=target, window=window, member=member), *forcing_vars)
if forcing_vars else None
],
rp=workspace.return_period(**args),
......@@ -366,11 +366,11 @@ def fit_var(config: ConfigBase,
param_to_read = param
# Don't look for Pr_sum and T_ave in forcing directory.
if param in config.forcing_rp_vars() and stat is None:
assert window == 1
if param in config.forcing_rp_vars() and stat is None:
assert window == 1
assert basis is None
infile = config.workspace().forcing(yearmon=input_range)
infile = config.workspace().forcing(yearmon=input_range, window=window)
else:
infile = config.workspace().results(yearmon=input_range, window=window, basis=basis)
......
......@@ -157,7 +157,7 @@ class ConfigBase(metaclass=abc.ABCMeta):
'RO_mm',
'Sa',
'Sm',
'Snowpack',
#'Snowpack',
'Ws'
]
......@@ -169,6 +169,21 @@ class ConfigBase(metaclass=abc.ABCMeta):
assert False
@classmethod
def state_vars(*, basis: Optional[Basis]=None) -> List[str]:
"""
Provides a list of state variables for which return periods should be calculated
"""
if not basis:
return['Snowpack']
if basis==Basis.BASIN:
return['Bt_RO']
assert False
@classmethod
def lsm_integrated_vars(cls, basis: Optional[Basis]=None) -> Dict[str, List[str]]:
"""
......@@ -232,6 +247,22 @@ class ConfigBase(metaclass=abc.ABCMeta):
assert False
@classmethod
def forcing_integrated_stats(cls, basis: Optional[Basis]=None) -> Dict[str, List[str]]:
"""
Provides a dictionary whose keys are stat names and whose values are a list of variables
two which that stat should be applied. It can be thought of as the inverse of forcing_integrated_vars()
"""
integrated_stats = {}
for var, varstats in cls.forcing_integrated_vars(basis=basis).items():
for stat in varstats:
if stat not in integrated_stats:
integrated_stats[stat] = []
integrated_stats[stat].append(var)
return integrated_stats
@classmethod
def forcing_integrated_var_names(cls, basis: Optional[Basis]=None) -> List[str]:
"""
......
......@@ -60,8 +60,8 @@ def monthly_observed(config: Config, yearmon: str, meta_steps: Dict[str, Step])
# Compute return periods
for window in [1] + config.integration_windows():
steps += compute_return_periods(config.workspace(),
result_vars=config.lsm_rp_vars() + config.forcing_rp_vars() if window == 1 else config.lsm_integrated_var_names() + config.forcing_integrated_var_names(),
forcing_vars=config.forcing_rp_vars() if window==1 else None,
result_vars=config.lsm_rp_vars() + config.forcing_rp_vars() if window == 1 else config.lsm_integrated_var_names(),
forcing_vars=config.forcing_rp_vars() if window==1 else config.forcing_integrated_var_names(),
yearmon=yearmon,
window=window)
......
......@@ -422,7 +422,11 @@ class DefaultWorkspace:
return self.make_path('state', sector=sector, yearmon=yearmon, member=member, target=target, window=None, method=method)
def forcing(self, *, yearmon: str, member: Optional[str]=None, target: Optional[str]=None) -> str:
def forcing(self, *,
yearmon: str,
window: int,
member: Optional[str]=None,
target: Optional[str]=None) -> str:
return self.make_path('forcing', yearmon=yearmon, member=member, target=target, window=None)
def results(self, *,
......
......@@ -61,8 +61,10 @@ def spinup(config, meta_steps):
# Time-integrate the variables
for window in config.integration_windows():
steps += time_integrate_forcing(config, window)
steps += time_integrate_results(config, window)
# Compute monthly fits (and then anomalies) over the fit period
for param in config.lsm_rp_vars() + config.forcing_rp_vars():
for month in all_months:
......@@ -209,7 +211,7 @@ def run_lsm_from_final_norm_state(config: Config) -> List[Step]:
# allow restarting in case of failure. But the runtime becomes dominated by the
# R startup and I/O, and takes about 5 seconds / iteration instead of 1 second /iteration.
run_lsm = wsim_lsm(
forcing=[config.workspace().forcing(yearmon=date_range(config.historical_yearmons()))],
forcing=[config.workspace().forcing(yearmon=date_range(config.historical_yearmons()), window=1)],
state=config.workspace().spinup_state(yearmon=initial_yearmon),
elevation=config.static_data().elevation(),
flowdir=config.static_data().flowdir(),
......@@ -265,7 +267,7 @@ def run_lsm_from_mean_spinup_state(config: Config) -> List[Step]:
run_lsm = wsim_lsm(
comment="LSM run from mean spinup state",
forcing=[config.workspace().forcing(yearmon=date_range(config.historical_yearmons()))],
forcing=[config.workspace().forcing(yearmon=date_range(config.historical_yearmons()), window=1)],
state=config.workspace().state(yearmon=first_timestep),
elevation=config.static_data().elevation(),
flowdir=config.static_data().flowdir(),
......@@ -288,6 +290,38 @@ def run_lsm_from_mean_spinup_state(config: Config) -> List[Step]:
]
def time_integrate_forcing(config:Config, window: int, *, basis: Optional[Basis]=None) -> List[Step]:
"""
Integrate forcing variables over the given time window
"""
yearmons_in = config.historical_yearmons()
yearmons_out = yearmons_in[window-1:]
integrate = wsim_integrate(
inputs= [read_vars(config.workspace().forcing(window=1, yearmon=date_range(yearmons_in)),
*config.forcing_integrated_vars(basis=basis).keys())
],
window=window,
stats=[stat + '::' + ','.join(varname) for stat, varname in config.forcing_integrated_stats(basis=basis).items()],
attrs=[attrs.integration_window(var='*', months=window)],
output=config.workspace().forcing(yearmon=date_range(yearmons_out),
window=window)
)
tag_name = config.workspace().tag('{}spinup_{}mo_forcing'.format((basis.value + '_' if basis else ''), window))
tag_steps = create_tag(name=tag_name, dependencies=integrate.targets)
integrate.replace_targets_with_tag_file(tag_name)
integrate.replace_dependencies(
config.workspace().tag('{}spinup_1mo_forcing'.format((basis.value + '_') if basis else '')))
return [
integrate,
*tag_steps
]
def time_integrate_results(config: Config, window: int, *, basis: Optional[Basis]=None) -> List[Step]:
"""
Integrate specified LSM results (and any included forcing variables) over the given time window
......@@ -297,10 +331,10 @@ def time_integrate_results(config: Config, window: int, *, basis: Optional[Basis
integrate = wsim_integrate(
inputs=[read_vars(config.workspace().results(window=1, yearmon=date_range(yearmons_in), basis=basis),
*{**config.lsm_integrated_vars(basis=basis), **config.forcing_integrated_vars(basis=basis)}.keys())
*config.lsm_integrated_vars(basis=basis).keys())
],
window=window,
stats=[stat + '::' + ','.join(varname) for stat, varname in config.all_integrated_stats(basis=basis).items()],
stats=[stat + '::' + ','.join(varname) for stat, varname in config.lsm_integrated_stats(basis=basis).items()],
attrs=[attrs.integration_window(var='*', months=window)],
output=config.workspace().results(yearmon=date_range(yearmons_out),
window=window,
......
......@@ -20,15 +20,12 @@
#' @param PET potential evapotranspiration [mm]
#' @param PETmE potential minus actual evapotranspiration [mm]
#' @param P_net net precipitation [mm]
#' @param Pr precipitation [mm]
#' @param RO_m3 runoff (taking account of detention) [m^3]
#' @param RO_mm runoff (taking account of detention) [mm]
#' @param Runoff_m3 runoff (not taking account of detention) [m^3]
#' @param Runoff_mm runoff (not taking account of detention) [mm]
#' @param Sa snow accumulation [mm]
#' @param Sm snowmelt [mm]
#' @param Snowpack snow water equivalent [mm]
#' @param T temperature [degrees C]
#' @param Ws average soil moisture [mm]
#' @param dWdt change in soil moisture [mm]
#' @param extent spatial extent of input matrices \code{(xmin, xmax, ymin, ymax)}
......@@ -44,15 +41,15 @@ make_results <- function(
PET,
PETmE,
P_net,
Pr,
# Pr,
RO_m3,
RO_mm,
Runoff_m3,
Runoff_mm,
Sa,
Sm,
Snowpack,
T,
# Snowpack,
# T,
Ws,
dWdt,
extent
......@@ -66,15 +63,15 @@ make_results <- function(
PET= PET,
PETmE= PETmE,
P_net= P_net,
Pr=Pr,
# Pr=Pr,
RO_m3= RO_m3,
RO_mm= RO_mm,
Runoff_m3= Runoff_m3,
Runoff_mm= Runoff_mm,
Sa= Sa,
Sm= Sm,
Snowpack=Snowpack,
T= T,
# Snowpack=Snowpack,
# T= T,
Ws= Ws,
dWdt= dWdt
)
......
......@@ -96,15 +96,15 @@ run <- function(static, state, forcing) {
P_net= P,
PET= E0,
PETmE= E0 - E,
Pr= forcing$Pr,
#Pr= forcing$Pr,
RO_mm= revised_runoff,
RO_m3= revised_runoff*area_m2/1000,
Runoff_mm= R,
Runoff_m3= R*area_m2/1000,
Sa= Sa,
Sm= ifelse(is.na(Sa), NA, Sm),
Snowpack= state$Snowpack,
T= forcing$T,
#Snowpack= state$Snowpack,
#T= forcing$T,
Ws= Ws_ave
)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment