revamp testing of PRMSSoilzone; add comparison notebook with full steps

autotest/test_control_read.py

@@ -1,7 +1,8 @@
+import pathlib as pl
+
 import numpy as np
 import pytest
 
-from pywatershed.constants import __pywatershed_root__ as pws_root
 from pywatershed.utils import ControlVariables, compare_control_files
 from utils import assert_or_print
 
@@ -43,7 +44,7 @@ def test_control_read(domain, control_keys):
 
 
 def test_control_compare():
-    common_dir = pws_root.parent / "test_data/common"
+    common_dir = pl.Path("../test_data/common")
     n_diffs = compare_control_files(
         common_dir / "control.single_hru",
         common_dir / "control.multi_hru",

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 = 1.0e-8
 
 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,37 @@ 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)
+        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

examples/compare_pws_prms.ipynb

@@ -7,18 +7,18 @@
    "source": [
     "# Compare PRMS and PWS\n",
     "\n",
-    "This notebook compares PWS (pywatershed) and PRMS runs in a controlled fashion based on how PRMS/NHM domains are setup for testing in the `test_data/` directory of the repo. The domains supplied with the pywatershed repo are `test_data/hru_1`, `test_data/drb_2yr`, and `test_data/ucb_2yr`. For this notebook the following file listing of any of these domain directories gives the critical to the functioning of this notebook. \n",
+    "This notebook compares PWS (pywatershed) and PRMS runs with input based on how PRMS/NHM domains are setup for testing in the `test_data/` directory of pywtaershed. The domains supplied with the pywatershed repo are `test_data/hru_1`, `test_data/drb_2yr`, and `test_data/ucb_2yr`. For this notebook the following file listing of any of these domain directories gives the critical to the functioning of this notebook.:\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",
-    "If you want to supply your own domain (e.g. CONUS), you would need all these files to be present in `test_data/some_domain`. Note that the `*.cbh` files must be pre-converted to netcdf for `prcp`, `tmin`, and `tmax`. The `control.test` file is a bit specific for being able to run sub-models and includes a nearly maximal amount of model output (time-inefficient for both PRMS and PWS). The stock control file 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. These of course can be modified in some fashion. For running a CONUS domain, 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).\n",
+    "If you want to supply your own domain (e.g. CONUS), you need all these files to be present in `test_data/some_domain`. Note that the `*.cbh` files must be pre-converted to netcdf for `prcp`, `tmin`, and `tmax`. 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 includes 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. These, of course, can be modified. For running a CONUS 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).\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 is included so that it can be tailored to your tastes."
+    "All of the code required for plotting below is included so that it can be further tailored to your tastes."
    ]
   },
   {
@@ -143,16 +143,26 @@
     "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 --domain={domain_name} --force\",\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,
@@ -162,7 +172,7 @@
    "source": [
     "if run_pws:\n",
     "    nhm_processes = [\n",
-    "        # pws.PRMSSolarGeometry,\n",
+    "        # pws.PRMSSolarGeometry,  # submodles are possible\n",
     "        # pws.PRMSAtmosphere,\n",
     "        # pws.PRMSCanopy,\n",
     "        # pws.PRMSSnow,\n",
@@ -205,6 +215,14 @@
     "    nhm.run(finalize=True)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "aa221ffe-eb0a-4332-94fd-9826710adc74",
+   "metadata": {},
+   "source": [
+    "## Compare outputs"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -537,7 +555,7 @@
    "outputs": [],
    "source": [
     "if pws.PRMSChannel in nhm_processes:\n",
-    "    plot_proc_stats(pws.PRMSChannel, \"RMSE\", 5)"
+    "    plot_proc_stats(pws.PRMSChannel, \"RMSE\", 4)"
    ]
   },
   {

test_data/common/control.multi_hru

@@ -238,7 +238,7 @@ nhruOutON_OFF
 nhruOutVars
 1
 1
-94
+101
 ####
 nhruOut_freq
 1
@@ -251,10 +251,11 @@ nhruOut_format
 4
 ####
 nhruOutVar_names
-94
+101
 4
 albedo
 cap_infil_tot
+cap_waterin
 contrib_fraction
 dprst_area_clos
 dprst_area_clos_max
@@ -270,6 +271,7 @@ dprst_vol_clos_frac
 dprst_vol_frac
 dprst_vol_open
 dprst_vol_open_frac
+dunnian_flow
 freeh2o
 gwres_flow
 gwres_in
@@ -313,8 +315,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
@@ -336,6 +341,7 @@ ssres_flow
 ssres_in
 ssres_stor
 ssr_to_gw
+swale_actet
 swrad
 tavgc
 tavgf
@@ -345,6 +351,7 @@ tmaxf
 tminc
 tminf
 transp_on
+unused_potet
 hru_outflow
 hru_streamflow_out
 ####

test_data/common/control.single_hru

@@ -238,7 +238,7 @@ nhruOutON_OFF
 nhruOutVars
 1
 1
-92
+99
 ####
 nhruOut_freq
 1
@@ -251,25 +251,27 @@ nhruOut_format
 5
 ####
 nhruOutVar_names
-92
+99
 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
@@ -313,8 +315,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
@@ -336,6 +341,7 @@ ssres_flow
 ssres_in
 ssres_stor
 ssr_to_gw
+swale_actet
 swrad
 tavgc
 tavgf
@@ -345,6 +351,7 @@ tmaxf
 tminc
 tminf
 transp_on
+unused_potet
 ####
 nhruOutBaseFileName
 1