diff --git a/.gitignore b/.gitignore index 4c7f1c3f..3c7b339c 100644 --- a/.gitignore +++ b/.gitignore @@ -62,6 +62,7 @@ generated/ # example output data and gis data examples/0*/ examples/snow_errors +examples/compare_pws_prms examples/runoff_errors examples/pynhm_gis examples/pywatershed_gis diff --git a/autotest/test_control_read.py b/autotest/test_control_read.py index bf1d8e06..fbc9a962 100644 --- a/autotest/test_control_read.py +++ b/autotest/test_control_read.py @@ -1,7 +1,9 @@ +import pathlib as pl + import numpy as np import pytest -from pywatershed.utils import ControlVariables +from pywatershed.utils import ControlVariables, compare_control_files from utils import assert_or_print @@ -19,11 +21,9 @@ def control_keys(): def test_control_read(domain, control_keys): control_file = domain["control_file"] print(f"parsing...'{control_file}'") - control = ControlVariables.load(control_file) - - # check control data answers = domain["test_ans"]["control_read"] + for key, value in answers.items(): if key in ( "start_time", @@ -32,13 +32,23 @@ def test_control_read(domain, control_keys): answers[key] = np.datetime64(value) elif key in ("initial_deltat",): answers[key] = np.timedelta64(int(value), "h") + results = { key: val for key, val in control.control.items() if key in answers.keys() } assert_or_print(results, answers, print_ans=domain["print_ans"]) - print(f"success parsing...'{control_file}'") + return + +def test_control_compare(): + common_dir = pl.Path("../test_data/common") + n_diffs = compare_control_files( + common_dir / "control.single_hru", + common_dir / "control.multi_hru", + silent=True, + ) + assert n_diffs == 5 return diff --git a/autotest/test_model.py b/autotest/test_model.py index 8e241429..64a9b62d 100644 --- a/autotest/test_model.py +++ b/autotest/test_model.py @@ -305,8 +305,8 @@ def test_model(domain, model_args, tmp_path): 9: { "PRMSChannel": { "seg_outflow": { - "drb_2yr": 1517.232887980279, - "hru_1": 13.696918669514927, + "drb_2yr": 1553.1874672413599, + "hru_1": 13.696710376067216, "ucb_2yr": 1694.5697712423928, }, }, @@ -314,9 +314,9 @@ def test_model(domain, model_args, tmp_path): 99: { "PRMSChannel": { "seg_outflow": { - "drb_2yr": 2350.499659332901, - "hru_1": 22.874414994530095, - "ucb_2yr": 733.2293013532435, + "drb_2yr": 2362.7940777644653, + "hru_1": 22.877787915898086, + "ucb_2yr": 733.0192586668293, }, }, }, diff --git a/autotest/test_prms_atmosphere.py b/autotest/test_prms_atmosphere.py index 5874c7cc..35b1aa6b 100644 --- a/autotest/test_prms_atmosphere.py +++ b/autotest/test_prms_atmosphere.py @@ -1,6 +1,3 @@ -from warnings import warn - -import numpy as np import pytest from pywatershed.atmosphere.prms_atmosphere import PRMSAtmosphere @@ -8,28 +5,15 @@ from pywatershed.base.control import Control from pywatershed.base.parameters import Parameters from pywatershed.parameters import PrmsParameters +from utils_compare import compare_in_memory, compare_netcdfs -params = ["params_sep", "params_one"] - - -# @pytest.fixture(scope="function", params=params) -# def control(domain, request): -# if request.param == "params_one": -# params = PrmsParameters.load(domain["param_file"]) -# dis = None - -# else: -# # channel needs both hru and seg dis files -# dis_hru_file = domain["dir"] / "parameters_dis_hru.nc" -# dis_data = Parameters.merge( -# Parameters.from_netcdf(dis_hru_file, encoding=False), -# ) -# dis = {"dis_hru": dis_data} +# compare in memory (faster) or full output files? or both! +do_compare_output_files = False +do_compare_in_memory = True +rtol = 1.0e-5 +atol = 1.0e-5 # why is this relatively low accuracy? -# param_file = domain["dir"] / "parameters_PRMSAtmosphere.nc" -# params = {"PRMSAtmosphere": PrmsParameters.from_netcdf(param_file)} - -# return Control.load(domain["control_file"], params=params, dis=dis) +params = ["params_sep", "params_one"] @pytest.fixture(scope="function") @@ -55,31 +39,11 @@ def parameters(domain, request): def test_compare_prms(domain, control, discretization, parameters, tmp_path): + comparison_var_names = PRMSAtmosphere.get_variables() + output_dir = domain["prms_output_dir"] cbh_dir = domain["cbh_inputs"]["prcp"].parent.resolve() - # get the answer data - comparison_var_names = [ - "tmaxf", - "tminf", - "hru_ppt", - "hru_rain", - "hru_snow", - "swrad", - "potet", - "transp_on", - "tmaxc", - "tavgc", - "tminc", - "prmx", - "pptmix", - "orad_hru", - ] - ans = {} - for key in comparison_var_names: - nc_pth = output_dir / f"{key}.nc" - ans[key] = adapter_factory(nc_pth, variable_name=key, control=control) - input_variables = {} for key in PRMSAtmosphere.get_inputs(): dir = "" @@ -96,53 +60,35 @@ def test_compare_prms(domain, control, discretization, parameters, tmp_path): netcdf_output_dir=tmp_path, ) - all_success = True - for istep in range(control.n_times): - control.advance() - atm.advance() - atm.calculate(1.0) - - # compare along the way - for key, val in ans.items(): - val.advance() - - for key in ans.keys(): - a1 = ans[key].current - a2 = atm[key].current - - tol = 1e-5 - if key == "swrad": - tol = 5e-4 - warn(f"using tol = {tol} for variable {key}") - if key == "tavgc": - tol = 1e-5 - warn(f"using tol = {tol} for variable {key}") - - success_a = np.allclose(a2, a1, atol=tol, rtol=0.00) - success_r = np.allclose(a2, a1, atol=0.00, rtol=tol) - success = False - if (not success_a) and (not success_r): - diff = a2 - a1 - diffratio = abs(diff / a2) - if (diffratio < 1e-6).all(): - success = True - continue - all_success = False - diffmin = diff.min() - diffmax = diff.max() - abs_diff = abs(diff) - absdiffmax = abs_diff.max() - wh_absdiffmax = np.where(abs_diff)[0] - print(f"time step {istep}") - print(f"output variable {key}") - print(f"prms {a1.min()} {a1.max()}") - print(f"pywatershed {a2.min()} {a2.max()}") - print(f"diff {diffmin} {diffmax}") - print(f"absdiffmax {absdiffmax}") - print(f"wh_absdiffmax {wh_absdiffmax}") - assert success - - atm.finalize() - - if not all_success: - raise Exception("pywatershed results do not match prms results") + if do_compare_in_memory: + answers = {} + for var in comparison_var_names: + var_pth = output_dir / f"{var}.nc" + answers[var] = adapter_factory( + var_pth, variable_name=var, control=control + ) + + # check the advance/calculate the state + tmaxf_id = id(atm.tmaxf) + + for ii in range(control.n_times): + control.advance() + atm.advance() + if ii == 0: + atm.output() + atm.calculate(1.0) + + compare_in_memory(atm, answers, atol=atol, rtol=rtol) + assert id(atm.tmaxf) == tmaxf_id + + if do_compare_output_files: + compare_netcdfs( + comparison_var_names, + tmp_path, + output_dir, + atol=atol, + rtol=rtol, + print_var_max_errs=False, + ) + + return diff --git a/autotest/test_prms_canopy.py b/autotest/test_prms_canopy.py index f87bd2e1..f022f03e 100644 --- a/autotest/test_prms_canopy.py +++ b/autotest/test_prms_canopy.py @@ -1,12 +1,17 @@ import pathlib as pl -import numpy as np import pytest from pywatershed.base.adapter import adapter_factory from pywatershed.base.control import Control from pywatershed.hydrology.prms_canopy import PRMSCanopy, has_prmscanopy_f from pywatershed.parameters import Parameters, PrmsParameters +from utils_compare import compare_in_memory, compare_netcdfs + +# compare in memory (faster) or full output files? or both! +do_compare_output_files = False +do_compare_in_memory = True +rtol = atol = 1.0e-5 calc_methods = ("numpy", "numba", "fortran") params = ("params_sep", "params_one") @@ -44,7 +49,6 @@ def test_compare_prms( ) tmp_path = pl.Path(tmp_path) - output_dir = domain["prms_output_dir"] # get the answer data comparison_var_names = [ @@ -56,19 +60,20 @@ def test_compare_prms( "hru_intcpstor", "hru_intcpevap", "intcp_changeover", + "intcp_form", + # "intcp_transp_on", # not a prms variable + "hru_intcpstor_change", + "hru_intcpstor_old", ] - ans = {} - for key in comparison_var_names: - nc_pth = output_dir / f"{key}.nc" - ans[key] = adapter_factory(nc_pth, variable_name=key, control=control) - # setup the canopy + output_dir = domain["prms_output_dir"] + input_variables = {} for key in PRMSCanopy.get_inputs(): nc_pth = output_dir / f"{key}.nc" input_variables[key] = nc_pth - cnp = PRMSCanopy( + canopy = PRMSCanopy( control=control, discretization=discretization, parameters=parameters, @@ -77,34 +82,37 @@ def test_compare_prms( calc_method=calc_method, ) - all_success = True + if do_compare_output_files: + nc_parent = tmp_path / domain["domain_name"] + canopy.initialize_netcdf(nc_parent) + + if do_compare_in_memory: + answers = {} + for var in comparison_var_names: + var_pth = output_dir / f"{var}.nc" + answers[var] = adapter_factory( + var_pth, variable_name=var, control=control + ) + for istep in range(control.n_times): control.advance() - cnp.advance() - cnp.calculate(1.0) - - # compare along the way - atol = 1.0e-5 - for key, val in ans.items(): - val.advance() - for key in ans.keys(): - a1 = ans[key].current - a2 = cnp[key] - success = np.isclose(a2, a1, atol=atol).all() - if not success: - all_success = False - diff = a1 - a2 - diffmin = diff.min() - diffmax = diff.max() - print(f"time step {istep}") - print(f"output variable {key}") - print(f"prms {a1.min()} {a1.max()}") - print(f"pywatershed {a2.min()} {a2.max()}") - print(f"diff {diffmin} {diffmax}") - - cnp.finalize() - - if not all_success: - raise Exception("pywatershed results do not match prms results") + canopy.advance() + canopy.calculate(1.0) + canopy.output() + if do_compare_in_memory: + compare_in_memory( + canopy, answers, atol=atol, rtol=rtol, skip_missing_ans=True + ) + + canopy.finalize() + + if do_compare_output_files: + compare_netcdfs( + comparison_var_names, + tmp_path / domain["domain_name"], + output_dir, + atol=atol, + rtol=rtol, + ) return diff --git a/autotest/test_prms_canopy_runoff.py b/autotest/test_prms_canopy_runoff.py deleted file mode 100644 index 493e3b2c..00000000 --- a/autotest/test_prms_canopy_runoff.py +++ /dev/null @@ -1,145 +0,0 @@ -import pathlib as pl - -import numpy as np -import pytest - -from pywatershed.base.adapter import adapter_factory -from pywatershed.base.control import Control -from pywatershed.hydrology.prms_canopy import PRMSCanopy -from pywatershed.hydrology.prms_runoff import PRMSRunoff -from pywatershed.parameters import PrmsParameters - - -@pytest.fixture(scope="function") -def params(domain): - return PrmsParameters.load(domain["param_file"]) - - -@pytest.fixture(scope="function") -def control(domain): - return Control.load(domain["control_file"]) - - -def test_canopy_runoff(domain, control, params, tmp_path): - tmp_path = pl.Path(tmp_path) - - # get the answer data - - comparison_var_names = [ - "infil", - "dprst_stor_hru", - "hru_impervstor", - "sroff", - ] - output_dir = domain["prms_output_dir"] - - # Read PRMS output into ans for comparison with pywatershed results - ans = {} - for key in comparison_var_names: - nc_pth = output_dir / f"{key}.nc" - ans[key] = adapter_factory(nc_pth, variable_name=key, control=control) - - # instantiate canopy - input_variables = {} - for key in PRMSCanopy.get_inputs(): - nc_pth = output_dir / f"{key}.nc" - input_variables[key] = nc_pth - - canopy = PRMSCanopy( - control=control, - discretization=None, - parameters=params, - **input_variables, - ) - - # instantiate runoff - input_variables = {} - for key in PRMSRunoff.get_inputs(): - nc_pth = output_dir / f"{key}.nc" - if "soil_moist" in str(nc_pth): - nc_pth = output_dir / "soil_moist_prev.nc" - input_variables[key] = nc_pth - input_variables["net_ppt"] = None - input_variables["net_rain"] = None - input_variables["net_snow"] = None - runoff = PRMSRunoff( - control=control, - discretization=None, - parameters=params, - **input_variables, - ) - - # wire up output from canopy as input to runoff - runoff.set_input_to_adapter( - "net_ppt", - adapter_factory(canopy.net_ppt, "net_ppt", control=control), - ) - runoff.set_input_to_adapter( - "net_rain", - adapter_factory(canopy.net_rain, "net_rain", control=control), - ) - runoff.set_input_to_adapter( - "net_snow", - adapter_factory(canopy.net_snow, "net_snow", control=control), - ) - - all_success = True - for istep in range(control.n_times): - control.advance() - canopy.advance() - runoff.advance() - - canopy.calculate(1.0) - runoff.calculate(1.0) - - # advance the answer, which is being read from a netcdf file - for key, val in ans.items(): - val.advance() - - # make a comparison check with answer - check = True - failfast = True - detailed = True - if check: - atol = 1.0e-5 - success = check_timestep_results( - runoff, istep, ans, atol, detailed - ) - if not success: - all_success = False - if failfast: - assert success, "stopping..." - - runoff.finalize() - - # check at the end and error if one or more steps didn't pass - if not all_success: - raise Exception("pywatershed results do not match prms results") - - return - - -def check_timestep_results(storageunit, istep, ans, atol, detailed=False): - all_success = True - for key in ans.keys(): - a1 = ans[key].current - a2 = storageunit[key] - success = np.isclose(a2, a1, atol=atol).all() - if not success: - all_success = False - diff = a1 - a2 - diffmin = diff.min() - diffmax = diff.max() - if True: - print(f"time step {istep}") - print(f"output variable {key}") - print(f"prms {a1.min()} {a1.max()}") - print(f"pywatershed {a2.min()} {a2.max()}") - print(f"diff {diffmin} {diffmax}") - if detailed: - idx = np.where(np.abs(diff) > atol)[0] - for i in idx: - print( - f"hru {i} prms {a1[i]} pywatershed {a2[i]} diff {diff[i]}" - ) - return all_success diff --git a/autotest/test_prms_channel.py b/autotest/test_prms_channel.py index e21bdae1..21af0266 100644 --- a/autotest/test_prms_channel.py +++ b/autotest/test_prms_channel.py @@ -9,13 +9,11 @@ from pywatershed.parameters import PrmsParameters from utils_compare import compare_in_memory, compare_netcdfs -# compare in memory (faster) or full output files? +# compare in memory (faster) or full output files? or both! do_compare_output_files = True do_compare_in_memory = True rtol = atol = 1.0e-7 -fail_fast = False - calc_methods = ("numpy", "numba", "fortran") params = ("params_sep", "params_one") @@ -57,9 +55,8 @@ def test_compare_prms( ) tmp_path = pl.Path(tmp_path) - - # load csv files into dataframes output_dir = domain["prms_output_dir"] + input_variables = {} for key in PRMSChannel.get_inputs(): nc_path = output_dir / f"{key}.nc" diff --git a/autotest/test_prms_runoff.py b/autotest/test_prms_runoff.py index f0fe0de6..8a7a7937 100644 --- a/autotest/test_prms_runoff.py +++ b/autotest/test_prms_runoff.py @@ -1,12 +1,17 @@ import pathlib as pl -import numpy as np import pytest from pywatershed.base.adapter import adapter_factory from pywatershed.base.control import Control from pywatershed.hydrology.prms_runoff import PRMSRunoff from pywatershed.parameters import Parameters, PrmsParameters +from utils_compare import compare_in_memory, compare_netcdfs + +# compare in memory (faster) or full output files? or both! +do_compare_output_files = False +do_compare_in_memory = True +rtol = atol = 1.0e-10 calc_methods = ("numpy", "numba") params = ("params_sep", "params_one") @@ -40,36 +45,16 @@ def test_compare_prms( ): tmp_path = pl.Path(tmp_path) - # get the answer data - - comparison_var_names = [ - "infil", - "infil_hru", - "dprst_stor_hru", - "dprst_seep_hru", - "hru_impervstor", - "sroff", - "dprst_evap_hru", - "hru_impervevap", - "dprst_insroff_hru", - "dprst_stor_hru", - "contrib_fraction", - "hru_sroffp", - "hru_sroffi", - # "hru_impervstor_change", - "dprst_sroff_hru", - "dprst_insroff_hru", - # "dprst_stor_hru_change", - ] - output_dir = domain["prms_output_dir"] + comparison_var_names = set(PRMSRunoff.get_variables()) + # TODO: get rid of this exception + comparison_var_names -= set( + [ + "dprst_area_open", + ] + ) - # Read PRMS output into ans for comparison with pywatershed results - ans = {} - for key in comparison_var_names: - nc_pth = output_dir / f"{key}.nc" - ans[key] = adapter_factory(nc_pth, variable_name=key, control=control) + output_dir = domain["prms_output_dir"] - # instantiate runoff input_variables = {} for key in PRMSRunoff.get_inputs(): nc_pth = output_dir / f"{key}.nc" @@ -79,66 +64,45 @@ def test_compare_prms( control=control, discretization=discretization, parameters=parameters, - calc_method=calc_method, **input_variables, budget_type="error", + calc_method=calc_method, ) - all_success = True + if do_compare_output_files: + nc_parent = tmp_path / domain["domain_name"] + runoff.initialize_netcdf(nc_parent) + # test that init netcdf twice raises a warning + with pytest.warns(UserWarning): + runoff.initialize_netcdf(nc_parent) + + if do_compare_in_memory: + answers = {} + for var in comparison_var_names: + var_pth = output_dir / f"{var}.nc" + answers[var] = adapter_factory( + var_pth, variable_name=var, control=control + ) + for istep in range(control.n_times): control.advance() runoff.advance() runoff.calculate(1.0) - - # advance the answer, which is being read from a netcdf file - for key, val in ans.items(): - val.advance() - - # make a comparison check with answer - check = True - failfast = True - detailed = True - if check: - atol = 1.0e-10 - success = check_timestep_results( - runoff, istep, ans, atol, detailed + runoff.output() + if do_compare_in_memory: + compare_in_memory( + runoff, answers, atol=atol, rtol=rtol, skip_missing_ans=True ) - if not success: - all_success = False - if failfast: - assert success, "stopping..." runoff.finalize() - # check at the end and error if one or more steps didn't pass - if not all_success: - raise Exception("pywatershed results do not match prms results") + if do_compare_output_files: + compare_netcdfs( + comparison_var_names, + tmp_path / domain["domain_name"], + output_dir, + atol=atol, + rtol=rtol, + ) return - - -def check_timestep_results(storageunit, istep, ans, atol, detailed=False): - all_success = True - for key in ans.keys(): - a1 = ans[key].current - a2 = storageunit[key] - success = np.isclose(a1, a2, atol=atol, rtol=atol).all() - if not success: - all_success = False - diff = a1 - a2 - diffmin = diff.min() - diffmax = diff.max() - if True: - print(f"time step {istep}") - print(f"output variable {key}") - print(f"prms {a1.min()} {a1.max()}") - print(f"pywatershed {a2.min()} {a2.max()}") - print(f"diff {diffmin} {diffmax}") - if detailed: - idx = np.where(np.abs(diff) > atol)[0] - for i in idx: - print( - f"hru {i} prms {a1[i]} pywatershed {a2[i]} diff {diff[i]}" - ) - asdf - return all_success diff --git a/autotest/test_prms_snow.py b/autotest/test_prms_snow.py index 73c9a3e7..fb76f25a 100644 --- a/autotest/test_prms_snow.py +++ b/autotest/test_prms_snow.py @@ -1,13 +1,17 @@ import pathlib as pl -import numpy as np import pytest from pywatershed.base.adapter import adapter_factory from pywatershed.base.control import Control -from pywatershed.constants import epsilon32, zero from pywatershed.hydrology.prms_snow import PRMSSnow from pywatershed.parameters import Parameters, PrmsParameters +from utils_compare import compare_in_memory, compare_netcdfs + +# compare in memory (faster) or full output files? or both! +do_compare_output_files = False +do_compare_in_memory = True +rtol = atol = 1.0e-3 calc_methods = ("numpy", "numba") params = ("params_sep", "params_one") @@ -41,49 +45,51 @@ def test_compare_prms( domain, control, discretization, parameters, tmp_path, calc_method ): tmp_path = pl.Path(tmp_path) - output_dir = domain["prms_output_dir"] # get the answer data comparison_var_names = [ - # # "ai", - # "albedo", - # # "frac_swe", - # # "freeh2o", - # # "iasw", - # # "int_alb", + # 'ai', + # 'albedo', + # 'frac_swe', + # 'freeh2o', + # 'freeh2o_change', + # 'freeh2o_prev', + # 'iasw', + # 'int_alb', "iso", - # # "lso", - # # "lst", - # # "mso", - # # "pk_def", - # # "pk_den", - # # "pk_depth", - # # "pk_ice", - # # "pk_precip", - # # "pk_temp", - # # "pksv", - # "pkwater_ante", -- dosent respect the iso memory eval + # 'lso', + # 'lst', + # 'mso', + # 'newsnow', + # 'pk_def', + # 'pk_den', + # 'pk_depth', + # 'pk_ice', + # 'pk_ice_change', + # 'pk_ice_prev', + # 'pk_precip', + # 'pk_temp', + # 'pksv', + # 'pkwater_ante', "pkwater_equiv", - # "pptmix_nopack", - # # "pss", - # "pst", - # # "salb", - # # "scrv", - # # "slst", + # 'pkwater_equiv_change', + # 'pptmix_nopack', + # 'pss', + # 'pst', + # 'salb', + # 'scrv', + # 'slst', "snow_evap", - # "snowcov_area", # issues on 12-23-1979 in drb - # "snowcov_areasv", - # "snowmelt", # issues on 12-22-1979 in drb - # #"snsv", + # 'snowcov_area', + # 'snowcov_areasv', + # 'snowmelt', + # 'snsv', "tcal", + "through_rain", ] - ans = {} - for key in comparison_var_names: - nc_pth = output_dir / f"{key}.nc" - ans[key] = adapter_factory(nc_pth, variable_name=key, control=control) + output_dir = domain["prms_output_dir"] - # setup the snow input_variables = {} for key in PRMSSnow.get_inputs(): nc_path = output_dir / f"{key}.nc" @@ -94,137 +100,41 @@ def test_compare_prms( discretization, parameters, **input_variables, - budget_type="warn", + budget_type="error", calc_method=calc_method, ) - all_success = True - iso_censor_mask = None - for istep in range(control.n_times): - # for istep in range(10): + if do_compare_output_files: + nc_parent = tmp_path / domain["domain_name"] + snow.initialize_netcdf(nc_parent) - print("\n") - print(f"solving time: {control.current_time}") + if do_compare_in_memory: + answers = {} + for var in comparison_var_names: + var_pth = output_dir / f"{var}.nc" + answers[var] = adapter_factory( + var_pth, variable_name=var, control=control + ) + for istep in range(control.n_times): control.advance() snow.advance() - - snow.calculate(float(istep)) - - # compare along the way - for key, val in ans.items(): - val.advance() - - # pkwater_equiv --- - # first check pkwater_equiv and assert they are with - # their own tolerance (approx float error) - key = "pkwater_equiv" - pkwe_ans = ans[key].current - pkwe = snow[key] - pkwater_absdiff = abs(pkwe - pkwe_ans) - atol = 5e-2 # 5 hundreds of an inch - - close = np.isclose(pkwe, pkwe_ans, atol=atol, rtol=zero) - if iso_censor_mask is None: - assert close.all() - else: - assert (close[~iso_censor_mask]).all() - - # if one of them is "zero" and the other is not, - # then there can be differences in other variables. - # so we'll exclude this case from the comparison of - # other variables as long as the SWEs pass the above test - pkwater_equiv_one_zero = (pkwe <= atol) | (pkwe_ans <= atol) - if iso_censor_mask is not None: - pkwater_equiv_one_zero = np.where( - iso_censor_mask, True, pkwater_equiv_one_zero + snow.calculate(1.0) + snow.output() + if do_compare_in_memory: + compare_in_memory( + snow, answers, atol=atol, rtol=rtol, skip_missing_ans=True ) - print( - f"pkwater_equiv_one_zero.sum(): " f"{pkwater_equiv_one_zero.sum()}" - ) - - # iso --- - # iso has a memory of packwater equiv which can cause about - # a 20% change in albedo based on pkwater differences in the - # noise... AND those differences in albedo can translate into - # longer term differences in other variables. The solution is - # to flag (memory) when iso gets out of snyc over reasonable/small - # pkwater equiv differences no not evaluate those points and - # after pk water has gone to zero for both pywatershed and prms. - # we'll want to report the number of points being censored this - # way at each time. this censor mask has to be used for pkwater - # above - key = "iso" - iso_ans = ans[key].current - iso = snow[key] - iso_absdiff = abs(iso - iso_ans) - - # anytime at least one of them is zero AND there - # are iso differences, then we add these points to the - # iso_censor list - iso_censor_mask_new = pkwater_equiv_one_zero & (iso_absdiff > 0) - - if iso_censor_mask is None: - iso_censor_mask = iso_censor_mask_new - else: - iso_censor_mask = iso_censor_mask | iso_censor_mask_new - - # reset points where near zero and have matching iso - pkwater_both_zero_and_iso = ( - (pkwe <= epsilon32) & (pkwe_ans <= epsilon32) & (iso == iso_ans) - ) - iso_censor_mask = np.where( - pkwater_both_zero_and_iso, False, iso_censor_mask - ) - print(f"iso_censor_mask.sum(): {iso_censor_mask.sum()}") - assert np.where(~iso_censor_mask, iso - iso_ans == 0, True).all() - - censor_comp = iso_censor_mask | pkwater_equiv_one_zero - - # The rest of the variables to be checked - atol = 1e-5 - for key in ans.keys(): - if key in ["pkwater_equiv", "iso"]: - continue - - print(key) - a1 = ans[key].current - a2 = snow[key] - - success = np.allclose(a1, a2, atol=atol, rtol=zero) - atol = 1e-3 - success_prelim = np.isclose(a1, a2, atol=atol) - success = np.where(~censor_comp, success_prelim, True).all() - - if not success: - all_success = False - diff = a1 - a2 - diffmin = diff.min() - diffmax = diff.max() - print(f"time step {istep}") - print(f"output variable: {key}") - print(f"prms {a1.min()} {a1.max()}") - print(f"pywatershed {a2.min()} {a2.max()}") - print(f"diff {diffmin} {diffmax}") - - zz = abs(a2 - a1) - max_diff = np.where(~censor_comp, zz, zero).max() - wh_max_diff = np.where(zz == max_diff) - print(f"wh_max_diff: {wh_max_diff}") - print(f"max diff: {zz[wh_max_diff]}") - print(f"prms: {a1[wh_max_diff]}") - print(f"pywatershed: {a2[wh_max_diff]}") - print( - f"pkwater_absdiff[wh_max_diff]: " - f"{pkwater_absdiff[wh_max_diff]}" - ) - print(f"pkwe[wh_max_diff]: {pkwe[wh_max_diff]}") - print(f"pkwe_ans[wh_max_diff]: {pkwe_ans[wh_max_diff]}") - assert success snow.finalize() - if not all_success: - raise Exception("pywatershed results do not match prms results") + if do_compare_output_files: + compare_netcdfs( + comparison_var_names, + tmp_path / domain["domain_name"], + output_dir, + atol=atol, + rtol=rtol, + ) return diff --git a/autotest/test_prms_soilzone.py b/autotest/test_prms_soilzone.py index a296bbd2..a3c7acef 100644 --- a/autotest/test_prms_soilzone.py +++ b/autotest/test_prms_soilzone.py @@ -1,12 +1,17 @@ import pathlib as pl -import numpy as np import pytest from pywatershed.base.adapter import adapter_factory from pywatershed.base.control import Control from pywatershed.hydrology.prms_soilzone import PRMSSoilzone from pywatershed.parameters import Parameters, PrmsParameters +from utils_compare import compare_in_memory, compare_netcdfs + +# compare in memory (faster) or full output files? or both! +do_compare_output_files = False +do_compare_in_memory = True +rtol = atol = 5.0e-6 calc_methods = ("numpy", "numba") params = ("params_sep", "params_one") @@ -39,39 +44,26 @@ def test_compare_prms( domain, control, discretization, parameters, tmp_path, calc_method ): tmp_path = pl.Path(tmp_path) + + comparison_var_names = list( + set(PRMSSoilzone.get_variables()) + # these are not prms variables per se + # the hru ones have non-hru equivalents being checked + # soil_zone_max and soil_lower_max would be nice to check but + # prms5.2.1 wont write them as hru variables + - set( + [ + "perv_actet_hru", + "soil_lower_change_hru", + "soil_lower_max", + "soil_rechr_change_hru", + "soil_zone_max", # not a prms variable? + ] + ) + ) + output_dir = domain["prms_output_dir"] - # get the answer data - comparison_var_names = [ - "cap_infil_tot", - "hru_actet", - "perv_actet", - "potet_lower", - "potet_rechr", - "recharge", - "slow_flow", - "slow_stor", - "soil_lower", - # "soil_lower_ratio", - "soil_moist", - "soil_moist_prev", - "soil_moist_tot", - "soil_rechr", - "soil_to_gw", - "soil_to_ssr", - "ssr_to_gw", - "ssres_flow", - "ssres_in", - "ssres_stor", - ### "unused_potet", - ] - - ans = {} - for key in comparison_var_names: - nc_pth = output_dir / f"{key}.nc" - ans[key] = adapter_factory(nc_pth, variable_name=key, control=control) - - # setup the soilzone input_variables = {} for key in PRMSSoilzone.get_inputs(): nc_path = output_dir / f"{key}.nc" @@ -85,58 +77,38 @@ def test_compare_prms( budget_type="error", calc_method=calc_method, ) - all_success = True - for istep in range(control.n_times): - # for istep in range(10): - control.advance() + if do_compare_output_files: + nc_parent = tmp_path / domain["domain_name"] + soil.initialize_netcdf(nc_parent) - # print("\n") - # print(control.current_time) + if do_compare_in_memory: + answers = {} + for var in comparison_var_names: + var_pth = output_dir / f"{var}.nc" + answers[var] = adapter_factory( + var_pth, variable_name=var, control=control + ) + for istep in range(control.n_times): + control.advance() soil.advance() - soil.calculate(float(istep)) - - # compare along the way - atol = 1.0e-4 - for key, val in ans.items(): - val.advance() - for key in ans.keys(): - a1 = ans[key].current.data - a2 = soil[key] - success = np.isclose(a1, a2, atol=atol).all() - # if success: - # print(f"SUCCESS output variable: {key}") - if not success: - all_success = False - # asdf - diff = a1 - a2 - diffmin = diff.min() - diffmax = diff.max() - # print(f"time step {istep}") - print(f"fail output variable: {key}") - print(f"prms {a1.min()} {a1.max()}") - print(f"pywatershed {a2.min()} {a2.max()}") - print(f"diff {diffmin} {diffmax}") - - # if istep == 15: - # asdf + soil.calculate(1.0) + soil.output() + if do_compare_in_memory: + compare_in_memory( + soil, answers, atol=atol, rtol=rtol, skip_missing_ans=True + ) soil.finalize() - if not all_success: - raise Exception("pywatershed results do not match prms results") - - # gw.output() - - # gw.finalize() - - # assert_error = False - # for key, (base, compare) in output_compare.items(): - # success, diff = NetCdfCompare(base, compare).compare() - # if not success: - # print(f"comparison for {key} failed: maximum error {diff}") - # assert_error = True - # assert not assert_error, "comparison failed" + if do_compare_output_files: + compare_netcdfs( + comparison_var_names, + tmp_path / domain["domain_name"], + output_dir, + atol=atol, + rtol=rtol, + ) return diff --git a/autotest/test_prms_solar_geom.py b/autotest/test_prms_solar_geom.py index cb2bd971..076c272b 100644 --- a/autotest/test_prms_solar_geom.py +++ b/autotest/test_prms_solar_geom.py @@ -1,4 +1,3 @@ -import numpy as np import pytest from pywatershed.atmosphere.prms_solar_geometry import PRMSSolarGeometry @@ -8,8 +7,11 @@ from pywatershed.parameters import PrmsParameters from utils_compare import compare_in_memory, compare_netcdfs -# in this case we'll compare netcdf files and in memory -atol = rtol = np.finfo(np.float32).resolution +# compare in memory (faster) or full output files? or both! +do_compare_output_files = False +do_compare_in_memory = True + +rtol = atol = 1.0e-10 params = ("params_sep", "params_one") @@ -43,6 +45,7 @@ def test_compare_prms( domain, control, discretization, parameters, tmp_path, from_prms_file ): output_dir = domain["prms_output_dir"] + prms_soltab_file = domain["prms_run_dir"] / "soltab_debug" if from_prms_file: from_prms_file = prms_soltab_file @@ -56,33 +59,37 @@ def test_compare_prms( from_prms_file=from_prms_file, netcdf_output_dir=tmp_path, ) - solar_geom.output() - solar_geom.finalize() - - compare_netcdfs( - PRMSSolarGeometry.get_variables(), - tmp_path, - output_dir, - atol=atol, - rtol=rtol, - ) - answers = {} - for var in PRMSSolarGeometry.get_variables(): - var_pth = output_dir / f"{var}.nc" - answers[var] = adapter_factory( - var_pth, variable_name=var, control=control + if do_compare_in_memory: + answers = {} + for var in PRMSSolarGeometry.get_variables(): + var_pth = output_dir / f"{var}.nc" + answers[var] = adapter_factory( + var_pth, variable_name=var, control=control + ) + + sunhrs_id = id(solar_geom.soltab_sunhrs) + + # Though the data is all calculate on the initial advance, + # we step through it using the timeseries array. + # we only need to output at time 0 + for ii in range(control.n_times): + control.advance() + solar_geom.advance() + if ii == 0: + solar_geom.output() + solar_geom.calculate(1.0) + + compare_in_memory(solar_geom, answers, atol=atol, rtol=rtol) + assert id(solar_geom.soltab_sunhrs) == sunhrs_id + + if do_compare_output_files: + compare_netcdfs( + PRMSSolarGeometry.get_variables(), + tmp_path, + output_dir, + atol=atol, + rtol=rtol, ) - # check the advance/calculate the state - sunhrs_id = id(solar_geom.soltab_sunhrs) - - for ii in range(control.n_times): - control.advance() - solar_geom.advance() - solar_geom.calculate(1.0) - - compare_in_memory(solar_geom, answers, atol=atol, rtol=rtol) - assert id(solar_geom.soltab_sunhrs) == sunhrs_id - return diff --git a/autotest/utils_compare.py b/autotest/utils_compare.py index bfe5e747..1a31985d 100644 --- a/autotest/utils_compare.py +++ b/autotest/utils_compare.py @@ -15,6 +15,8 @@ def assert_allclose( strict: bool = False, also_check_w_np: bool = True, error_message: str = "Comparison unsuccessful (default message)", + print_max_errs: bool = False, + var_name: str = "", ): """Reinvent np.testing.assert_allclose to get useful diagnostincs in debug @@ -29,6 +31,7 @@ def assert_allclose( handling of scalars mentioned in the Notes section is disabled. also_check_w_np: first check using np.testing.assert_allclose using the same options. + print_max_errs: bool=False. Print max abs and rel err for each var. """ if also_check_w_np: @@ -61,6 +64,11 @@ def assert_allclose( close = abs_close | rel_close + if print_max_errs: + sp = "" if len(var_name) == 0 else " " + print(f"{var_name}{sp}max abs err: {abs_diff.max()}") + print(f"{var_name}{sp}max rel err: {rel_abs_diff.max()}") + assert close.all() @@ -73,10 +81,18 @@ def compare_in_memory( strict: bool = False, also_check_w_np: bool = True, error_message: str = None, + skip_missing_ans: bool = False, ): # TODO: docstring for var in process.get_variables(): + if var not in answers.keys(): + if skip_missing_ans: + continue + else: + msg = f"Variable '{var}' not found in the answers provided." + raise KeyError(msg) + answers[var].advance() if isinstance(process[var], pws.base.timeseries.TimeseriesArray): @@ -106,6 +122,7 @@ def compare_netcdfs( strict: bool = False, also_check_w_np: bool = True, error_message: str = None, + print_var_max_errs: bool = False, ): # TODO: docstring # TODO: improve error message @@ -130,4 +147,6 @@ def compare_netcdfs( strict=strict, also_check_w_np=also_check_w_np, error_message=error_message, + print_max_errs=print_var_max_errs, + var_name=var, ) diff --git a/doc/conf.py b/doc/conf.py index 4e4a047d..efce9f6e 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -97,14 +97,13 @@ # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. # -# html_theme = "pydata_sphinx_theme" html_theme = "sphinx_book_theme" html_title = "pywatershed" html_context = { - "github_user": "pydata", - "github_repo": "xarray", + "github_user": "EC-USGS", + "github_repo": "pywatershed", "github_version": "main", "doc_path": "doc", } diff --git a/doc/index.rst b/doc/index.rst index 1d7f200e..0ef4a6d8 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -4,12 +4,14 @@ pywatershed: a hydrologic model in Python #################################################### -Welcome to `pywatershed` docs! +Welcome to the `pywatershed` docs! -Browse the API reference, developer info, and index using the table of -contents on the left. +To learn more about this project, see `About `_. -| For introductory example notebooks, look in the `examples/ `_ directory in the repository. Numbered starting at 00, these are meant to be completed in order. Though no notebook outputs are saved in Github, these notebooks can easily navigated to and run in WholeTale containers (free but sign-up or log-in required): +Please browse the API reference, developer info, and +index using the table of contents on the left. + +| For introductory example notebooks, look in the `examples/ `_ directory in the repository. Numbered starting at 00, these are meant to be completed in order. Though no notebook outputs are saved in Github, these notebooks can be easily found and run in WholeTale containers (free but sign-up or log-in required): .. image:: https://raw.githubusercontent.com/whole-tale/wt-design-docs/master/badges/wholetale-explore.svg :target: https://dashboard.wholetale.org @@ -24,6 +26,12 @@ We value your feedback! To view the repository, suggest an edit to this document Thank you for your interest. +.. toctree:: + :hidden: + :caption: About + + About + .. toctree:: :hidden: :caption: API Reference diff --git a/doc/whats-new.rst b/doc/whats-new.rst index 16a5ca9c..3f62481c 100644 --- a/doc/whats-new.rst +++ b/doc/whats-new.rst @@ -23,7 +23,7 @@ New features (:pull:`232`) By `James McCreight `_. - Conda feedstock for pywatershed ``_. By `Joseph Hughes `_. - + Breaking changes ~~~~~~~~~~~~~~~~ @@ -41,14 +41,22 @@ Bug fixes ~~~~~~~~~ - Resolve issues with different ways of specifiying netcdf output options. (:pull:`230`) By `James McCreight `_. - +- PRMSSoilzone remove soil_moist_prev because soil_moist is not prognotic and + PRMSRunoff was needing it in the advance and not getting the correct value. + PRMSRunoff now depends on soil_lower_prev and soil_rechr_prev instead. + (:pull:`244`) By `James McCreight `_. Documentation ~~~~~~~~~~~~~ +- Add about section for version 1.0 to describe how pywatershed matches PRMS' + NHM configuration and how to perform the comparison. + (:pull:`244`) By `James McCreight `_. Internal changes ~~~~~~~~~~~~~~~~ +- Refactor tests against PRMS for consistency, flexibility, and thoroughness. + (:pull:`244`) By `James McCreight `_. .. _whats-new.0.2.1: diff --git a/examples/03_compare_pws_prms.ipynb b/examples/03_compare_pws_prms.ipynb new file mode 100644 index 00000000..b0bb3a08 --- /dev/null +++ b/examples/03_compare_pws_prms.ipynb @@ -0,0 +1,607 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "2d33f0d3-b07c-4fcb-8fd8-1d23014ec3cb", + "metadata": {}, + "source": [ + "# Compare pywatershed and PRMS\n", + "\n", + "This notebook compares pywatershed (PWS) and PRMS outputs. As part of release of v1.0 this notebook is meant to give users insight on how wimilar results from pywatershed and PRMS are. They are identical or very close in the vast majority of cases. This notebook makes it easy to find when they are not by providing statistics at individual HRUs and timeseries for all HRUs.\n", + "\n", + "Note that this notebook requires an editable install of pywatershed (`pip install -e` in the pywatershed repository root) for the requisite data. PRMS/NHM domains which may be used are in the `test_data/` directory of pywatershed (`hru_1`, `drb_2yr`, and `ucb_2yr`) but any other domain may be used. \n", + "\n", + "## Notes on setting up other domains\n", + "You may want to supply your own domain and see how pywatershed works on it. Here are notes on doing so. Domains must supply the correct, required files in `test_data/your_domain` which are given in this listing:\n", + "\n", + "```\n", + "control.test prcp.cbh sf_data tmax.nc tmin.nc\n", + "myparam.param prcp.nc tmax.cbh tmin.cbh\n", + "```\n", + "\n", + "The `*.cbh` files must be pre-converted to netcdf for `prcp`, `tmin`, and `tmax` and how to do this can be found near the top of notebook 02. The `control.test` and `myparam.param` files are used by both PRMS and PWS. The `control.test` files in the repo are specific for being able to run sub-models and include a nearly maximal amount of model output (time-inefficient for both PRMS and PWS). The stock control files can be found in `test_data/common` there is a file for single-hru domains and multi-hru domains and these are identical (as appropriate) for the domains included in the repository. For running a large domain, for example, it is desirable to reduce the total amount of output (but this may not allow for PWS sub-models to be run as PRMS dosent necessarily supply all the required fields). So you may modify the `control.test` file but take careful note of what options are available in pywatershed as currently only NHM configuration is available.\n", + "\n", + "The runs of PRMS use double precision binaries produced by the `prms_src/prms5.2.1` source code in the pywatershed repository. The procedure used below is exactly as done in CI for running regression tests against PRMS.\n", + "\n", + "All of the code required for plotting below is included so that it can be further tailored to your tastes.\n", + "\n", + "## Imports, etc" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7275fe95-8e0b-45f9-837c-ad6f7d38ced7", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + } + }, + "outputs": [], + "source": [ + "# auto-format the code in this notebook\n", + "%load_ext jupyter_black" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "816e410d-c806-4467-8d56-e6add8c7f516", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + } + }, + "outputs": [], + "source": [ + "import pathlib as pl\n", + "from platform import processor\n", + "from pprint import pprint\n", + "from shutil import rmtree\n", + "import subprocess\n", + "from sys import platform\n", + "import warnings\n", + "\n", + "import hvplot.pandas # noqa\n", + "import hvplot.xarray # noqa\n", + "import numpy as np\n", + "import pandas as pd\n", + "import pywatershed as pws\n", + "import xarray as xr\n", + "\n", + "repo_root = pws.constants.__pywatershed_root__.parent\n", + "nb_output_dir = pl.Path(\"./03_compare_pws_prms\")" + ] + }, + { + "cell_type": "markdown", + "id": "06ff55cb-6298-4bbc-b654-46568db7f1fa", + "metadata": {}, + "source": [ + "## Configuration\n", + "\n", + "Specify what you want!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4950eb10-25c0-403c-80ef-857cd003bc5c", + "metadata": {}, + "outputs": [], + "source": [ + "domain_name: str = \"drb_2yr\" # must be present in test_data/domain_name\n", + "calc_method: str = \"numba\"\n", + "budget_type: str = None\n", + "\n", + "run_prms: bool = True ## always forced/overwrite\n", + "\n", + "run_pws: bool = True # run if the output does not exist on disk\n", + "force_pws_run: bool = True # if it exists on disk, re-run it and overwrite?" + ] + }, + { + "cell_type": "markdown", + "id": "7fb91db9-3933-420a-a25c-166f3aa85981", + "metadata": {}, + "source": [ + "## Run PRMS" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "007b24ce-c8fe-4409-9880-13746acf113b", + "metadata": {}, + "outputs": [], + "source": [ + "domain_dir = repo_root / f\"test_data/{domain_name}\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7d873087-b3aa-42c6-b16f-7ee8d213e849", + "metadata": {}, + "outputs": [], + "source": [ + "# use pytest to run the domains as in CI\n", + "if run_prms:\n", + " print(f\"PRMS running domain in {repo_root / f'test_data' / domain_name}\")\n", + " subprocess.run(\n", + " f\"pytest -s -n=2 run_prms_domains.py --domain={domain_name} -vv --force\",\n", + " shell=True,\n", + " cwd=repo_root / \"test_data/generate\",\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "256e81e0-ab88-408f-bf59-85531d5377a2", + "metadata": {}, + "outputs": [], + "source": [ + "# Convert PRMS output to netcdf as in CI\n", + "if run_prms:\n", + " if \"conus\" in domain_name:\n", + " nproc = 2 # memory bound for CONUS\n", + " conv_only = \"::make_netcdf_files\"\n", + " else:\n", + " nproc = 8 # processor bound otherwise\n", + " conv_only = \"\"\n", + "\n", + " subprocess.run(\n", + " f\"pytest -n={nproc} convert_prms_output_to_nc.py{conv_only} --domain={domain_name} --force\",\n", + " shell=True,\n", + " cwd=repo_root / \"test_data/generate\",\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "120e0be5-f30b-4152-9a40-2f7bfa775e79", + "metadata": {}, + "source": [ + "## Run pywatershed" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "529247fe-7f0b-4783-90a4-a4ddf55c1489", + "metadata": {}, + "outputs": [], + "source": [ + "if run_pws:\n", + " nhm_processes = [\n", + " pws.PRMSSolarGeometry, # submodles are possible\n", + " pws.PRMSAtmosphere,\n", + " pws.PRMSCanopy,\n", + " pws.PRMSSnow,\n", + " pws.PRMSRunoff,\n", + " pws.PRMSSoilzone,\n", + " pws.PRMSGroundwater,\n", + " pws.PRMSChannel,\n", + " ]\n", + "\n", + " if len(nhm_processes) == 8:\n", + " input_dir = domain_dir\n", + " run_dir = nb_output_dir / f\"{domain_name}_full_nhm\"\n", + " else:\n", + " input_dir = domain_dir / \"output\"\n", + " run_dir = nb_output_dir / f\"{domain_name}_subset_nhm\"\n", + "\n", + " control = pws.Control.load(domain_dir / \"control.test\")\n", + " params = pws.parameters.PrmsParameters.load(domain_dir / \"myparam.param\")\n", + "\n", + " if run_dir.exists():\n", + " if force_pws_run:\n", + " rmtree(run_dir)\n", + " else:\n", + " raise RuntimeError(\"run directory exists\")\n", + "\n", + " print(f\"PWS writing output to {run_dir}\")\n", + "\n", + " control.options = control.options | {\n", + " \"input_dir\": input_dir,\n", + " \"budget_type\": budget_type,\n", + " \"calc_method\": calc_method,\n", + " \"netcdf_output_dir\": run_dir,\n", + " }\n", + "\n", + " nhm = pws.Model(\n", + " nhm_processes,\n", + " control=control,\n", + " parameters=params,\n", + " )\n", + " nhm.run(finalize=True)" + ] + }, + { + "cell_type": "markdown", + "id": "aa221ffe-eb0a-4332-94fd-9826710adc74", + "metadata": {}, + "source": [ + "## Compare outputs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "97cc606f-4d69-439f-84cb-00b78ae13605", + "metadata": {}, + "outputs": [], + "source": [ + "def compare_var_timeseries(var_name, rmse_min=None):\n", + " \"\"\"Plots compare timeseries a PWS and PRMS variable for all locations in domain (scrollable).\n", + "\n", + " Args:\n", + " var_name: string name of variable\n", + " rmse_min: only plot locations which exceed this minimum rmse between PRMS and PWS.\n", + "\n", + " \"\"\"\n", + " from textwrap import fill\n", + "\n", + " var_meta = pws.meta.find_variables(var_name)[var_name]\n", + " ylabel = f\"{fill(var_meta['desc'], 40)}\\n({var_meta['units']})\"\n", + "\n", + " prms_file = domain_dir / f\"output/{var_name}.nc\"\n", + " if not prms_file.exists():\n", + " return None\n", + " prms_var = xr.open_dataarray(prms_file)\n", + " pws_var = xr.open_dataarray(run_dir / f\"{var_name}.nc\")\n", + "\n", + " if rmse_min is not None:\n", + " if \"time\" in prms_var.dims:\n", + " time_dim = \"time\"\n", + " else:\n", + " time_dim = \"doy\"\n", + "\n", + " rmse = np.sqrt((pws_var - prms_var).mean(dim=time_dim) ** 2)\n", + " mask_ge_min = rmse >= rmse_min\n", + " n_mask = len(np.where(mask_ge_min)[0])\n", + " print(f\"There are {n_mask} locations with RMSE > {rmse_min}\")\n", + " if n_mask == 0:\n", + " return None\n", + " prms_var = prms_var.where(mask_ge_min, drop=True)\n", + " pws_var = pws_var.where(mask_ge_min, drop=True)\n", + "\n", + " comp_ds = xr.merge(\n", + " [\n", + " prms_var.rename(\"prms\"),\n", + " pws_var.rename(\"pws\"),\n", + " ]\n", + " )\n", + " var_meta = pws.meta.find_variables(var_name)[var_name]\n", + " space_coord = list(comp_ds.coords)\n", + " for t_coord in [\"doy\", \"time\"]:\n", + " if t_coord in space_coord:\n", + " space_coord.remove(t_coord)\n", + "\n", + " display(\n", + " comp_ds.hvplot(\n", + " frame_width=800,\n", + " frame_height=500,\n", + " groupby=space_coord,\n", + " # title=title,\n", + " ylabel=ylabel,\n", + " group_label=\"Model\",\n", + " )\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0958ac67-a965-4990-8aba-41731846e3bf", + "metadata": {}, + "outputs": [], + "source": [ + "def calc_stat_location(var_name, stat_name):\n", + " \"\"\"Calculate a statistic location-wise (over time).\n", + "\n", + " Args:\n", + " var_name: str for the variable of interest\n", + " stat_name: one of [\"rmse\", \"rrmse\"]\n", + " \"\"\"\n", + " prms_file = domain_dir / f\"output/{var_name}.nc\"\n", + " if not prms_file.exists():\n", + " print(f\"PRMS file '{prms_file}' DNE, skipping.\")\n", + " return None\n", + " prms = xr.open_dataarray(prms_file, decode_timedelta=False)\n", + " pws_file = run_dir / f\"{var_name}.nc\"\n", + " assert pws_file.exists()\n", + " nhm_after = xr.open_dataarray(pws_file, decode_timedelta=False)\n", + " if \"time\" in prms.dims:\n", + " time_dim = \"time\"\n", + " else:\n", + " time_dim = \"doy\"\n", + " if stat_name.lower() == \"rmse\":\n", + " stat = np.sqrt((nhm_after - prms).mean(dim=time_dim) ** 2)\n", + " elif stat_name.lower() == \"rrmse\":\n", + " stat = np.sqrt(((nhm_after - prms) / prms).mean(dim=time_dim) ** 2)\n", + " return stat.to_dataframe().melt(ignore_index=False)\n", + "\n", + "\n", + "def box_jitter_plot(\n", + " df, subplot_width: int = 400, stat_name: str = \"Statistic\"\n", + "):\n", + " \"\"\"Box/violin-plot of a dataframe.\n", + "\n", + " Args:\n", + " df: a pd.Dataframe\n", + " subplot_width: int for how wide the subplots should be\n", + " stat_name: str of the statisitc name\n", + " \"\"\"\n", + " from textwrap import fill\n", + "\n", + " var_name = df.variable.iloc[0]\n", + " var_meta = pws.meta.find_variables(var_name)[var_name]\n", + " ylabel = (\n", + " f\"{stat_name} of\\n{fill(var_meta['desc'], 40)}\\n({var_meta['units']})\"\n", + " )\n", + " coord = df.index.name\n", + "\n", + " box = df.hvplot.violin(y=\"value\", by=\"variable\", legend=False)\n", + " jitter = df.hvplot.scatter(\n", + " y=\"value\",\n", + " x=\"variable\",\n", + " hover_cols=[coord],\n", + " )\n", + " return (box * jitter).opts(\n", + " width=subplot_width,\n", + " # xlabel=f\"over {coord}s\",\n", + " xlabel=\"\",\n", + " ylabel=ylabel,\n", + " )\n", + "\n", + "\n", + "def plot_proc_stats(\n", + " proc, stat_name: str = \"RMSE\", ncols: int = 5, subplot_width: int = 300\n", + "):\n", + " \"\"\"Plot pywatershed process stats.\n", + "\n", + " For a process (e.g. pws.PRMSRunoff), make box/violin plots of its stats for each of its (available) variables\n", + "\n", + " Args:\n", + " proc: a pws.Process subclass\n", + " stat_name: string of the statistic desired to be passed to box_jitter_plot\n", + " ncols: int number of columns in the plot\n", + " subplot_widt: int width of the subplots\n", + "\n", + " \"\"\"\n", + " var_stats = []\n", + " for var_name in proc.get_variables():\n", + " var_stats += [calc_stat_location(var_name, stat_name)]\n", + "\n", + " var_plots = [\n", + " box_jitter_plot(vv, subplot_width=subplot_width, stat_name=stat_name)\n", + " for vv in var_stats\n", + " if vv is not None\n", + " ]\n", + " if len(var_plots) == 0:\n", + " return None\n", + "\n", + " plot = var_plots[0]\n", + " for vv in var_plots[1:]:\n", + " plot += vv\n", + "\n", + " plot = plot.opts(shared_axes=False)\n", + " if len(var_plots) > 1:\n", + " plot = plot.cols(ncols)\n", + "\n", + " with warnings.catch_warnings():\n", + " warnings.simplefilter(\"ignore\")\n", + " display(plot)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "acadd2d2-3652-4443-8860-fa2cea2103db", + "metadata": {}, + "outputs": [], + "source": [ + "if pws.PRMSSolarGeometry in nhm_processes:\n", + " plot_proc_stats(pws.PRMSSolarGeometry, \"RMSE\", 5)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "858a8820-f086-4c9d-827b-159137aeedb6", + "metadata": {}, + "outputs": [], + "source": [ + "if pws.PRMSSolarGeometry in nhm_processes:\n", + " compare_var_timeseries(\"soltab_potsw\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "96057162-c53c-4124-9e6f-4ea41ead03a3", + "metadata": {}, + "outputs": [], + "source": [ + "if pws.PRMSAtmosphere in nhm_processes:\n", + " plot_proc_stats(pws.PRMSAtmosphere, \"RMSE\", 4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "88f6d07e-0c2c-4038-88e2-43bf3525fc94", + "metadata": {}, + "outputs": [], + "source": [ + "if pws.PRMSAtmosphere in nhm_processes:\n", + " compare_var_timeseries(\"tmaxf\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0e78b252-0e0a-49a9-bc15-d3acde820606", + "metadata": {}, + "outputs": [], + "source": [ + "if pws.PRMSCanopy in nhm_processes:\n", + " plot_proc_stats(pws.PRMSCanopy, \"RMSE\", 4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3066a8a2-d122-4100-bb33-2f672afa3dc1", + "metadata": {}, + "outputs": [], + "source": [ + "if pws.PRMSCanopy in nhm_processes:\n", + " compare_var_timeseries(\"intcp_stor\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d7d47572-1688-48f4-a916-17f63038a2fc", + "metadata": {}, + "outputs": [], + "source": [ + "if pws.PRMSSnow in nhm_processes:\n", + " plot_proc_stats(pws.PRMSSnow, \"RMSE\", 4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11a12175-0e7f-480d-a647-d58aff2494fc", + "metadata": {}, + "outputs": [], + "source": [ + "if pws.PRMSSnow in nhm_processes:\n", + " compare_var_timeseries(\"pkwater_equiv\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d3e1f0c3-2ed9-4a05-813a-4586fc8e3bd1", + "metadata": {}, + "outputs": [], + "source": [ + "if pws.PRMSRunoff in nhm_processes:\n", + " plot_proc_stats(pws.PRMSRunoff, \"RMSE\", 4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1f49fce8-112b-4dcc-8f67-26eaec71ca35", + "metadata": {}, + "outputs": [], + "source": [ + "if pws.PRMSRunoff in nhm_processes:\n", + " compare_var_timeseries(\"contrib_fraction\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "50de17cc-bdee-4d34-bf16-7289d0425bfa", + "metadata": {}, + "outputs": [], + "source": [ + "if pws.PRMSSoilzone in nhm_processes:\n", + " plot_proc_stats(pws.PRMSSoilzone, \"RMSE\", 4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d3248cbf-e4aa-44a2-9259-68b3b1faea67", + "metadata": {}, + "outputs": [], + "source": [ + "if pws.PRMSSoilzone in nhm_processes:\n", + " compare_var_timeseries(\"soil_rechr\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2c72ef1f-984a-48bd-b417-fb9754c09a50", + "metadata": {}, + "outputs": [], + "source": [ + "if pws.PRMSGroundwater in nhm_processes:\n", + " plot_proc_stats(pws.PRMSGroundwater, \"RMSE\", 4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0569deac-c07f-453d-9f61-588edc7ad355", + "metadata": {}, + "outputs": [], + "source": [ + "if pws.PRMSGroundwater in nhm_processes:\n", + " compare_var_timeseries(\"gwres_flow_vol\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "59b4b8e5-f944-496f-9f85-800677159272", + "metadata": {}, + "outputs": [], + "source": [ + "if pws.PRMSChannel in nhm_processes:\n", + " plot_proc_stats(pws.PRMSChannel, \"RMSE\", 4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8a230917-77ea-487f-9919-86c05ebcb964", + "metadata": {}, + "outputs": [], + "source": [ + "if pws.PRMSChannel in nhm_processes:\n", + " compare_var_timeseries(\"seg_outflow\") # , rmse_min=0.01)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c9b6a847-e7e8-412c-adea-fd63a516279b", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1d152230-b6cc-434a-9331-90799083dab7", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "language_info": { + "codemirror_mode": { + "name": "ipython" + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pywatershed/__init__.py b/pywatershed/__init__.py index 538b6170..d47f0389 100644 --- a/pywatershed/__init__.py +++ b/pywatershed/__init__.py @@ -18,13 +18,7 @@ from .hydrology.prms_snow import PRMSSnow from .hydrology.prms_soilzone import PRMSSoilzone from .hydrology.starfit import Starfit -from .utils import ( - ControlVariables, - NetCdfCompare, - NetCdfRead, - NetCdfWrite, - Soltab, -) +from .utils import ControlVariables, NetCdfRead, NetCdfWrite, Soltab from .utils.csv_utils import CsvFile from .version import __version__ diff --git a/pywatershed/hydrology/prms_canopy.f90 b/pywatershed/hydrology/prms_canopy.f90 index bbce57c7..82cdc7bf 100644 --- a/pywatershed/hydrology/prms_canopy.f90 +++ b/pywatershed/hydrology/prms_canopy.f90 @@ -18,7 +18,6 @@ pure subroutine calc_canopy( & hru_snow, & ! in intcp_changeover, & ! in intcp_evap, & ! in - intcp_form_in, & ! in - NOT USED AS INPUT? intcp_stor, & ! in intcp_transp_on, & ! in net_ppt, & ! in @@ -44,6 +43,7 @@ pure subroutine calc_canopy( & OFF, & ! in ACTIVE, & ! in intcp_evap_out, & ! out + intcp_form_out, & ! out intcp_stor_out, & ! out net_rain_out, & ! out net_snow_out, & ! out @@ -64,7 +64,7 @@ pure subroutine calc_canopy( & ! Input Vectors integer(kind=4), intent(in), dimension(nhru) :: & - intcp_form_in, intcp_transp_on + intcp_transp_on integer(kind=8), intent(in), dimension(nhru) :: & cov_type, hru_type @@ -76,7 +76,7 @@ pure subroutine calc_canopy( & ! Output vectors integer(kind=4), intent(out), dimension(nhru) :: & - intcp_transp_on_out + intcp_transp_on_out, intcp_form_out real(kind=8), intent(out), dimension(nhru) :: & intcp_evap_out, intcp_stor_out, net_rain_out, net_snow_out, net_ppt_out, & hru_intcpstor_out, hru_intcpevap_out, intcp_changeover_out @@ -88,8 +88,7 @@ pure subroutine calc_canopy( & netrain, netsnow, cov, stor_max_rain, & intcpstor, intcpevap, changeover, extra_water, diff, & epan_coef, evrn, evsn, zz, dd, & - last, intcpstor_in, netrain_in, netsnow_in, & - intcp_form + last, intcpstor_in, netrain_in, netsnow_in ! initialize output values ? ! intcp_evap_out = intcpevap @@ -114,11 +113,9 @@ pure subroutine calc_canopy( & stor_max_rain = wrain_intcp(ii) end if - intcp_form = RAIN - ! intcp_form(ii) = RAIN + intcp_form_out(ii) = RAIN if (hru_snow(ii) > zero) then - intcp_form = SNOW - ! intcp_form(ii) = SNOW + intcp_form_out(ii) = SNOW end if intcpstor = intcp_stor(ii) @@ -233,7 +230,7 @@ pure subroutine calc_canopy( & ! IF (( evrn 0) then intcpstor = zz diff --git a/pywatershed/hydrology/prms_canopy.py b/pywatershed/hydrology/prms_canopy.py index da6acb51..f8cdbc81 100644 --- a/pywatershed/hydrology/prms_canopy.py +++ b/pywatershed/hydrology/prms_canopy.py @@ -7,7 +7,7 @@ from ..base.adapter import adaptable from ..base.conservative_process import ConservativeProcess from ..base.control import Control -from ..constants import CovType, HruType, numba_num_threads, zero +from ..constants import CovType, HruType, nan, numba_num_threads, zero from ..parameters import Parameters try: @@ -119,7 +119,7 @@ def get_init_values() -> dict: "net_snow": zero, "intcp_changeover": zero, "intcp_evap": zero, - "intcp_form": 0, # could make boolean but would have to make the RAIN/SNOW match + "intcp_form": nan, # = RAIN could make boolean but would have to make the RAIN/SNOW match "intcp_stor": zero, "intcp_transp_on": 0, # could make boolean "hru_intcpevap": zero, @@ -290,6 +290,7 @@ def _calculate(self, time_length): if self._calc_method.lower() != "fortran": ( self.intcp_evap[:], + self.intcp_form[:], self.intcp_stor[:], self.net_rain[:], self.net_snow[:], @@ -310,7 +311,6 @@ def _calculate(self, time_length): hru_snow=self.hru_snow, intcp_changeover=self.intcp_changeover, intcp_evap=self.intcp_evap, - intcp_form=self.intcp_form, intcp_stor=self.intcp_stor, intcp_transp_on=self.intcp_transp_on, net_ppt=self.net_ppt, @@ -341,6 +341,7 @@ def _calculate(self, time_length): else: ( self.intcp_evap[:], + self.intcp_form[:], self.intcp_stor[:], self.net_rain[:], self.net_snow[:], @@ -360,7 +361,6 @@ def _calculate(self, time_length): self.hru_snow, self.intcp_changeover, self.intcp_evap, - self.intcp_form, self.intcp_stor, self.intcp_transp_on, self.net_ppt, @@ -407,7 +407,6 @@ def _calculate_numpy( hru_snow, intcp_changeover, intcp_evap, - intcp_form, intcp_stor, intcp_transp_on, net_ppt, @@ -440,6 +439,8 @@ def _calculate_numpy( # actual inputs # Keep the f90 call signature consistent with the args in # python/numba. + + intcp_form = np.full_like(hru_rain, np.nan, dtype="int32") for i in prange(nhru): netrain = hru_rain[i] netsnow = hru_snow[i] @@ -596,6 +597,7 @@ def _calculate_numpy( return ( intcp_evap, + intcp_form, intcp_stor, net_rain, net_snow, diff --git a/pywatershed/hydrology/prms_canopy.pyf b/pywatershed/hydrology/prms_canopy.pyf index e962ce95..c6d7f6aa 100644 --- a/pywatershed/hydrology/prms_canopy.pyf +++ b/pywatershed/hydrology/prms_canopy.pyf @@ -6,7 +6,7 @@ python module prms_canopy_f ! in module canopy ! in :prms_canopy_f:prms_canopy.f90 real(kind=8), parameter,optional :: zero=0.0 real(kind=8), parameter,optional :: one=1.0 - subroutine calc_canopy(nhru,cov_type,covden_sum,covden_win,hru_intcpstor,hru_intcpevap,hru_ppt,hru_rain,hru_snow,intcp_changeover,intcp_evap,intcp_form_in,intcp_stor,intcp_transp_on,net_ppt,net_rain,net_snow,pkwater_ante,potet,potet_sublim,snow_intcp,srain_intcp,transp_on,wrain_intcp,time_length,hru_type,nearzero,dnearzero,baresoil,grasses,land,lake,rain,snow,off,active,intcp_evap_out,intcp_stor_out,net_rain_out,net_snow_out,net_ppt_out,hru_intcpstor_out,hru_intcpevap_out,intcp_changeover_out,intcp_transp_on_out) ! in :prms_canopy_f:prms_canopy.f90:canopy + subroutine calc_canopy(nhru,cov_type,covden_sum,covden_win,hru_intcpstor,hru_intcpevap,hru_ppt,hru_rain,hru_snow,intcp_changeover,intcp_evap,intcp_stor,intcp_transp_on,net_ppt,net_rain,net_snow,pkwater_ante,potet,potet_sublim,snow_intcp,srain_intcp,transp_on,wrain_intcp,time_length,hru_type,nearzero,dnearzero,baresoil,grasses,land,lake,rain,snow,off,active,intcp_evap_out,intcp_form_out,intcp_stor_out,net_rain_out,net_snow_out,net_ppt_out,hru_intcpstor_out,hru_intcpevap_out,intcp_changeover_out,intcp_transp_on_out) ! in :prms_canopy_f:prms_canopy.f90:canopy !integer(kind=4) intent(in) :: nhru integer, optional,intent(in),check(len(cov_type)>=nhru),depend(cov_type) :: nhru=len(cov_type) integer(kind=8) dimension(nhru),intent(in) :: cov_type @@ -19,7 +19,6 @@ python module prms_canopy_f ! in real(kind=8) dimension(nhru),intent(in) :: hru_snow real(kind=8) dimension(nhru),intent(in) :: intcp_changeover real(kind=8) dimension(nhru),intent(in) :: intcp_evap - integer(kind=4) dimension(nhru),intent(in) :: intcp_form_in real(kind=8) dimension(nhru),intent(in) :: intcp_stor integer(kind=4) dimension(nhru),intent(in) :: intcp_transp_on real(kind=8) dimension(nhru),intent(in) :: net_ppt @@ -45,6 +44,7 @@ python module prms_canopy_f ! in integer(kind=4) intent(in) :: off integer(kind=4) intent(in) :: active real(kind=8) dimension(nhru),intent(out) :: intcp_evap_out + integer(kind=4) dimension(nhru),intent(out) :: intcp_form_out real(kind=8) dimension(nhru),intent(out) :: intcp_stor_out real(kind=8) dimension(nhru),intent(out) :: net_rain_out real(kind=8) dimension(nhru),intent(out) :: net_snow_out diff --git a/pywatershed/hydrology/prms_channel.py b/pywatershed/hydrology/prms_channel.py index e158c4d4..76d0b734 100644 --- a/pywatershed/hydrology/prms_channel.py +++ b/pywatershed/hydrology/prms_channel.py @@ -227,8 +227,13 @@ def _initialize_channel_data(self) -> None: ) # JLM: This is a bad idea and should throw an error rather than edit # inputs in place during run + # should also be done before computing velocity + mask_too_flat = self.seg_slope < 1e-7 + if mask_too_flat.any(): + msg = "seg_slope < 1.0e-7, set to 1.0e-4" + warn(msg, UserWarning) self.seg_slope = np.where( - self.seg_slope < 1e-7, 0.0001, self.seg_slope + self.seg_slope < 1e-7, 1.0e-4, self.seg_slope ) # not in prms6 # initialize Kcoef to 24.0 for segments with zero velocities diff --git a/pywatershed/hydrology/prms_runoff.py b/pywatershed/hydrology/prms_runoff.py index 62565615..2fe7a12b 100644 --- a/pywatershed/hydrology/prms_runoff.py +++ b/pywatershed/hydrology/prms_runoff.py @@ -42,8 +42,8 @@ class PRMSRunoff(ConservativeProcess): control: a Control object discretization: a discretization of class Parameters parameters: a parameter object of class Parameters - soil_moist_prev: Previous storage of capillary reservoir for each - HRU + soil_lower_prev: Previous storage of lower reservoir for each HRU + soil_rechr_prev: Previous storage of recharge reservoir for each HRU net_ppt: Precipitation (rain and/or snow) that falls through the canopy for each HRU net_rain: Rain that falls through canopy for each HRU @@ -72,7 +72,8 @@ def __init__( control: Control, discretization: Parameters, parameters: Parameters, - soil_moist_prev: adaptable, + soil_lower_prev: adaptable, + soil_rechr_prev: adaptable, net_ppt: adaptable, net_rain: adaptable, net_snow: adaptable, @@ -164,7 +165,8 @@ def get_parameters() -> tuple: @staticmethod def get_inputs() -> tuple: return ( - "soil_moist_prev", + "soil_lower_prev", + "soil_rechr_prev", "net_rain", "net_ppt", "net_snow", @@ -197,7 +199,7 @@ def get_init_values() -> dict: "hru_impervstor_change": zero, "dprst_vol_frac": zero, "dprst_vol_clos": zero, - "dprst_vol_open": nan, + "dprst_vol_open": zero, "dprst_vol_clos_frac": zero, "dprst_vol_open_frac": zero, "dprst_area_clos": zero, @@ -457,7 +459,8 @@ def _calculate(self, time_length, vectorized=False): potet=self.potet, snow_evap=self.snow_evap, hru_intcpevap=self.hru_intcpevap, - soil_moist_prev=self.soil_moist_prev, + soil_lower_prev=self.soil_lower_prev, + soil_rechr_prev=self.soil_rechr_prev, soil_moist_max=self.soil_moist_max, carea_max=self.carea_max, smidx_coef=self.smidx_coef, @@ -543,7 +546,8 @@ def _calculate_numpy( potet, snow_evap, hru_intcpevap, - soil_moist_prev, + soil_lower_prev, + soil_rechr_prev, soil_moist_max, carea_max, smidx_coef, @@ -601,27 +605,30 @@ def _calculate_numpy( ): dprst_chk = 0 infil[:] = 0.0 + + soil_moist_prev = soil_lower_prev + soil_rechr_prev + for k in prange(nhru): # TODO: remove duplicated vars # TODO: move setting constants outside the loop. # cdl i = Hru_route_order(k) i = k - runoff = 0.0 + runoff = zero hruarea = hru_area[i] perv_area = hru_perv[i] perv_frac = hru_frac_perv[i] - srp = 0.0 - sri = 0.0 - hru_sroffp[i] = 0.0 - contrib_fraction[i] = 0.0 + srp = zero + sri = zero + hru_sroffp[i] = zero + contrib_fraction[i] = zero hruarea_imperv = hru_imperv[i] - imperv_frac = 0.0 - if hruarea_imperv > 0.0: + imperv_frac = zero + if hruarea_imperv > zero: imperv_frac = hru_percent_imperv[i] - hru_sroffi[i] = 0.0 - imperv_evap[i] = 0.0 - hru_impervevap[i] = 0.0 + hru_sroffi[i] = zero + imperv_evap[i] = zero + hru_impervevap[i] = zero avail_et = potet[i] - snow_evap[i] - hru_intcpevap[i] availh2o = intcp_changeover[i] + net_rain[i] @@ -731,7 +738,7 @@ def _calculate_numpy( # runoff and srunoff # Compute runoff for pervious and impervious area, and depression # storage area - srunoff = 0.0 + srunoff = zero if hru_type[i] == LAND: runoff = runoff + srp * perv_area + sri * hruarea_imperv srunoff = runoff / hruarea @@ -926,10 +933,13 @@ def compute_infil( # snowpack was small and was lost to sublimation. # if net_snow < NEARZERO and net_rain > 0.0: if cond6 and through_rain > 0.0: + ## if net_snow < 1.0e-6 and net_rain > 0.0: # cond3 & cond4 & cond6 & cond1 # this is through_rain's top/most narrow case + avail_water = avail_water + through_rain infil = infil + through_rain + # double_counting += 1 # if double_counting > 1: # print("cond4") @@ -1066,9 +1076,9 @@ def dprst_comp( dprst_srp = dprst_srp + dprst_srp_clos / hru_area dprst_vol_clos = dprst_vol_clos + dprst_srp_clos srp = srp - dprst_srp / perv_frac - if srp < 0.0: - if srp < -NEARZERO: - srp = 0.0 + if srp < zero: + srp = zero + # if srp < -NEARZERO: this should raise an error if sri > 0.0: tmp = sri * imperv_frac * sro_to_dprst_imperv * hru_area diff --git a/pywatershed/hydrology/prms_soilzone.py b/pywatershed/hydrology/prms_soilzone.py index b293f88e..9c49a199 100644 --- a/pywatershed/hydrology/prms_soilzone.py +++ b/pywatershed/hydrology/prms_soilzone.py @@ -168,7 +168,6 @@ def get_init_values() -> dict: "soil_lower_ratio": zero, "soil_lower_max": nan, # completely set later "soil_moist": nan, # sm_climateflow - "soil_moist_prev": zero, # sm_climateflow "soil_moist_tot": nan, # completely set later "soil_rechr": nan, # sm_climateflow "soil_rechr_change": nan, # sm_climateflow @@ -214,6 +213,7 @@ def _set_initial_conditions(self): # changes/checks on params should throw errors so # parameter values remain the responsibility of the users # and their deficiencies are transparent? + self.hru_area_imperv = self.hru_percent_imperv * self.hru_area self.hru_area_perv = self.hru_area - self.hru_area_imperv @@ -287,32 +287,72 @@ def _set_initial_conditions(self): # JLM: some of these should just be runtime/value errors # JLM: These are for "ACTIVE and non-lake" hrus.... # JLM check that. - self.soil_moist_max = np.where( - self.soil_moist_max < 0.00001, 0.00001, self.soil_moist_max - ) - self.soil_rechr_max = np.where( - self.soil_rechr_max < 0.00001, 0.00001, self.soil_rechr_max - ) + + mask = self.soil_moist_max < 1.0e-5 + if mask.any(): + msg = "soil_moist_max < 1.0e-5, set to 1.0e-5" + warn(msg, UserWarning) + self.soil_moist_max = np.where(mask, 1.0e-5, self.soil_moist_max) + + mask = self.soil_rechr_max < 1.0e-5 + if mask.any(): + msg = "soil_rechr_max < 1.0e-5, set to 1.0e-5" + warn(msg, UserWarning) + self.soil_rechr_max = np.where(mask, 1.0e-5, self.soil_rechr_max) + + mask = self.soil_rechr_max > self.soil_moist_max + if mask.any(): + msg = ( + "soil_rechr_max > soil_moist_max, " + "soil_rechr_max set to soil_moist_max" + ) + warn(msg, UserWarning) self.soil_rechr_max = np.where( - self.soil_rechr_max > self.soil_moist_max, + mask, self.soil_moist_max, self.soil_rechr_max, ) + + mask = self.soil_rechr > self.soil_rechr_max + if mask.any(): + msg = ( + "soil_rechr_init > soil_rechr_max, " + "setting soil_rechr_init to soil_rechr_max" + ) + warn(msg, UserWarning) self.soil_rechr = np.where( - self.soil_rechr > self.soil_rechr_max, + mask, self.soil_rechr_max, self.soil_rechr, ) + + mask = self.soil_moist > self.soil_moist_max + if mask.any(): + msg = ( + "soil_moist_init > soil_moist_max, " + "setting soil_moist to soil_moist max" + ) + warn(msg, UserWarning) self.soil_moist = np.where( - self.soil_moist > self.soil_moist_max, + mask, self.soil_moist_max, self.soil_moist, ) - self.soil_rechr = np.where( - self.soil_rechr > self.soil_moist, self.soil_moist, self.soil_rechr - ) + + mask = self.soil_rechr > self.soil_moist + if mask.any(): + msg = "soil_rechr > soil_moist, setting soil_rechr to soil_moist" + warn(msg, UserWarning) + self.soil_rechr = np.where(mask, self.soil_moist, self.soil_rechr) + + mask = self.ssres_stor > self._sat_threshold + if mask.any(): + msg = ( + "ssres_stor > _sat_threshold, " + "setting ssres_stor to _sat_threshold" + ) self.ssres_stor = np.where( - self.ssres_stor > self._sat_threshold, + mask, self._sat_threshold, self.ssres_stor, ) @@ -448,7 +488,6 @@ def _calculate(self, simulation_time): self.cap_waterin[:], self.soil_moist[:], self.soil_rechr[:], - self.soil_moist_prev[:], self.hru_actet[:], self.cap_infil_tot[:], self.slow_stor[:], @@ -532,7 +571,6 @@ def _calculate(self, simulation_time): soil_lower_ratio=self.soil_lower_ratio, soil_moist=self.soil_moist, soil_moist_max=self.soil_moist_max, - soil_moist_prev=self.soil_moist_prev, soil_moist_tot=self.soil_moist_tot, soil_rechr=self.soil_rechr, soil_rechr_change=self.soil_rechr_change, @@ -591,7 +629,6 @@ def _calculate_numpy( slow_stor, soil_moist, soil_moist_max, - soil_moist_prev, soil_rechr, soil_to_gw, soil_to_ssr, @@ -661,8 +698,9 @@ def _calculate_numpy( _snow_free[:] = one - snowcov_area - # Do this here and not in advance as this is not an individual storage - soil_moist_prev[:] = soil_moist + # we dont track soil_moist_prev as it's not prognostic + # soil_moist_prev = soil_rechr and soil_lower + # soil_moist_prev[:] = soil_moist # JLM: ET calculations to be removed from soilzone. hru_actet = hru_impervevap + hru_intcpevap + snow_evap @@ -1005,7 +1043,6 @@ def _calculate_numpy( cap_waterin, soil_moist, soil_rechr, - soil_moist_prev, hru_actet, cap_infil_tot, slow_stor, diff --git a/pywatershed/utils/__init__.py b/pywatershed/utils/__init__.py index 1b713d7b..c41bcf0b 100644 --- a/pywatershed/utils/__init__.py +++ b/pywatershed/utils/__init__.py @@ -1,7 +1,7 @@ from .cbh_utils import cbh_file_to_netcdf -from .control import ControlVariables +from .control import ControlVariables, compare_control_files from .csv_utils import CsvFile -from .netcdf_utils import NetCdfCompare, NetCdfRead, NetCdfWrite +from .netcdf_utils import NetCdfRead, NetCdfWrite from .prms5_file_util import PrmsFile from .prms5util import ( Soltab, diff --git a/pywatershed/utils/control.py b/pywatershed/utils/control.py index bde156c7..524e2f74 100644 --- a/pywatershed/utils/control.py +++ b/pywatershed/utils/control.py @@ -1,5 +1,8 @@ import pathlib as pl from typing import Union +from warnings import warn + +import numpy as np from .prms5_file_util import PrmsFile @@ -54,3 +57,54 @@ def load(control_file: fileish) -> "ControlVariables": return ControlVariables( PrmsFile(control_file, file_type="control").get_data() ) + + +def compare_control_files(file0, file1, silent=False): + """Compare the contents of two control files. + + Args: + file0: str or pathlib.Path for the first PRMS control file + file1: str or pathlib.Path for the second PRMS control file + silent: suppress informative warnings (if you only want the return + value) + + Returns: + A bash-style integer where zero indicates the files are identical and + all other values indicate the files are not identical and the total + number of keys with differences. + + """ + + return_code = 0 + + ctl0 = ControlVariables.load(file0).control + ctl1 = ControlVariables.load(file1).control + + ctl0_keys_only = set(ctl0.keys()) - set(ctl1.keys()) + ctl1_keys_only = set(ctl1.keys()) - set(ctl0.keys()) + + if len(ctl0_keys_only): + return_code = len(ctl0_keys_only) + if not silent: + warn(f"Keys only in control file0: {sorted(ctl0_keys_only)}") + if len(ctl1_keys_only): + return_code = len(ctl1_keys_only) + if not silent: + warn(f"Keys only in control file1: {sorted(ctl1_keys_only)}") + + common_keys = set(ctl0.keys()).intersection(set(ctl1.keys())) + + for ck in common_keys: + try: + np.testing.assert_equal(ctl0[ck], ctl1[ck]) + except AssertionError: + return_code += 1 + if not silent: + msg = ( + f"Keys '{ck}' differ\n" + f"file0: {ctl0[ck]}\n" + f"file1: {ctl1[ck]}\n" + ) + warn(msg) + + return return_code diff --git a/pywatershed/utils/netcdf_utils.py b/pywatershed/utils/netcdf_utils.py index 6ca6744e..fafa3d01 100644 --- a/pywatershed/utils/netcdf_utils.py +++ b/pywatershed/utils/netcdf_utils.py @@ -375,7 +375,7 @@ def __init__( clobber: bool = True, zlib: bool = True, complevel: int = 4, - chunk_sizes: dict = {"time": 30, "hruid": 0}, + chunk_sizes: dict = {"time": 1, "hruid": 0}, ): if isinstance(variables, dict): group_variables = [] @@ -611,91 +611,3 @@ def add_all_data( self.variables[name][:, :] = data[:, :] return - - -class NetCdfCompare(Accessor): - """Compare variables in two NetCDF files - - Args: - base_pth: path to base NetCDF file - compare_pth: path to NetCDF file that will be compared to - the base NetCDF file - verbose: boolean determining if information should be written - to stdout - """ - - def __init__( - self, - base_pth: fileish, - compare_pth: fileish, - verbose: bool = False, - ) -> None: - self._base = NetCdfRead(base_pth) - self._compare = NetCdfRead(compare_pth) - self._verbose = verbose - if not self.__validate_comparison: - raise KeyError( - f"Variables in {compare_pth} do not exist in {base_pth}" - ) - - def __del__(self): - if self._base.dataset.isopen(): - self._base.dataset.close() - if self._compare.dataset.isopen(): - self._compare.dataset.close() - - def compare( - self, - atol: float = ATOL, - ) -> (bool, dict): - """ - - Args: - atol: absolute tolerance to use in numpy allclose (default - is single precision machine precisions ~1e-7) - - Returns: - success: boolean indicating if comparisons for all variables - in the NetCDF file passed - failures: dictionary with the maximum absolute difference for - each variable in the NetCDF file. - - """ - return self.__compare(atol=atol) - - def __compare( - self, - atol: float = ATOL, - ): - success = True - failures = {} - for variable in self._compare.variables: - base_arr = self._base.get_data(variable) - compare_arr = self._compare.get_data(variable) - if not np.allclose(compare_arr, base_arr, atol=atol): - success = False - diff = np.abs(compare_arr - base_arr) - failures[variable] = ( - np.max(diff), - atol, - np.argmax(np.max(diff, axis=0)), - ) - return success, failures - - @property - def __validate_comparison(self) -> bool: - """Validate the variables in the comparison NetCDF file relative - to the base NetCDF file - - - Returns: - valid: boolean indicating if all variables in the comparison - NetCDF file are in the base NetCDF file - - """ - valid = False - for variable in self._compare.variables: - if variable in self._base.variables: - valid = True - break - return valid diff --git a/test_data/common/control.multi_hru b/test_data/common/control.multi_hru index effd9667..60756075 100644 --- a/test_data/common/control.multi_hru +++ b/test_data/common/control.multi_hru @@ -238,7 +238,7 @@ nhruOutON_OFF nhruOutVars 1 1 -83 +102 #### nhruOut_freq 1 @@ -251,16 +251,27 @@ nhruOut_format 4 #### nhruOutVar_names -83 +102 4 albedo cap_infil_tot +cap_waterin contrib_fraction +dprst_area_clos +dprst_area_clos_max +dprst_area_open +dprst_area_open_max dprst_seep_hru dprst_evap_hru dprst_insroff_hru dprst_sroff_hru dprst_stor_hru +dprst_vol_clos +dprst_vol_clos_frac +dprst_vol_frac +dprst_vol_open +dprst_vol_open_frac +dunnian_flow freeh2o gwres_flow gwres_in @@ -278,9 +289,12 @@ hru_snow hru_sroffi hru_sroffp hru_storage +imperv_evap +imperv_stor infil intcp_changeover intcp_evap +intcp_form intcp_stor iso net_ppt @@ -302,8 +316,11 @@ potet_rechr pptmix pptmix_nopack pref_flow +pref_flow_in pref_flow_infil +pref_flow_max pref_flow_stor +pref_flow_thrsh prmx pst recharge @@ -325,6 +342,7 @@ ssres_flow ssres_in ssres_stor ssr_to_gw +swale_actet swrad tavgc tavgf @@ -334,6 +352,7 @@ tmaxf tminc tminf transp_on +unused_potet hru_outflow hru_streamflow_out #### diff --git a/test_data/common/control.single_hru b/test_data/common/control.single_hru index 74ea87b3..5dce61f1 100644 --- a/test_data/common/control.single_hru +++ b/test_data/common/control.single_hru @@ -238,7 +238,7 @@ nhruOutON_OFF nhruOutVars 1 1 -81 +100 #### nhruOut_freq 1 @@ -251,16 +251,27 @@ nhruOut_format 5 #### nhruOutVar_names -81 +100 4 albedo cap_infil_tot +cap_waterin contrib_fraction +dprst_area_clos +dprst_area_clos_max +dprst_area_open +dprst_area_open_max +dprst_seep_hru dprst_evap_hru dprst_insroff_hru -dprst_seep_hru -dprst_stor_hru dprst_sroff_hru +dprst_stor_hru +dprst_vol_clos +dprst_vol_clos_frac +dprst_vol_frac +dprst_vol_open +dprst_vol_open_frac +dunnian_flow freeh2o gwres_flow gwres_in @@ -278,9 +289,12 @@ hru_snow hru_sroffi hru_sroffp hru_storage +imperv_evap +imperv_stor infil intcp_changeover intcp_evap +intcp_form intcp_stor iso net_ppt @@ -302,8 +316,11 @@ potet_rechr pptmix pptmix_nopack pref_flow +pref_flow_in pref_flow_infil +pref_flow_max pref_flow_stor +pref_flow_thrsh prmx pst recharge @@ -325,6 +342,7 @@ ssres_flow ssres_in ssres_stor ssr_to_gw +swale_actet swrad tavgc tavgf @@ -334,6 +352,7 @@ tmaxf tminc tminf transp_on +unused_potet #### nhruOutBaseFileName 1 diff --git a/test_data/drb_2yr/control.test b/test_data/drb_2yr/control.test index effd9667..60756075 100644 --- a/test_data/drb_2yr/control.test +++ b/test_data/drb_2yr/control.test @@ -238,7 +238,7 @@ nhruOutON_OFF nhruOutVars 1 1 -83 +102 #### nhruOut_freq 1 @@ -251,16 +251,27 @@ nhruOut_format 4 #### nhruOutVar_names -83 +102 4 albedo cap_infil_tot +cap_waterin contrib_fraction +dprst_area_clos +dprst_area_clos_max +dprst_area_open +dprst_area_open_max dprst_seep_hru dprst_evap_hru dprst_insroff_hru dprst_sroff_hru dprst_stor_hru +dprst_vol_clos +dprst_vol_clos_frac +dprst_vol_frac +dprst_vol_open +dprst_vol_open_frac +dunnian_flow freeh2o gwres_flow gwres_in @@ -278,9 +289,12 @@ hru_snow hru_sroffi hru_sroffp hru_storage +imperv_evap +imperv_stor infil intcp_changeover intcp_evap +intcp_form intcp_stor iso net_ppt @@ -302,8 +316,11 @@ potet_rechr pptmix pptmix_nopack pref_flow +pref_flow_in pref_flow_infil +pref_flow_max pref_flow_stor +pref_flow_thrsh prmx pst recharge @@ -325,6 +342,7 @@ ssres_flow ssres_in ssres_stor ssr_to_gw +swale_actet swrad tavgc tavgf @@ -334,6 +352,7 @@ tmaxf tminc tminf transp_on +unused_potet hru_outflow hru_streamflow_out #### diff --git a/test_data/generate/conftest.py b/test_data/generate/conftest.py index 31468dd7..b0ddd77d 100644 --- a/test_data/generate/conftest.py +++ b/test_data/generate/conftest.py @@ -4,6 +4,7 @@ from fnmatch import fnmatch from platform import processor from typing import List +from warnings import warn import pytest @@ -64,16 +65,28 @@ def scheduler_active(): def enforce_scheduler(test_dir): + """Enforce the use of scheduler + + Args: + test_dir: the domain directory to run or schedule + + Return: + True if must use a scheduler, False if not + """ + if scheduler_active(): - return None + return False glob_match = list( fnmatch(str(test_dir), gg) for gg in domain_globs_schedule ) if any(glob_match): - raise RuntimeError( + msg = ( f"Domain '{test_dir}' must be scheduled (use --force to override)" ) - return None + warn(msg, RuntimeWarning) + return True + + return False def collect_simulations( @@ -91,7 +104,9 @@ def collect_simulations( # optionally enforce scheduler if not force: - enforce_scheduler(test_dir) + schedule = enforce_scheduler(test_dir) + if schedule: + continue # if control file is found, add simulation ctrl_file = next( diff --git a/test_data/generate/convert_prms_output_to_nc.py b/test_data/generate/convert_prms_output_to_nc.py index af560a6d..3f5c436b 100644 --- a/test_data/generate/convert_prms_output_to_nc.py +++ b/test_data/generate/convert_prms_output_to_nc.py @@ -1,4 +1,3 @@ -from pathlib import Path from filelock import FileLock import pytest diff --git a/test_data/generate/prms_diagnostic_variables.py b/test_data/generate/prms_diagnostic_variables.py index 12610eb3..92cccfa7 100644 --- a/test_data/generate/prms_diagnostic_variables.py +++ b/test_data/generate/prms_diagnostic_variables.py @@ -13,10 +13,10 @@ "gwres_stor": pws.PRMSGroundwater, "dprst_stor_hru": pws.PRMSRunoff, "hru_impervstor": pws.PRMSRunoff, + "hru_intcpstor": pws.PRMSCanopy, "pref_flow_stor": pws.PRMSSoilzone, "slow_stor": pws.PRMSSoilzone, "soil_lower": pws.PRMSSoilzone, - "soil_moist": pws.PRMSSoilzone, "soil_rechr": pws.PRMSSoilzone, "ssres_stor": pws.PRMSSoilzone, "freeh2o": pws.PRMSSnow, @@ -25,6 +25,9 @@ prev_rename = { "gwres_stor": "gwres_stor_old", + "hru_impervstor": "hru_impervstor_old", + "hru_intcpstor": "hru_intcpstor_old", + "dprst_stor_hru": "dprst_stor_hru_old", } change_rename = {} @@ -166,8 +169,6 @@ def diagnose_final_vars_to_nc( new variable which requires much/all of the PRMS output to be already present in netcdf output format. This is why this is a final diagnostic - Currently only: "through_rain" - Args: var_name: str name of the variable to create, not a PRMS variable. output_dir: where the netcdf file will be written, also where to look diff --git a/test_data/hru_1/control.test b/test_data/hru_1/control.test index 74ea87b3..5dce61f1 100644 --- a/test_data/hru_1/control.test +++ b/test_data/hru_1/control.test @@ -238,7 +238,7 @@ nhruOutON_OFF nhruOutVars 1 1 -81 +100 #### nhruOut_freq 1 @@ -251,16 +251,27 @@ nhruOut_format 5 #### nhruOutVar_names -81 +100 4 albedo cap_infil_tot +cap_waterin contrib_fraction +dprst_area_clos +dprst_area_clos_max +dprst_area_open +dprst_area_open_max +dprst_seep_hru dprst_evap_hru dprst_insroff_hru -dprst_seep_hru -dprst_stor_hru dprst_sroff_hru +dprst_stor_hru +dprst_vol_clos +dprst_vol_clos_frac +dprst_vol_frac +dprst_vol_open +dprst_vol_open_frac +dunnian_flow freeh2o gwres_flow gwres_in @@ -278,9 +289,12 @@ hru_snow hru_sroffi hru_sroffp hru_storage +imperv_evap +imperv_stor infil intcp_changeover intcp_evap +intcp_form intcp_stor iso net_ppt @@ -302,8 +316,11 @@ potet_rechr pptmix pptmix_nopack pref_flow +pref_flow_in pref_flow_infil +pref_flow_max pref_flow_stor +pref_flow_thrsh prmx pst recharge @@ -325,6 +342,7 @@ ssres_flow ssres_in ssres_stor ssr_to_gw +swale_actet swrad tavgc tavgf @@ -334,6 +352,7 @@ tmaxf tminc tminf transp_on +unused_potet #### nhruOutBaseFileName 1 diff --git a/test_data/ucb_2yr/control.test b/test_data/ucb_2yr/control.test index effd9667..60756075 100644 --- a/test_data/ucb_2yr/control.test +++ b/test_data/ucb_2yr/control.test @@ -238,7 +238,7 @@ nhruOutON_OFF nhruOutVars 1 1 -83 +102 #### nhruOut_freq 1 @@ -251,16 +251,27 @@ nhruOut_format 4 #### nhruOutVar_names -83 +102 4 albedo cap_infil_tot +cap_waterin contrib_fraction +dprst_area_clos +dprst_area_clos_max +dprst_area_open +dprst_area_open_max dprst_seep_hru dprst_evap_hru dprst_insroff_hru dprst_sroff_hru dprst_stor_hru +dprst_vol_clos +dprst_vol_clos_frac +dprst_vol_frac +dprst_vol_open +dprst_vol_open_frac +dunnian_flow freeh2o gwres_flow gwres_in @@ -278,9 +289,12 @@ hru_snow hru_sroffi hru_sroffp hru_storage +imperv_evap +imperv_stor infil intcp_changeover intcp_evap +intcp_form intcp_stor iso net_ppt @@ -302,8 +316,11 @@ potet_rechr pptmix pptmix_nopack pref_flow +pref_flow_in pref_flow_infil +pref_flow_max pref_flow_stor +pref_flow_thrsh prmx pst recharge @@ -325,6 +342,7 @@ ssres_flow ssres_in ssres_stor ssr_to_gw +swale_actet swrad tavgc tavgf @@ -334,6 +352,7 @@ tmaxf tminc tminf transp_on +unused_potet hru_outflow hru_streamflow_out ####