parameter edits

autotest/test_cbh_to_netcdf.py

@@ -5,7 +5,7 @@
 import xarray as xr
 
 from pywatershed.parameters import PrmsParameters
-from pywatershed.utils.cbh_utils import cbh_files_to_df, cbh_files_to_netcdf
+from pywatershed.utils.cbh_utils import cbh_file_to_netcdf, cbh_files_to_df
 from utils import assert_or_print
 
 var_cases = ["prcp", "rhavg", "tmax", "tmin"]
@@ -30,7 +30,7 @@ def params_and_none(request, domain):
             "prcp": 0.12495362248866715,
         },
         "df_concat": 40.7829059976932,
-        "np_dict": {
+        "file_to_netcdf": {
             "prcp": 0.12495362248866718,
             "rhavg": 62.738591579267386,
             "tmax": 60.26121597238987,
@@ -43,7 +43,7 @@ def params_and_none(request, domain):
             "prcp": 0.13502103505843072,
         },
         "df_concat": 34.968545798553144,
-        "np_dict": {
+        "file_to_netcdf": {
             "prcp": 0.13502103505843072,
             "tmax": 61.986048747913195,
             "tmin": 42.78456761268781,
@@ -55,7 +55,7 @@ def params_and_none(request, domain):
             "prcp": 0.046485653521159784,
         },
         "df_concat": 34.05471072235577,
-        "np_dict": {
+        "file_to_netcdf": {
             "prcp": 0.046485653521159805,
             "rhavg": 50.87569199962628,
             "tmax": 56.74147989347377,
@@ -79,29 +79,32 @@ def test_cbh_files_to_df(domain, var, params):
     return
 
 
-def test_cbh_files_to_netcdf(domain, params, tmp_path):
-    nc_file = tmp_path / "cbh_files_to_netcdf.nc"
+def test_cbh_file_to_netcdf(domain, params, tmp_path):
     input_files_dict = domain["cbh_inputs"]
-    _ = cbh_files_to_netcdf(input_files_dict, params, nc_file)
-
-    answers = answer_key[domain["domain_name"]]["np_dict"]
-    results_ds = xr.open_dataset(nc_file)
     results = {}
+    for var, cbh_file in input_files_dict.items():
+        nc_file = tmp_path / f"{var}.nc"
+        _ = cbh_file_to_netcdf(cbh_file, params, nc_file)
+        results[var] = xr.open_dataset(nc_file)[var]
+
+    results["time"] = results[var].time
+
+    answers = answer_key[domain["domain_name"]]["file_to_netcdf"]
     for var in answers.keys():
         if var == "time":
             results[var] = (
-                results_ds[var]
+                results[var]
                 .values.astype("datetime64[s]")
                 .astype(np.int32)
                 .mean()
             )
         else:
-            results[var] = results_ds[var].mean().mean().values.tolist()
+            results[var] = results[var].mean().mean().values.tolist()
 
     assert_or_print(
         results,
         answers,
-        "files_to_np_dict",
+        "file_to_netcdf",
         print_ans=domain["print_ans"],
         close=True,
     )

autotest/test_parameters.py

@@ -0,0 +1,32 @@
+from pprint import pprint
+
+from numpy.testing import assert_equal
+
+from pywatershed.base.parameters import _set_dict_read_write
+from pywatershed.parameters import Parameters, PrmsParameters
+from utils import assert_dicts_equal
+
+
+def test_param_dd_param(domain):
+    # round trip from read-only to read-write to read-only
+    # use a PRMS Parameter file for now
+    domain_dir = domain["param_file"].parent
+    params = Parameters.from_netcdf(
+        domain_dir / "parameters_PRMSGroundwater.nc"
+    )
+    # These both copy by default
+    param_dd = params.to_dd()
+    params2 = Parameters(**param_dd.data)
+
+    assert_dicts_equal(params.data, params2.data)
+
+    param_dd.data_vars["gwflow_coef"] *= 4
+    params3 = Parameters(**param_dd.data)
+
+    try:
+        assert_dicts_equal(params.data, params3.data)
+        assert False
+    except AssertionError:
+        pass
+
+    return

autotest/utils.py

@@ -1,3 +1,5 @@
+from types import MappingProxyType
+
 import numpy as np
 
 print_ans = False
@@ -41,3 +43,22 @@ def assert_or_print(
     # Always fail if printing answers
     assert not print_ans
     return
+
+
+def assert_dicts_equal(dic1, dic2):
+    assert set(dic1.keys()) == set(dic2.keys())
+
+    # add cases as needed
+    for kk, vv in dic1.items():
+        if isinstance(vv, (dict, MappingProxyType)):
+            assert_dicts_equal(dic1[kk], dic2[kk])
+        elif isinstance(vv, np.ndarray):
+            np.testing.assert_equal(dic1[kk], dic2[kk])
+        else:
+            if (
+                isinstance(vv, float)
+                and np.isnan(dic1[kk])
+                and np.isnan(dic2[kk])
+            ):
+                continue
+            assert dic1[kk] == dic2[kk]

doc/api/utils.rst

@@ -5,11 +5,8 @@
 Utils
 ##########
 
-Utilities to define PRMS variables, dimensions, and parameters.
-
 .. autosummary::
    :toctree: generated/
 
     ControlVariables
-    PrmsParameters
-
+    utils.cbh_file_to_netcdf

doc/whats-new.rst

@@ -19,7 +19,9 @@ v0.3.0 (unreleased)
 
 New features
 ~~~~~~~~~~~~
-- `Conda feedstock for pywatershed <https://github.com/conda-forge/staged-recipes/pull/23428>`_.
+- Example notebook of how to edit Parameters with associated bug fixes to do so.
+  (:pull:`232`) By `James McCreight <https://github.com/jmccreight>`_.
+- Conda feedstock for pywatershed `<https://github.com/conda-forge/staged-recipes/pull/23428>`_.
   By `Joseph Hughes <https://github.com/jdhughes-usgs>`_.
 
 
@@ -37,6 +39,8 @@ Performance
 
 Bug fixes
 ~~~~~~~~~
+- Resolve issues with different ways of specifiying netcdf output options.
+  (:pull:`230`) By `James McCreight <https://github.com/jmccreight>`_.
 
 
 Documentation

examples/02_prms_legacy_models.ipynb

@@ -104,17 +104,17 @@
     "    rmtree(cbh_nc_dir)\n",
     "cbh_nc_dir.mkdir(parents=True)\n",
     "\n",
-    "cbh_files = {\n",
-    "    \"prcp\": domain_dir / \"prcp.cbh\",\n",
-    "    \"tmax\": domain_dir / \"tmax.cbh\",\n",
-    "    \"tmin\": domain_dir / \"tmin.cbh\",\n",
-    "}\n",
+    "cbh_files = [\n",
+    "    domain_dir / \"prcp.cbh\",\n",
+    "    domain_dir / \"tmax.cbh\",\n",
+    "    domain_dir / \"tmin.cbh\",\n",
+    "]\n",
     "\n",
     "params = pws.parameters.PrmsParameters.load(domain_dir / \"myparam.param\")\n",
     "\n",
-    "for kk, vv in cbh_files.items():\n",
-    "    out_file = cbh_nc_dir / f\"{kk}.nc\"\n",
-    "    pws.utils.cbh_files_to_netcdf({kk: vv}, params, out_file)"
+    "for cbh_file in cbh_files:\n",
+    "    out_file = cbh_nc_dir / cbh_file.with_suffix(\".nc\").name\n",
+    "    pws.utils.cbh_file_to_netcdf(cbh_file, params, out_file)"
    ]
   },
   {
@@ -437,17 +437,10 @@
    "outputs": [],
    "source": [
     "submodel = pws.Model(\n",
-    "    submodel_processes,\n",
+    "    [pws.PRMSSoilzone, pws.PRMSGroundwater, pws.PRMSChannel],\n",
     "    control=control,\n",
     "    parameters=params,\n",
-    ")\n",
-    "\n",
-    "pws.analysis.ModelGraph(\n",
-    "    submodel,\n",
-    "    hide_variables=not show_params,\n",
-    "    show_params=show_params,\n",
-    "    process_colors=palette,\n",
-    ").SVG(verbose=True, dpi=48)"
+    ")"
    ]
   },
   {
@@ -468,19 +461,6 @@
     "%%time\n",
     "submodel.run(finalize=True)"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "bd25d1a7-3ca3-4d5e-9623-e5c0d3f20b53",
-   "metadata": {
-    "editable": true,
-    "slideshow": {
-     "slide_type": ""
-    }
-   },
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {

examples/param_edits.ipynb

@@ -0,0 +1,143 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "f932e1c0-24ac-4be8-a886-8554cd1e4412",
+   "metadata": {},
+   "source": [
+    "# Parameter Editing\n",
+    "This notebook gives a quick recipe for how to do calibration or sensitivity analysis with pywatershed. \n",
+    "It is a design feature that parameters, more specifically the Parameter class, are read-only because \n",
+    "it should be the case that parameters supplied are used and the code is not opaquely modifying these.\n",
+    "\n",
+    "As a consequence, one has to make the Parameter class editable. Below this is accomplished by doing\n",
+    "`the_parameters.to_dd()` which returns a DatasetDict which is editable. One has to know something about\n",
+    "how DatasetDicts are constructed to edit effectively, information can be found in the [documentation](https://pywatershed.readthedocs.io/en/main/api/generated/pywatershed.base.DatasetDict.html#pywatershed.base.DatasetDict). \n",
+    "The edited DatasetDict can be made a Parameters object again by `Parameters(**param_dict.data)`, as shown below. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "73380b42-3336-4c6c-99cb-f797841a132a",
+   "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": "769aa543-dfa5-43e9-84f6-0830ad0cb1d0",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    }
+   },
+   "outputs": [],
+   "source": [
+    "import pathlib as pl\n",
+    "from pprint import pprint\n",
+    "\n",
+    "import numpy as np\n",
+    "import pywatershed as pws\n",
+    "import xarray as xr"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "01da8dca-622f-44c2-9041-8a43826d5d6f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "domain_dir = pws.constants.__pywatershed_root__ / \"data/drb_2yr\"\n",
+    "nb_output_dir = pl.Path(\"./param_edits\")\n",
+    "nb_output_dir.mkdir(exist_ok=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8b08b96b-9bcc-4b1b-a2e9-231949b2dd0c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# A legacy PRMS parameter file\n",
+    "params = pws.parameters.PrmsParameters.load(domain_dir / \"myparam.param\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cd179039-a91d-4057-8e11-4ba49cb387bd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "param_list = []\n",
+    "param_files = []\n",
+    "for ii in range(11):\n",
+    "    param_dict = params.to_dd()  # copies by default\n",
+    "    multiplier = ii * 0.05 + 0.75\n",
+    "    print(\"multiplier = \", multiplier)\n",
+    "    param_dict.data_vars[\"K_coef\"] *= multiplier\n",
+    "    param_file_name = nb_output_dir / f\"perturbed_params_{str(ii).zfill(3)}.nc\"\n",
+    "    param_files += [param_file_name]\n",
+    "    # These could avoid export to netcdf4 if just using in memory\n",
+    "    # could store in a list like: param_list.append(pws.Parameters(**param_dict.data))\n",
+    "    pws.Parameters(**param_dict.data).to_netcdf(\n",
+    "        param_file_name, use_xr=True\n",
+    "    )  # using xarray, more work necessary for nc4 export"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "04cb7d5a-4979-4859-872c-6ccd6f6f1fc3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# this provides a check that the values from file are what we expect\n",
+    "for ff in param_files:\n",
+    "    # the problem arises on the read with xarray default decoding\n",
+    "    # but we can just open the netcdf file as Parameters\n",
+    "    # ds = xr.open_dataset(ff, decode_times=False, decode_timedelta=False)\n",
+    "    # k_coef = ds[\"K_coef\"]\n",
+    "    new_params = pws.Parameters.from_netcdf(ff)\n",
+    "    k_coef = new_params.data_vars[\"K_coef\"]\n",
+    "    multipliers = k_coef / params.data_vars[\"K_coef\"]\n",
+    "    assert (multipliers - multipliers[0] < 1e-15).all()\n",
+    "    print(multipliers[0])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8dc58917-2b29-474a-bb56-1d401c68b4ef",
+   "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
+}