feat: Add SOLUS soil dataset download and bedrock datasets

[JoshCederstrom]

This pull request adds the USDA SOLUS100 dataset as a third gridded soil data source and fixes the ROSETTA-based soil parameter workflow. The two changes are related: the parameter fixes apply to every source processed through ROSETTA, and SOLUS brings a depth-to-bedrock product that enables a new spatially variable initial groundwater tool for spin-up simulations.

The ROSETTA fixes address two stacked problems. First, rosetta-soil 0.3 (released 02/2026) changed rosetta() to return alpha, n, and Ksat in linear units instead of log10, so the back transform silently corrupted every derived grid in any environment installed after that release. The workflow now requests geometric-mean estimates (estimate_type="geo", matching the log10-mean convention of ROSETTA and of rosetta-soil < 0.3) 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 is not the ideal van Genuchten to Brooks-Corey conversion and overestimates the magnitude. It is now converted following Morel-Seytoux et al. (1996, WRR 32(5)), which preserves the effective capillary drive so that infiltration behavior is insensitive to the choice of retention model. As an independent check, the converted psib for a loam agrees within ~6% of the Rawls & Brakensiek (1982) literature value for that texture class. Pixels are also now predicted in large batched calls rather than one rosetta() call per pixel, which substantially speeds up the soil workflow on large watersheds.

Functions

The SOLUS tools live on the Soil class and mirror the structure of the existing ISRIC/POLARIS paths.

• run_soil_workflow(watershed, output_dir, source='SOLUS'): the existing workflow gains 'SOLUS' as a source option. Downloads the property grids, fills nodata, estimates the tRIBS parameter grids with ROSETTA (model code 3: sand, silt, clay, bulk density), computes the conductivity decay parameter f, and writes the soil classification map, table, and scgrid.gdf exactly as for the other sources.
• get_solus_grids(bbox, depths, variables): downloads SOLUS100 (100 m) soil property rasters via Cloud Optimized GeoTIFF windowed reads, only the pixels covering the domain are fetched over HTTP, never the CONUS-wide file and reprojects them to the project CRS. Currently mapped for the top three depth layers (0-5, 5-15, 15-30 cm).
• get_solus_bedrock(bbox): downloads the SOLUS100 depth-to-bedrock raster the same way, converts cm to m, writes it as an ASCII grid, and sets the BEDROCKFILE option.
• generate_initial_groundwater(bedrock_path, fraction, filename=None): builds a spatially variable initial water table by scaling a bedrock depth raster by a user-specified fraction (e.g. 0.7 places the initial water table 70% of the way from the surface to bedrock). Intended to simplify setting up the initial depth to bedrock raster, pairs naturally with get_solus_bedrock().
• process_solus_soil(grid_input, output=None, ks_only=False): converts downloaded SOLUS property grids into the tRIBS parameter grids (Ks, theta_r, theta_s, psib, m) through the same shared ROSETTA prediction and van Genuchten to Brooks-Corey conversion path used by the ISRIC workflow.

Usage

from pytRIBS.classes import Soil

soil = Soil()

# Full soil workflow from SOLUS100
soil.run_soil_workflow(watershed, 'data/model/soil', source='SOLUS')

# Depth-to-bedrock and a spatially variable initial water table for spin-up runs
bbox = list(watershed.bounds)  # [minx, miny, maxx, maxy] in the project CRS
bedrock = soil.get_solus_bedrock(bbox)
soil.generate_initial_groundwater(bedrock, fraction=0.7,
                                  filename='data/model/soil/initial_gw.asc')

Technical Changes

• pytRIBS/soil/soil.py
  • Added _rosetta_to_tribs_params(), the shared prediction core used by both the ISRIC and SOLUS paths so the ROSETTA handling and vG-to-BC conversions cannot drift apart. It batches pixels through rosetta() in chunks (rosetta materializes a 1000-member bootstrap ensemble per call, so chunking bounds memory on large watersheds), requests estimate_type="geo", converts psib via Morel-Seytoux et al. (1996) eq. 17.
  • Refactored process_raw_soil() (ISRIC) onto the shared core, removing the per-pixel loop and the now incorrect 10** back-transforms.
  • Added get_solus_grids(), get_solus_bedrock(), process_solus_soil(), and generate_initial_groundwater(). SOLUS bulk density is delivered scaled by 100 and is converted to g/cm3 before prediction — without this, ROSETTA silently rejects the out-of-range values and falls back to its texture-only model.
  • Added _fill_solus_nonsoil(), which repairs SOLUS non-soil cells after the standard nodata fill and writes the nonsoil_mask.asc footprint described above.
  • create_soil_map() accepts the SOLUS texture grid types (sandtotal/claytotal, already in %, no unit conversion).
  • compute_ks_decay() gained a min_f floor (default 1e-4) applied after curve fitting, preventing near-zero decay values that tRIBS cannot handle numerically where Ks is nearly uniform with depth.
  • run_soil_workflow() wires in the SOLUS branches (download, processing, soil map) and builds the decay-fit inputs with file-existence checks; parameter rasters are written with controlled decimal precision (grids at 2 decimals, groundwater at 0, f at 5), substantially reducing file sizes.
• pytRIBS/shared/inout.py: write_ascii() now detects ASCII grid header lines by content instead of assuming a fixed 6-line header, fixing a crash in the decimal-formatting path when GDAL writes dx/dy headers in place of cellsize.
• pyproject.toml: rosetta-soil>=0.1.2 raised to rosetta-soil>=0.3.