rainonsnow

rainonsnow is a research repository for exploring rain-on-snow hydrology in the Alps, focussing on extreme runoff as a consequence.

The repo consists of:

• A code pipeline for bringing in raster data from Google Earth Engine and analyzing it.
• A set of Shiny apps for exploring the data and models.
• An R package for the distributional learning model and helper functions.

Current Workflow

Scripts are numbered in pipeline order (scripts/1-* … scripts/7-*). Intermediate files land under derived/ (and yearly NetCDF under derived/eo/).

1. scripts/1-download_data-eo.py — Download hourly ERA5-Land from Google Earth Engine to derived/eo/era5_land_hourly_alps_<year>.nc (bands and years from inputs/data_specifications.yaml).
2. scripts/2-tablify_spatial_eo.r — Raster NetCDF → tabular hourly series; writes derived/era5_land_hourly_alps_all.rds.
3. scripts/3-pot_spatial_eo.r — Peaks over threshold per cell (derived/era5_land_hourly_alps_peaks.rds, thresholds, metadata).
4. scripts/4-distributional_learning.r — Quantile regression forest on POT peak hours (models and predictions used by later apps).
5. scripts/5-runoff_marginals.r — Marginal return-level tables for runoff (derived/era5_land_hourly_alps_dl_*).
6. scripts/6-drivers_joint_distribution.r — Joint rainfall–snowmelt model per cell.
7. scripts/7-likeliest_rain_snow.r — Conditional rain–snow structure given extreme runoff (optional follow-on).

Dependency management

Python dependencies are managed with uv (pyproject.toml; run uv sync so the lockfile matches the resolved environment). R package dependencies are declared in DESCRIPTION; install them from the repo root with something like devtools::install_deps() (or your preferred workflow). While editing package code, devtools::load_all() attaches the package from source.

Python Setup

This repository uses uv for the Python environment.

uv sync

The Python dependencies are declared in pyproject.toml.

Data specifications (inputs/data_specifications.yaml)

Pipeline settings for ERA5-Land export and tabular cleanup live in inputs/data_specifications.yaml:

• earth_engine.project_id — Google Cloud project passed to ee.Initialize(project=...).
• download — Earth Engine collection_id, first_year / last_year, and variables (hourly band names exported from the collection).
• tablify.epsilon_mm — Near-zero cutoff for rain and melt (mm/h equivalent) in script 2 before other steps.

Authentication is separate from this file: run once per machine/user:

earthengine authenticate

Run the downloader with:

uv run python scripts/1-download_data-eo.py

If data_specifications.yaml is missing, scripts 1 and 2 fall back to the same defaults as the template file.

Input Controls

YAML files under inputs/ configure scripts and apps (beyond data_specifications.yaml above):

• distributional_learning.yaml — Response and predictors for dl_rqforest in script 4.
• rain_snow_joint_model.yaml — Joint marginal + copula options for script 6 (key fit_joint_rain_snow_cells), and likeliest_rain_snow settings for script 7 (conditional rain–snow surface given runoff).
• pot_metadata.yaml — POT quantile / min-gap for script 3 and apps/3-pot-explorer.

R Setup

In addition to the analysis stream of this repository, this repository is also structured like an R package, with source files in R/, documentation in man/, and package metadata in DESCRIPTION and NAMESPACE. This is useful so that custom functions can be more easily accessed by the analysis scripts, and more easily used by anyone developing this analysis.

To make these functions available to the analysis scripts, run:

devtools::load_all()

If you make a change to the R package, run devtools::document() to update the package documentation.

Run The Analysis

Run the scripts in order (each writes inputs for the next) from the repo root.

uv python scripts/1-download_data-eo.py
Rscript scripts/2-tablify_spatial_eo.r
Rscript scripts/3-pot_spatial_eo.r
Rscript scripts/4-distributional_learning.r
Rscript scripts/5-runoff_marginals.r
Rscript scripts/6-drivers_joint_distribution.r

Shiny apps

