feat: Add test suite, unit tests on PR, and weekly live-service tests

.github/workflows/canary.yml

@@ -0,0 +1,46 @@
+name: canary
+
+# Weekly canary run.
+#  1. Installs the latest releases of all (unpinned) dependencies and runs the
+#     unit suite, so a breaking dependency release (e.g. rosetta-soil changing
+#     output conventions) is flagged.
+#  2. Makes small live requests against the external services pytRIBS depends
+#     on (ISRIC SoilGrids, SOLUS, POLARIS, NASA Giovanni).
+#
+# The Giovanni test needs Earthdata Login credentials stored as repository
+# secrets EARTHDATA_USERNAME and EARTHDATA_PASSWORD.
+
+on:
+  schedule:
+    - cron: '0 12 * * 1'  # Mondays 12:00 UTC
+  workflow_dispatch:
+
+jobs:
+  canary:
+    name: Weekly dependency and live-service canary
+    runs-on: ubuntu-latest
+    timeout-minutes: 45
+
+    steps:
+      - name: Checkout repository
+        uses: actions/checkout@v5
+
+      - name: Set up Python 3.13
+        uses: actions/setup-python@v6
+        with:
+          python-version: "3.13"
+
+      - name: Install pytRIBS with latest dependency releases
+        run: |
+          python -m pip install --upgrade pip
+          pip install -e ".[test]"
+          pip list
+
+      - name: Run unit tests against latest dependencies
+        run: pytest -v
+
+      - name: Run live-service canaries
+        env:
+          EARTHDATA_USERNAME: ${{ secrets.EARTHDATA_USERNAME }}
+          EARTHDATA_PASSWORD: ${{ secrets.EARTHDATA_PASSWORD }}
+        run: pytest -m network -v

.github/workflows/install_and_test.yml

@@ -33,7 +33,10 @@ jobs:
       - name: Install pytRIBS
         run: |
           python -m pip install --upgrade pip
-          pip install -e .
+          pip install -e ".[test]"
 
       - name: Smoke test import
         run: python -c "import pytRIBS; print('pytRIBS imported successfully')"
+
+      - name: Run unit tests
+        run: pytest -v

CHANGELOG.md