Interactive apps live under apps/*, with folder names numbered to match the scripts they support (e.g. apps/4a-distributional-learning-fit and apps/4b-distributional-learning-diagnostics for script 4). Run them from the repository root so paths such as derived/ and devtools::load_all() resolve:

shiny::runApp("apps/<app-folder>")

Typical dependencies include shiny, ggplot2, tidyverse, sf, and rnaturalearth; individual apps may need yaml, probaverse, distionary, famish, rvinecopulib, and others used in the analysis scripts.

POT explorer (apps/3-pot-explorer)

Inspect peaks over threshold (POT) extraction: hourly runoff for one grid cell and year with the POT threshold and marked peaks; event timing scatter of all POT peaks for that cell (runoff vs day of year, all years pooled). Pick the cell from the map or sidebar; optional re-run of script 3 from the sidebar.

Data: derived/era5_land_hourly_alps_all.rds (scripts/2-tablify_spatial_eo.r) and derived/era5_land_hourly_alps_peaks.rds (scripts/3-pot_spatial_eo.r).

shiny::runApp("apps/3-pot-explorer")

Distributional learning — fit (apps/4a-distributional-learning-fit)

Tune dl_rqforest hyperparameters (hints under each control), fit one cell as an in-memory preview, then run script 4 for all cells (writes inputs/distributional_learning.yaml and derived outputs). Model: hourly runoff ~ rainfall + snowmelt on POT peaks (fixed).

Data: derived/era5_land_hourly_alps_peaks.rds (script 3).

shiny::runApp("apps/4a-distributional-learning-fit")

Distributional learning — diagnostics (apps/4b-distributional-learning-diagnostics)

Explore fitted models: P–P calibration, skill versus marginal, map, and conditional runoff CDFs at clicked rain–snow pairs. Fitted models load at startup (CDF works for every cell); load P–P / skill when you want calibration plots (from derived/era5_land_hourly_alps_dl_diagnostics.rds when present; otherwise predictions RDS or rebuild from models).

Data: peaks + derived/era5_land_hourly_alps_dl_rqforest_models.rds; derived/era5_land_hourly_alps_dl_diagnostics.rds after script 4 (optional derived/era5_land_hourly_alps_dl_predictions.rds for fallback recompute).

shiny::runApp("apps/4b-distributional-learning-diagnostics")

Return-level explorer (apps/5-return-level-explorer)

Map of marginal runoff return levels by cell, frequency–magnitude curves (forest mixture vs GP tail), and rain–snow likelihood surfaces at a chosen return period.

Data: peaks from script 3; script 4 models; precomputed marginal return levels from script 5 (derived/era5_land_hourly_alps_dl_marginal_return_levels.rds or the bundle described in the app header), or derived/era5_land_hourly_alps_dl_predictions.rds as a slower fallback.

shiny::runApp("apps/5-return-level-explorer")

Joint rainfall–snowmelt explorer (apps/6-joint-rain-snow-explorer)

Marginal fits and copula per cell: Gaussian-score diagnostics, joint density contours, marginal histograms, and frequency–magnitude curves for rainfall and snowmelt. Optional re-run of joint fitting from the sidebar (writes inputs/rain_snow_joint_model.yaml and runs scripts/6-drivers_joint_distribution.r).

Data: derived/era5_land_hourly_alps_all.rds and joint output from script 6 (derived/era5_land_hourly_alps_joint_rain_snow.rds); fitting options and script 7 settings share inputs/rain_snow_joint_model.yaml.

shiny::runApp("apps/6-joint-rain-snow-explorer")

Runoff marginals explorer (apps/5-runoff-marginals-explorer)

Side-by-side frequency–magnitude curves: distributional-learning marginals (Random Forest mixture vs GP conversion) and naive POT-only marginals (distionary::dst_empirical vs famish::fit_dst_gp on peak runoff). Map cell selection; optional matched y-axis limits and log return level on both panels.

Data: derived/era5_land_hourly_alps_peaks.rds and derived/era5_land_hourly_alps_dl_return_levels.rds from scripts/5-runoff_marginals.r.

shiny::runApp("apps/5-runoff-marginals-explorer")

R Package Demonstration

library(rainonsnow)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(probaverse)
#> ── Attaching core probaverse packages ──────────────────────────────────────────
#> ✔ distionary   0.1.0   Create and Evaluate Probability Distributions
#> ✔ distplyr     0.2.0   Manipulate and Combine Probability Distributions
#> ✔ famish       0.2.0   Flexibly Tune Families of Probability Distributions

The R package is titled rainonsnow and provides a small distributional learning interface for modelling the distribution of a target variable given some predictors.

For a simple workflow, fit a model with dl_rqforest() and then call predict(), using the mtcars dataset from the stats package.

df <- as_tibble(mtcars)
model <- dl_rqforest(
  data = df,
  yname = "hp",
  xnames = c("wt", "drat", "gear")
)

df <- mutate(df, distribution = predict(model), .before = everything())
df
#> # A tibble: 32 × 12
#>    distribution   mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear
#>    <list>       <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#>  1 <dst>         21       6  160    110  3.9   2.62  16.5     0     1     4
#>  2 <dst>         21       6  160    110  3.9   2.88  17.0     0     1     4
#>  3 <dst>         22.8     4  108     93  3.85  2.32  18.6     1     1     4
#>  4 <dst>         21.4     6  258    110  3.08  3.22  19.4     1     0     3
#>  5 <dst>         18.7     8  360    175  3.15  3.44  17.0     0     0     3
#>  6 <dst>         18.1     6  225    105  2.76  3.46  20.2     1     0     3
#>  7 <dst>         14.3     8  360    245  3.21  3.57  15.8     0     0     3
#>  8 <dst>         24.4     4  147.    62  3.69  3.19  20       1     0     4
#>  9 <dst>         22.8     4  141.    95  3.92  3.15  22.9     1     0     4
#> 10 <dst>         19.2     6  168.   123  3.92  3.44  18.3     1     0     4
#> # ℹ 22 more rows
#> # ℹ 1 more variable: carb <dbl>

The predictions are distributions. Take a look at the first distribution using the probaverse, for example, and plot its cdf.

plot(df$distribution[[1]], n = 1000)

The package also includes a null model to handle failures gracefully. For example, if you ask for na_action = "null" and the training data contain missing values, or if the training fails, dl_rqforest() returns a dl_null object instead of failing:

df2 <- as_tibble(mtcars)
df2$wt[2] <- NA_real_

dl_rqforest(
  data = df2,
  yname = "hp",
  xnames = c("wt", "drat", "gear")
)
#> <dl_rqforest>
#> response: hp
#> predictors: 3
#> training rows: 31

Predicting on these objects always returns a null distribution with the number of rows of the data:

predict(dl_null(), newdata = tibble(x = 1:2))
#> [[1]]
#> Null distribution (NA) 
#> 
#> [[2]]
#> Null distribution (NA)

Status

This repository is in an active research/prototyping state, so scripts, paths, and interfaces may change as the workflow evolves.

License

This repository is licensed under the MIT License - see the LICENSE file for details.