@@ -22,9 +22,11 @@ The v1.0.0 changes listed below are abbreviated. For specific details refer to t
 * **SOLUS100 Soil Data Source:** Added `'SOLUS'` as a third source option in `run_soil_workflow()` alongside ISRIC and POLARIS. New methods `get_solus_grids()` and `get_solus_bedrock()` download the USDA SOLUS100 (100 m) soil property and depth-to-bedrock rasters via Cloud Optimized GeoTIFF windowed reads, fetching only the watershed extent. SOLUS parameters are estimated with ROSETTA model code 3 (sand, silt, clay, bulk density) through the same prediction and van Genuchten to Brooks-Corey conversion path as ISRIC. The workflow scales SOLUS bulk density from its delivered g/cm3 x 100 encoding, and repairs SOLUS non-soil pixels (mapped rock outcrop and water, encoded as undeclared zeros) by interpolating from neighboring soil pixels and re-normalizing textures, so parameter grids and the soil classification map are continuous. The non-soil footprint is preserved as `nonsoil_mask.asc` with a printed warning, since interpolated soil properties potentially misrepresent rock outcrop cells unless the soil column is also constrained. ([#51](https://github.com/tRIBS-Model/pytRIBS/pull/51))
 * **Bedrock-Based Initial Groundwater:** New `generate_initial_groundwater()` method computes a spatially variable initial water table by scaling a bedrock depth raster (e.g. from `get_solus_bedrock()`) by a user-specified fraction, replacing guesswork uniform initializations for spin-up runs. ([#51](https://github.com/tRIBS-Model/pytRIBS/pull/51))
 * **Ks Decay Floor:** `compute_ks_decay()` gained a `min_f` parameter (default 1e-4) applied after curve fitting, preventing near-zero decay values that tRIBS cannot handle numerically where Ks is nearly uniform with depth. ([#51](https://github.com/tRIBS-Model/pytRIBS/pull/51))
+* **Test Suite:** Added a pytest suite covering the major workflows (ROSETTA soil parameters with known-answer checks, soil classification, NLDAS unit conversions, file format round trips, evaluation metrics, and class construction), run on every pull request. A new weekly workflow additionally tests against the latest dependency releases and makes small live requests to the external data services (ISRIC, SOLUS, POLARIS, NASA Giovanni) to catch upstream API changes early. ([#52](https://github.com/tRIBS-Model/pytRIBS/pull/52))
 
 ### Fixed
 * **ROSETTA Soil Parameter Grids:** Fixed two stacked issues in `process_raw_soil()`. First, `rosetta-soil` 0.32 (released 03/2026) changed `rosetta()` to return `alpha`, `n`, and `Ksat` in linear units rather than log10; pytRIBS's back-transform then silently corrupted every derived grid in any environment installed after that release. The workflow now requests geometric-mean estimates and the dependency is pinned to `rosetta-soil>=0.3`. Second, the Brooks-Corey air-entry pressure `psib` was previously approximated as `-1/alpha`, which overestimates the value depending on soil type, it is now converted from the van Genuchten parameters following Morel-Seytoux et al. (1996, WRR 32(5)), preserving the effective capillary drive so infiltration behavior is insensitive to the retention model choice. Pixels are also now predicted in large batched calls rather than one `rosetta()` call per pixel, substantially speeding up the soil workflow. ([#51](https://github.com/tRIBS-Model/pytRIBS/pull/51))
+* **Class Construction from Input Files:** Constructing the `Soil`, `Land`, `Met`, or `Mesh` classes with the `input_file` argument raised a `TypeError`; all four constructors now read the input file the same way as `Model` and `Results`. ([#52](https://github.com/tRIBS-Model/pytRIBS/pull/52))
 
 ### Changed & Refactored
 * **Python 3.11+ Required:** The minimum supported Python version is raised from 3.10 to 3.11, driven by `rosetta-soil` >= 0.3 requiring 3.11. Python 3.10 reaches end-of-life in October 2026. ([#51](https://github.com/tRIBS-Model/pytRIBS/pull/51))

PR_DESCRIPTION.md

@@ -0,0 +1,41 @@
+This pull request adds a test suite and continuous-integration coverage to pytRIBS. Until now the GitHub Actions workflow only installed the package and imported it. Multiple failure modes are now covered by tests, in two tiers: fast unit tests that run on every pull request, and scheduled live-service tests that run weekly against the real external APIs.
+
+### Unit tests (run on every PR)
+
+* **`test_rosetta_soil_params.py`:** direct regression guard. Feeds four textbook textures (sand, loam, clay, silt loam) through `_rosetta_to_tribs_params` and asserts (1) literature-informed plausibility windows on Ks per texture, physical bounds on theta_r/theta_s/psib/m, and correct texture ordering and (2) "golden values" frozen from `rosetta-soil` 0.3.x at a loose ±25% tolerance to flag subtler prediction drift. Also verifies the ROSETTA model-code selection.
+* **`test_soil_map.py`:** USDA texture classification through `create_soil_map`, covering both the ISRIC (g/kg) and SOLUS (% mass) input conventions and the sequential class renumbering.
+* **`test_ks_decay.py`:** generates synthetic Ks-with-depth grids from a known decay parameter and asserts `compute_ks_decay` recovers it, floors uniform profiles at `min_f`, and propagates nodata.
+* **`test_met_conversions.py`:** hand-checked "golden values" for the NLDAS unit conversions (Pa to hPa, K to degC, specific humidity to RH, 10 m to 2 m wind log-profile scaling), the GMT timestamp shift, `.mdf`/`.sdf` round trips, and the observation-station validation errors.
+* **`test_soil_tables.py` / `test_ascii_io.py`:** round trips for the `.sdt` soil table and the ESRI ASCII raster writer/reader, including NaN-to-nodata replacement and cell-size preservation.
+* **`test_evaluate_metrics.py`:** known-answer checks for NSE, KGE, RMSE, and percent bias.
+* **`test_classes.py`:** constructs every top-level class bare and from a written `.in` file, locking in the constructor fix and the input-file write/read round trip.
+
+### Live-service canaries (weekly, `pytest -m network`)
+
+Small real requests with plausibility gates on the responses, so an upstream API change is flagged by a failed-workflow email instead of by a broken user run:
+
+* ISRIC SoilGrids WCS, SOLUS100 COG windowed reads, and POLARIS tile downloads over a ~2×2 km test domain, each asserting a readable, non-empty grid with in-range values.
+* NASA Giovanni NLDAS point timeseries, asserting all seven forcing variables return with physically plausible values. Because `get_nldas_point` swallows per-variable fetch errors, a missing column is the failure signal.
+
+### Continuous integration
+
+* **`install_and_test.yml`:** now installs `.[test]` and runs `pytest` (network tests excluded by marker) after the import smoke test, on Python 3.11 and 3.13. Because dependencies are unpinned, every PR run also exercises the newest dependency releases.
+* **`canary.yml` (new):** scheduled Mondays 12:00 UTC plus manual dispatch. Installs the latest release of every dependency and runs the full unit suite, so a breaking dependency release is caught within a week even with no pushes. The Giovanni test authenticates via the `EARTHDATA_USERNAME`/`EARTHDATA_PASSWORD` repository secrets.
+
+### Technical changes
+
+* **`tests/`** (new): nine test modules plus a `conftest.py` fixture that writes small georeferenced ESRI ASCII grids for the raster-based tests.
+* **`pytRIBS/classes.py`:** fixed the broken `input_file` path in the `Soil`, `Land`, `Met`, and `Mesh` constructors described above; each now also carries a populated `options` dictionary, consistent with `Results`.
+* **`pyproject.toml`:** added a `test` extra (`pytest>=8.0`) and pytest configuration — `testpaths`, the `network` marker, and `addopts = "-m 'not network'"` so the default `pytest` invocation never touches the network (an explicit `-m network` overrides it).
+* **`.github/workflows/install_and_test.yml`:** unit test step added.
+* **`.github/workflows/canary.yml`:** new scheduled canary workflow.
+
+### Running the tests
+
+```bash
+pip install -e ".[test]"
+
+pytest              # fast unit suite (network canaries excluded)
+pytest -m network   # live-service canaries (Giovanni needs Earthdata credentials)
+```
+

pyproject.toml

@@ -41,5 +41,16 @@ dependencies = [
 version = "1.0.0"
 requires-python = ">=3.11"
 
+[project.optional-dependencies]
+test = ["pytest>=8.0"]
+
 [project.urls]
-Repository = "https://github.com/tRIBS-Model/pytRIBS"
+Repository = "https://github.com/tRIBS-Model/pytRIBS"
+
+[tool.pytest.ini_options]
+testpaths = ["tests"]
+# Network canaries are excluded by default; run them with `pytest -m network`
+addopts = "-m 'not network'"
+markers = [
+    "network: tests that hit live external services (ISRIC, SOLUS, POLARIS, NASA Giovanni)",
+]

pytRIBS/classes.py

@@ -253,10 +253,11 @@ def __init__(self, input_file=None,meta=None):
         if meta is not None:
             self.meta=meta
 
+        # read_input_file assigns values into self.options in place
+        self.options = Infile.create_input_file()
         if input_file is not None:
-            options = Shared.read_input_file(input_file)
-        else:
-            options = Infile.create_input_file()
+            Shared.read_input_file(self, input_file)
+        options = self.options
 
         # Initialize attributes
         self.soilmapname = options['soilmapname']
@@ -305,10 +306,11 @@ def __init__(self, input_file=None,meta=None):
         if meta is not None:
             self.meta=meta
 
+        # read_input_file assigns values into self.options in place
+        self.options = Infile.create_input_file()
         if input_file is not None:
-            options = Shared.read_input_file(input_file)
-        else:
-            options = Infile.create_input_file()
+            Shared.read_input_file(self, input_file)
+        options = self.options
 
         # Initialize attributes
         self.landmapname = options['landmapname']
@@ -393,10 +395,11 @@ def __init__(self, preprocess_args=None, generate_mesh_args=None,
         if generate_mesh_args is not None:
             self.mesh_generator = GenerateMesh(*generate_mesh_args)
 
+        # read_input_file assigns values into self.options in place
+        self.options = Infile.create_input_file()
         if input_file is not None:
-            options = Shared.read_input_file(input_file)
-        else:
-            options = Infile.create_input_file()
+            Shared.read_input_file(self, input_file)
+        options = self.options
 
         # Initialize attributes
         self.pointfilename = options['pointfilename']
@@ -453,10 +456,11 @@ def __init__(self, input_file=None, meta=None):
         if meta is not None:
             self.meta = meta
 
+        # read_input_file assigns values into self.options in place
+        self.options = Infile.create_input_file()
         if input_file is not None:
-            options = Shared.read_input_file(input_file)
-        else:
-            options = Infile.create_input_file()
+            Shared.read_input_file(self, input_file)
+        options = self.options
 
         self.hydrometstations = options['hydrometstations']
         self.gaugestations = options['gaugestations']

tests/conftest.py

@@ -0,0 +1,33 @@
+import numpy as np
+import pytest
+import rasterio
+from rasterio.transform import from_origin
+
+from pytRIBS.shared.inout import InOut
+
+
+@pytest.fixture
+def make_ascii_raster(tmp_path):
+    """Factory that writes a small ESRI ASCII grid and returns its path.
+
+    Grids are placed in UTM zone 12N around a synthetic Arizona-like origin so
+    that any code touching the CRS or transform sees realistic values.
+    """
+
+    def _make(name, data, nodata=-9999.0, cellsize=100.0):
+        data = np.asarray(data, dtype='float32')
+        profile = {
+            'driver': 'AAIGrid',
+            'dtype': 'float32',
+            'count': 1,
+            'nodata': nodata,
+            'width': data.shape[1],
+            'height': data.shape[0],
+            'crs': rasterio.crs.CRS.from_epsg(32612),
+            'transform': from_origin(400000.0, 3900000.0, cellsize, cellsize),
+        }
+        path = str(tmp_path / name)
+        InOut.write_ascii({'data': data, 'profile': profile}, path)
+        return path
+
+    return _make

tests/test_ascii_io.py

@@ -0,0 +1,30 @@
+"""Round-trip tests for the ESRI ASCII raster reader/writer in InOut."""
+import numpy as np
+import pytest
+
+from pytRIBS.shared.inout import InOut
+
+
+def test_ascii_round_trip(make_ascii_raster):
+    data = [[1.5, 2.25], [3.75, 4.0]]
+    path = make_ascii_raster('grid.asc', data, cellsize=100.0)
+
+    raster = InOut.read_ascii(path)
+
+    assert raster['data'].tolist() == data
+    assert raster['profile']['nodata'] == -9999.0
+    assert raster['profile']['width'] == 2
+    assert raster['profile']['height'] == 2
+    # Cell size survives the header rewrite in write_ascii
+    assert raster['profile']['transform'].a == pytest.approx(100.0)
+
+
+def test_ascii_write_replaces_nan_with_nodata(make_ascii_raster):
+    data = [[1.0, np.nan], [np.nan, 4.0]]
+    path = make_ascii_raster('grid_nan.asc', data)
+
+    raster = InOut.read_ascii(path)
+
+    assert raster['data'][0, 1] == pytest.approx(-9999.0)
+    assert raster['data'][1, 0] == pytest.approx(-9999.0)
+    assert raster['data'][0, 0] == pytest.approx(1.0)

tests/test_classes.py

@@ -0,0 +1,59 @@
+"""Construction and input-file round-trip tests for the top-level pytRIBS classes."""
+import pytest
+
+from pytRIBS.classes import Land, Met, Mesh, Model, Soil
+
+
+@pytest.fixture
+def input_file(tmp_path):
+    """A tRIBS .in file written by Model with a handful of options set."""
+    m = Model()
+    m.options['startdate']['value'] = '06/01/2015/00/00'
+    m.options['runtime']['value'] = '240'
+    m.options['outfilename']['value'] = 'results/test'
+    m.options['soiltablename']['value'] = 'data/model/soil/soils.sdt'
+    m.options['soilmapname']['value'] = 'data/model/soil/soil_classes.soi'
+    m.options['landtablename']['value'] = 'data/model/land/land.ldt'
+    m.options['hydrometstations']['value'] = 'data/model/met/meteor/met.sdf'
+    m.options['gaugestations']['value'] = 'data/model/met/precip/precip.sdf'
+    m.options['pointfilename']['value'] = 'data/model/mesh/mesh.points'
+    path = str(tmp_path / 'test.in')
+    m.write_input_file(path)
+    return path
+
+
+def test_model_input_file_round_trip(input_file):
+    m = Model(input_file=input_file)
+    assert m.options['startdate']['value'] == '06/01/2015/00/00'
+    assert m.options['runtime']['value'] == '240'
+    assert m.options['outfilename']['value'] == 'results/test'
+
+
+def test_default_construction():
+    # Every top-level class must construct without an input file
+    assert Soil().soilmapname['value'] is None
+    assert Land().landtablename['value'] is None
+    assert Met().hydrometstations['value'] is None
+    assert Mesh().pointfilename['value'] is None
+
+
+def test_soil_from_input_file(input_file):
+    s = Soil(input_file=input_file)
+    assert s.soiltablename['value'] == 'data/model/soil/soils.sdt'
+    assert s.soilmapname['value'] == 'data/model/soil/soil_classes.soi'
+
+
+def test_land_from_input_file(input_file):
+    land = Land(input_file=input_file)
+    assert land.landtablename['value'] == 'data/model/land/land.ldt'
+
+
+def test_met_from_input_file(input_file):
+    met = Met(input_file=input_file)
+    assert met.hydrometstations['value'] == 'data/model/met/meteor/met.sdf'
+    assert met.gaugestations['value'] == 'data/model/met/precip/precip.sdf'
+
+
+def test_mesh_from_input_file(input_file):
+    mesh = Mesh(input_file=input_file)
+    assert mesh.pointfilename['value'] == 'data/model/mesh/mesh.points'

tests/test_evaluate_metrics.py

@@ -0,0 +1,26 @@
+"""Known-answer tests for the model evaluation metrics."""
+import numpy as np
+import pytest
+
+from pytRIBS.results.evaluate import Evaluate
+
+OBSERVED = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
+
+
+def test_perfect_simulation_scores():
+    simulated = OBSERVED.copy()
+    assert Evaluate.nash_sutcliffe(OBSERVED, simulated) == pytest.approx(1.0)
+    assert Evaluate.kling_gupta_efficiency(OBSERVED, simulated) == pytest.approx(1.0)
+    assert Evaluate.root_mean_squared_error(OBSERVED, simulated) == pytest.approx(0.0)
+    assert Evaluate.percent_bias(OBSERVED, simulated) == pytest.approx(0.0)
+
+
+def test_constant_offset_scores():
+    simulated = OBSERVED + 1.0  # uniform overestimation by 1
+    # NSE = 1 - 5/10
+    assert Evaluate.nash_sutcliffe(OBSERVED, simulated) == pytest.approx(0.5)
+    assert Evaluate.root_mean_squared_error(OBSERVED, simulated) == pytest.approx(1.0)
+    # Overestimation is negative PBIAS under this sign convention
+    assert Evaluate.percent_bias(OBSERVED, simulated) == pytest.approx(-100.0 / 3.0)
+    # r = 1, alpha = 1, beta = 4/3 -> KGE = 1 - 1/3
+    assert Evaluate.kling_gupta_efficiency(OBSERVED, simulated) == pytest.approx(2.0 / 3.0)

tests/test_ks_decay.py

@@ -0,0 +1,45 @@
+"""Tests for compute_ks_decay: recovery of the Ivanov et al. (2004) decay parameter f."""
+import numpy as np
+import pytest
+
+from pytRIBS.shared.inout import InOut
+from pytRIBS.soil.soil import SoilProcessor
+
+NODATA = -9999.0
+DEPTHS_MM = [1.0, 100.0, 400.0, 800.0]
+
+
+def ivanov_ks(k0, f, z):
+    """Ivanov et al. (2004) eqn 17: mean Ks over depth z given surface K0."""
+    z = np.asarray(z, dtype=float)
+    return k0 * (f * z) / (np.exp(f * z) - 1.0)
+
+
+def test_compute_ks_decay_recovers_known_f(make_ascii_raster, tmp_path):
+    k0 = 10.0  # mm/hr at the surface
+    f_slow, f_fast = 0.005, 0.02  # 1/mm
+
+    # Pixel layout:  (0,0) slow decay | (0,1) uniform profile (hits the floor)
+    #                (1,0) nodata     | (1,1) fast decay
+    grid_input = []
+    for depth in DEPTHS_MM:
+        data = np.array([
+            [ivanov_ks(k0, f_slow, depth), k0],
+            [k0, ivanov_ks(k0, f_fast, depth)],
+        ])
+        if depth == 100.0:
+            data[1, 0] = NODATA  # poison one layer of pixel (1,0)
+        path = make_ascii_raster(f'ks_{int(depth)}mm.asc', data, nodata=NODATA)
+        grid_input.append({'depth': depth, 'path': path})
+
+    out = str(tmp_path / 'f.asc')
+    SoilProcessor().compute_ks_decay(grid_input, output=out)
+
+    f_grid = InOut.read_ascii(out)['data']
+
+    assert f_grid[0, 0] == pytest.approx(f_slow, rel=0.15)
+    assert f_grid[1, 1] == pytest.approx(f_fast, rel=0.15)
+    # Uniform Ks profile must be floored at min_f, not left at ~0
+    assert f_grid[0, 1] == pytest.approx(1e-4, abs=1e-5)
+    # Nodata in any input layer propagates nodata to the output
+    assert f_grid[1, 0] == pytest.approx(NODATA)