full_documentation.md

Climate Data Processing Pipeline (data_flow)

Table of Contents

1. Purpose
2. Architecture
   • Repository structure
   • Entry point and execution flow
   • Two operating modes
   • R5 reference classes
   • Configuration
   • Supported variables
   • Outputs
3. Installation and Setup
4. Examples of Use
5. Pipeline Stage Details
   • Pre-processing (pp)
   • Gap-filling (gf)
   • Gap-filling validation (gf_valid)
   • Homogenisation (hg)
   • Climatologies (cl)
   • Gridding (gr)
   • Gridding validation (gr_valid)
   • Grid comparison (gr_comp)
6. R5 Class Field Reference
7. Lines of Improvement

Purpose

This repository contains an R-based pipeline for automated post-processing of quality-controlled meteorological station data from Spain. It is the second stage of the PTI+ Clima climate data workflow, running immediately after the quality_control repository. The pipeline takes the cleaned station observations produced by that first stage and transforms them into a set of gap-filled, homogenised daily time series and regular spatial grids — everything needed to underpin climate services, climate indices, and gridded climatological products.

More concretely, the pipeline applies six sequential processing steps to each meteorological variable:

1. Pre-processing (pp) — Crops observations to the analysis period, derives composite variables (e.g. thermal range from max and min temperature, relative humidity from dew-point), classifies stations as candidates or auxiliaries, and computes inter-station distance and correlation matrices.
2. Gap-filling (gf) — Estimates missing values in each candidate station series by regressing on neighbouring auxiliary series; supports additive-delta, multiplicative-ratio, log-ratio, direct, and quantile-mapping methods.
3. Gap-filling validation (gf_valid) — Leave-one-out cross-validation of the gap-filling step, producing per-station accuracy diagnostics.
4. Homogenisation (hg) — Detects and corrects systematic shifts in each candidate station series using the Standard Normal Homogeneity Test (SNHT; Alexandersson 1986), applied recursively at both daily and monthly time scales.
5. Climatologies (cl) — Computes per-station daily climatological statistics (mean, standard deviation, 5th / 95th percentiles, absolute extremes) using a 31-day centred moving window.
6. Gridding (gr) — Creates daily regular grids for mainland Spain plus the Balearic Islands (UTM zone 30N) and the Canary Islands (UTM zone 28N) separately, using universal kriging with station elevation and distance to the coast as co-variates.

Optional validation and comparison steps (gr_valid, gr_comp) can also be run to assess spatial-interpolation accuracy and to compare the produced grids against reference products.

The pipeline is designed to be run unattended, from the command line or as a scheduled job, without any user interaction once it has been configured. All per-variable settings live in a single YAML file, so processing a new variable or rerunning with modified parameters requires only an edit to that file. Outputs are daily NetCDF grids and per-stage .rds R objects that are fully self-contained and reproducible.

Relationship to quality_control: The QC pipeline must be run first. Either its .RData output files (local mode) or its database tables (database mode) provide the raw input for this pipeline. Neither repository can substitute for the other; together they form the complete climate station data processing chain.

Figure: The full PTI+ Clima data workflow. The quality_control step (left) produces cleaned station data, which data_flow then processes through pre-processing, gap-filling, homogenisation, climatologies, and gridding to produce daily gridded climate products.

Architecture

Repository structure

data_flow/
├── R/
│   ├── main.R               ← Entry point; orchestrates the full pipeline
│   ├── functions.R          ← All functions and R5 class definitions
│   ├── config.yml           ← Per-variable configuration options
│   ├── pp_report.Rmd        ← Pre-processing summary report template
│   ├── gf_report.Rmd        ← Gap-filling summary report template
│   ├── gf_validation.Rmd    ← Gap-filling validation report template
│   ├── hg_report.Rmd        ← Homogenisation summary report template
│   ├── gr_report.Rmd        ← Gridding summary report template
│   ├── gr_validation.Rmd    ← Gridding validation report template
│   └── gr_comparison.Rmd    ← Grid comparison report template
├── data_raw/
│   ├── grd_pen.rda          ← Support grid: mainland Spain + Balearics (UTM 30N)
│   └── grd_can.rda          ← Support grid: Canary Islands (UTM 28N)
└── data/                    ← Output directory (created at runtime)
    ├── pp/
    ├── gf/
    ├── hg/
    ├── cl/
    └── gr/

Entry point and execution flow

R/main.R is the sole entry point. It is driven by argparser for command-line argument parsing and config::get() for reading per-variable settings from R/config.yml. All pipeline logic is in R/functions.R, which defines both standalone helper functions and five R5 reference-class objects (pp, gf, hg, cl, gr). The main script sources this file, then executes conditional blocks for each requested stage — each block reads the .rds output of the previous stage, processes data, writes a new .rds, and renders an HTML report.

Two operating modes

Database mode (default): Station data are read from and written to a PostgreSQL database (aemet on dana-sc-database). Quality-controlled input data live in the quality_control schema; intermediate results are written to the data_flow schema. This is the intended production mode and requires network access to the database server.

Local mode (flag -l / --local): Data are read from and written to local .rds and .RData files. The QC output files produced by the quality_control repository must be placed in the directory specified by config$dir$qc_dir. This mode is useful for development, for running without database access, or for reproducing results from archived files.

R5 reference classes

Each pipeline stage is represented by an R5 reference-class object. These objects carry all data, metadata, configuration, and processing methods in a single self-contained unit. Because R5 objects are modified in place (using <<- assignment), methods do not need to return values; they update the object's fields directly. The five classes are:

• pp — pre-processed data; fields: sts (station metadata), dat (data matrix), dates, ovr / dst / dis (overlap, correlation-distance, Euclidean-distance matrices), config, args, info.
• gf — gap-filled data; inherits the pp fields and adds coef (regression coefficients) and stats.
• hg — homogenised data; adds fill (binary matrix flagging gap-filled values), brk_pts (detected break points), hg_scores, hg_stats. The initialize method accepts either a gf object or individual field arguments, which allows the class to be constructed from archived data without rerunning prior stages.
• cl — climatologies; adds cli (list of six 365 × N matrices: mean, sd, q05, q95, max, min).
• gr — gridded data; methods drive kriging and NetCDF writing across both spatial domains.

Configuration

All per-variable settings are in R/config.yml. There is a default: section with settings shared across all variables; individual variable sections (e.g. tmax:, pr:) selectively override these defaults. The main fields are:

Field	Description
dir$qc_dir	Path to quality-control output (input to pp)
dir$pp_dir … dir$gr_dir	Output directories for each stage
period_analysis$start/end	Analysis period (1961-01-01 to 2025-10-31)
var$name	Internal variable identifier (e.g. "tmax")
var$AEMET_name	Raw AEMET variable code(s) in the QC data
var$scale	"relative" (temperature-like) or "absolute" (precipitation-like)
var$unit	Physical unit string
var$factor	Scale factor applied to raw AEMET values
var$range	Physical plausibility range (lower, upper)
cand$min_years$active/inactive	Minimum data years for candidate classification
cand$exclude_auto	Logical: exclude automatic stations from candidates
aux$min_years	Minimum data years for auxiliary status
aux$distance	Distance metric for auxiliary selection: "correlation" or "euclidean"
aux$max_dist_km	Maximum distance (km) for auxiliary stations
gf$method	Gap-filling method: "diff", "ratio", "logratio", "direct", or "qmap"
gf$seasonal	Logical: use seasonally varying gap-filling coefficients
hg$correction	Homogenisation correction type: "daily", "monthly", "annual", or "none"
hg$method	Correction method: "diff" or "ratio"
hg$clevel_daily/monthly	SNHT confidence levels (%) for daily and monthly tests
hg$n_neighbours	Number of reference neighbours for SNHT

Supported variables

Config key	Description	AEMET source variable(s)	Unit
tmax	Maximum daily temperature	TMAX	°C
tmin	Minimum daily temperature	TMIN	°C
trange	Daily thermal range	TMAX, TMIN	°C
pr	Daily total precipitation	P	mm day⁻¹
hr	Mean daily relative humidity (4 synoptic obs.)	HU00, HU07, HU13, HU18	%
ssrd	Daily total global radiation	RGLODIA, TOTSOL	kJ m⁻² day⁻¹
ws	Mean daily wind speed (4 synoptic obs.)	VEL_00, VEL_07, VEL_13, VEL_18	km h⁻¹

Additional variables can be added by creating new sections in config.yml following the same structure as the existing entries.

Outputs

Stage	Output	Description
pp	<pp_dir>/<var>.rds	Pre-processed R5 object
pp	<pp_dir>/<var>.html	Summary HTML report
gf	<gf_dir>/<var>.rds	Gap-filled R5 object
gf	<gf_dir>/<var>.html	Summary HTML report
gf	<gf_dir>/<var>/<stn>.html	Per-station diagnostic reports
gf_valid	<gf_dir>/<var>_valid.html	Cross-validation report
hg	<hg_dir>/<var>.rds	Homogenised R5 object
hg	<hg_dir>/<var>.html	Summary HTML report
cl	<cl_dir>/<var>.rds	Climatology R5 object
gr	<gr_dir>/<var>_pen.nc	Daily grids, mainland + Balearics (NetCDF)
gr	<gr_dir>/<var>_can.nc	Daily grids, Canary Islands (NetCDF)
gr	<gr_dir>/<var>.html	Gridding summary report
gr_valid	<gr_dir>/<var>_valid.html	Gridding cross-validation report

The .rds objects are fully self-contained: they embed all data, metadata, configuration, and R5 methods, and can be loaded independently of the pipeline to reproduce any result.

Installation and Setup

Prerequisites

• R ≥ 4.2
• R packages (installed automatically via pacman when the pipeline first runs): tidyverse, argparser, config, chron, Rfast, snowfall, future, terra, sf, rnaturalearth, gstat, ncdf4, abind, stringi, RPostgres, dplyr, automap, FNN
• Pandoc — required for rendering R Markdown reports; ships with RStudio. Set the path in main.R if not using RStudio:

  Sys.setenv(RSTUDIO_PANDOC = "/Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools")

• For database mode: network access to dana-sc-database and read/write permissions on the aemet PostgreSQL database.
• For local mode: the .RData output files produced by a prior run of the quality_control pipeline for the same variable.

Installation steps

1. Clone the repository:

   git clone https://github.com/PTI-Clima/data_flow.git
   cd data_flow

2. Open data_flow.Rproj in RStudio, or set the working directory manually:

   setwd("path/to/data_flow")

3. (Optional) Verify R package dependencies by running:

   pacman::p_load(tidyverse, argparser, config, chron, Rfast, snowfall,
                  future, terra, sf, rnaturalearth, gstat, ncdf4, abind,
                  stringi, RPostgres, dplyr, automap, FNN)

   Any missing packages will be installed automatically on the first pipeline run.

4. Review R/config.yml:

   • In local mode: update dir$qc_dir to the folder containing the .RData QC output files, and update the dir$pp_dir, dir$gf_dir, dir$hg_dir, dir$cl_dir, and dir$gr_dir paths to appropriate writable locations.
   • In database mode: database credentials are read from a local ~/.pgpass file or environment variables; update default$db$host and default$db$port in config.yml if necessary.
   • Optionally adjust period_analysis$start / end if you are working with a different time window.

5. Verify the setup with a short trial run (see the next section).

Examples of Use

Basic command-line usage

Run the full pipeline for maximum temperature using the database:

Rscript R/main.R tmax

Run only gap-filling and homogenisation in local mode:

Rscript R/main.R pr -l -p gf hg

Run a quick trial with the first 100 stations and no reports (useful for verifying the setup):

Rscript R/main.R tmax -t -l -n --no_indiv_rep

Get the full help message:

Rscript R/main.R --help

The available arguments are:

Argument	Description
var	Variable to process (must match a section in config.yml)
-t / --trial	Trial mode: process only the first 100 stations
-l / --local	Local mode: read/write files instead of database
-n / --no_global_rep	Skip the global HTML report
--no_indiv_rep	Skip per-station HTML reports
-p / --procs	Stages to run, space-separated (default: all in order)
-v / --verbosity	Verbosity level: 0 = silent, >0 = status messages (default: 1)
-c / --ncores	CPU cores: 0 = all available, default 1

Processing all variables in sequence

For a full operational update of all configured variables, a shell script like the following can be run as a cron job or batch task:

#!/bin/bash
for var in tmax tmin trange pr hr ssrd ws; do
    Rscript R/main.R $var -l
done

On a Unix-like system, to run this every night at 03:00, add a cron entry:

0 3 * * * cd /path/to/data_flow && bash run_all.sh >> logs/run.log 2>&1

Running individual pipeline stages

To rerun only the gridding stage for precipitation after updating station data, without recomputing the earlier stages:

Rscript R/main.R pr -l -p gr

To run gap-filling and its validation together:

Rscript R/main.R tmax -l -p gf gf_valid

Interactive use in R

To run or debug interactively, you can source main.R inside an R session after manually setting argument values. The read_inline_args() function at the top of main.R accepts an arg_values list for this purpose:

# Simulate a command-line call from within an R session
source("R/functions.R")
args <- read_inline_args(list(
  var       = "tmax",
  local     = TRUE,
  procs     = c("pp", "gf"),
  verbosity = 1,
  ncores    = 1
))

You can then step through the pipeline blocks in main.R interactively, inspect intermediate R5 objects (PP, GF, etc.) in the environment, and debug individual methods.

Inspecting outputs

After a successful run, the .rds objects can be loaded directly in R:

PP <- readRDS("data/pp/tmax.rds")

# Station metadata
head(PP$sts)

# Data matrix: rows = days, columns = stations
dim(PP$dat)

# Check which stations are candidates
table(PP$sts$cand, useNA = "always")

For homogenised data:

HG <- readRDS("data/hg/tmax.rds")

# How many break points were detected?
length(HG$brk_pts)

# Summary statistics for the homogenisation corrections
head(HG$hg_stats)

NetCDF grids can be inspected with the terra or ncdf4 packages:

library(terra)
r <- rast("data/gr/tmax_pen.nc")
plot(r[[1]])  # map of the first day

Pipeline Stage Details

This section describes the algorithms, criteria, and diagnostics for each stage of the pipeline in detail. All configuration parameters referenced here are defined in R/config.yml; the Architecture section above contains the complete parameter table.

Pre-processing (pp)

The pre-processing step reads quality-controlled station data from the database or from local .RData files, crops observations to the configured analysis period, derives composite variables where needed, classifies each station as candidate, auxiliary, or excluded, and computes inter-station similarity matrices.

Derived variables

Some variables are not directly stored in the AEMET database but are computed during pre-processing:

• trange: daily thermal range, computed as tmax − tmin.
• hr: mean daily relative humidity, computed as the average of observations at 00:00, 07:00, 13:00, and 18:00 UTC.
• ssrd: daily total global radiation, selecting from RGLODIA and TOTSOL as available per station.
• ws: mean daily wind speed from observations at 00:00, 07:00, 13:00, and 18:00 UTC.

Station classification

Each station is classified into one of three categories:

• Candidate (cand == TRUE): stations with sufficient data to be gap-filled, homogenised, and used in gridding. Both active and inactive stations can be candidates if they meet the minimum year thresholds.
• Auxiliary (cand == FALSE): stations with enough data to serve as reference series for gap-filling and homogenisation, but not themselves targeted for correction.
• Excluded (cand == NA): stations with too little data, known quality problems, or other disqualifying characteristics.

The following filters are applied in sequence:

1. Minimum auxiliary years: stations with fewer than aux$min_years complete years of data are set to excluded.
2. Automatic station exclusion: if cand$exclude_auto = TRUE, all automatic weather stations are excluded.
3. Blacklist: station IDs listed in aux$exclude are excluded unconditionally.
4. Precipitation data quality: for pr, stations with coarse-resolution data (0.5 mm or 1 mm step) are excluded, as these indicate old rain gauges that cannot resolve small rainfall amounts.
5. Candidate threshold: stations that pass all exclusion filters are promoted to candidate if they have ≥ cand$min_years$active years (for active stations) or ≥ cand$min_years$inactive years (for inactive stations).
6. Whitelist: station IDs in aux$include are forced to candidate status regardless of other criteria.
7. Province / island coverage: the year threshold is relaxed iteratively (by one year at a time, up to a maximum of cand$min_years$inactive relaxations) for any spatial domain (province or island) that would otherwise have fewer than cand$min_stations_per_prov candidates. This ensures at least minimal spatial coverage in data-sparse regions.
8. Proximity deduplication: among candidate pairs closer than 1 km (a common situation where both a manual and an automatic sensor exist at the same observatory), one is demoted to auxiliary using priority rules: (i) prefer the active station; (ii) if both are active, prefer the automatic one; (iii) if both are automatic and active, retain the longer series.

Inter-station matrices

The step computes two pairwise matrices over all non-excluded stations:

• ovr: fractional overlap of non-NA data between every pair of stations, used to assess how many simultaneous observations are available for computing regression coefficients.
• dst: correlation-based distance (1 − |r|), computed from overlapping observations.
• dis: Euclidean distance in kilometres.

Pre-processing report

The HTML report (pp_report.Rmd) includes:

• Full configuration dump
• Station classification maps (candidates in red, auxiliaries in blue, excluded in grey) for both spatial domains
• Data resolution table showing the distribution of measurement precision by station
• Tabular station lists (candidate, auxiliary, excluded) with metadata
• Daily time series plots for all candidate stations and a data availability stacked-bar chart by year

Gap-filling (gf)

Gap-filling estimates missing values in each candidate station series using neighbouring auxiliary stations as predictors. Auxiliary stations are ranked by the aux$distance metric ("correlation" or "euclidean"), and the top-ranked ones that overlap with the gap period are used.

Gap-filling methods

The method is set by gf$method:

Method	Formula	Best suited for
diff	aux + mean(cand − aux) over overlap	Temperature-type variables
ratio	aux × mean(cand / aux)	Precipitation, radiation
logratio	exp(log(aux) + mean(log(cand / aux)))	Skewed positive variables
direct	Linear regression of cand on aux	General use
qmap	Quantile mapping from aux distribution to cand distribution	Distribution-preserving fill

Seasonal coefficients

When gf$seasonal = TRUE, the regression offset or ratio is computed separately for each calendar month rather than over the full overlapping period. This allows the filling to track seasonal relationships (e.g. using April-only temperature differences to fill April gaps), which is particularly important for variables with strong seasonal cycles.

Outlier truncation

After filling, values outside the domain-specific plausibility range are clipped. The pipeline uses a two-level guard: the configured physical range (var$range) is the absolute limit, while values beyond the historical domain-level extremes by more than var$outlier_tolerance percent are also flagged and clipped.

Gap-filling report

The HTML report (gf_report.Rmd) includes:

• Histograms and spatial maps of the percentage of completeness before and after gap-filling
• Outlier detection: counts of values above and below the domain historical extremes, separately for the peninsular and Canary Islands domains
• A station table listing the number of remaining missing days after filling
• Daily time series plots for all candidate stations, with links to per-station individual HTML reports that show the observed and filled series side by side

Gap-filling validation (gf_valid)

The validation step applies leave-one-out (LOO) cross-validation to assess gap-filling accuracy. For each candidate station in turn, its observed values are treated as missing and re-estimated from the remaining stations using the same algorithm. The resulting predicted/observed pairs are then scored with a comprehensive set of statistics.

S-mode and T-mode analyses

Statistics are computed in two complementary modes:

• T-mode (temporal, per station): all time steps for a single station are aggregated, revealing which stations are harder to fill.
• S-mode (spatial, per day): all stations for a single day are aggregated, revealing whether performance degrades during spatially coherent extreme events.

Both modes are broken down by calendar month to expose seasonal variation in performance.

Magnitude statistics

Statistic	Description	Ideal value
MAE	Mean absolute error	0
ME	Mean error (bias)	0
r	Pearson's correlation	1
rSD	Ratio of predicted to observed standard deviation	1
KGE	Kling-Gupta efficiency (blend of r, ME, rSD)	1
d	Willmott's index of agreement	1

Binary classification statistics (for pr and ws only — event/no-event at zero threshold)

Statistic	Description
PPV	Positive predictive value (precision): fraction of predicted-wet days that were observed wet
NPV	Negative predictive value: fraction of predicted-dry days that were observed dry
TPR	True positive rate (sensitivity / recall): fraction of observed-wet days correctly predicted
TNR	True negative rate (specificity): fraction of observed-dry days correctly predicted
F1	Harmonic mean of PPV and TPR
MCC	Matthews correlation coefficient: overall binary classification accuracy

Report contents

The validation report (gf_validation.Rmd) includes: predicted vs. observed scatter plots (hexbin density); Q-Q plots by month; split violin plots of observed and predicted distributions; tabular summary of all statistics by month; S-mode time series of 5th and 95th quantile bias over the analysis period; and spatial maps of per-station statistics.

Homogenisation (hg)

The homogenisation step detects and corrects systematic level shifts in each candidate station series. It uses the Standard Normal Homogeneity Test (SNHT; Alexandersson, 1986) applied recursively, first at the daily scale and optionally at the monthly scale.

Reference series construction

For each candidate station, a reference series Q is constructed from its hg$n_neighbours nearest neighbours (ranked by the same distance metric as gap-filling). Each neighbour is standardised to z-scores and the neighbours are averaged. The candidate is then expressed as a ratio or additive difference relative to this consensus, producing a series that captures only station-specific anomalies — genuine local climate signals cancel out in the reference.

SNHT algorithm (daily pass)

1. The SNHT statistic T(k) is computed at every position k in the test series Q. The position k* where T(k) is maximised is the candidate break point.
2. If T(k*) exceeds the critical value at the configured confidence level (hg$clevel_daily, one of 90, 92, 94, 95, 97.5, or 99 %), the break is accepted. Critical values follow Khaliq and Ouarda (2007), with an optional multiplicative correction factor.
3. A correction is computed for the segment before k* using the configured hg$method ("diff" for additive corrections, "ratio" for multiplicative). The granularity is set by hg$correction:
   • "daily": a smooth seasonal curve is fitted, allowing the correction to vary through the year.
   • "monthly": separate factors for each of the 12 calendar months.
   • "annual": a single factor for the full break period.
   • "none": detection only, no correction applied.
4. The correction is applied and the test is repeated on the corrected series. This recursive loop continues until no further significant break points are found.

Monthly pass

If hg$correction_monthly is not "none", an additional SNHT pass is performed on monthly aggregates (using hg$clevel_monthly) after the daily corrections are complete. This catches lower-frequency seasonal biases that may not be detectable at daily resolution.

Post-homogenisation station rejection

After all corrections, the method hg_update_cand() compares each candidate's long-term trend to the median trend of its five nearest neighbours. Stations whose trend difference exceeds the hg$rejection_threshold$annual and/or hg$rejection_threshold$monthly percentile cutoffs (computed from the distribution of absolute trend differences across all candidates) are demoted to auxiliary. This removes stations with residual inhomogeneities that the SNHT correction did not fully resolve.

Homogenisation report

The HTML report (hg_report.Rmd) includes:

• Configuration summary with correction settings in plain language
• Post-correction outlier check: counts of values above/below domain historical extremes
• Interactive summary DataTable (sortable, searchable, downloadable as CSV/Excel) with one row per candidate station, showing: retention status (retained/rejected, colour-coded green/red); number of break points at daily (N (d)) and monthly (N (m)) scales; maximum SNHT test statistic; annual and monthly trend differences, colour-coded by magnitude using a Spectral palette; and for precipitation, the wet-day proportion trend difference
• Per-row expandable diagnostic plots: the full corrected time series with break-point markers, a correction heatmap (magnitude by year × month), and where monthly correction was applied, a 4×3 grid of monthly trend panels showing candidate vs. regional trend
• Spatial maps of retained and rejected stations

Climatologies (cl)

(This step is ignored in the current status of the workflow, and climatologies and calculated in the index_calc routines.)

Gridding (gr)

The gridding step interpolates the homogenised daily station values onto a regular grid using universal kriging. Two spatial domains are processed separately using dedicated support grids:

• Peninsular domain (pen): mainland Spain, the Balearic Islands, and the Spanish North African territories. Grid is in UTM zone 30N, loaded from data_raw/grd_pen.rda.
• Canary Islands domain (can): UTM zone 28N, loaded from data_raw/grd_can.rda.

Universal kriging with external drift

For each day and each domain, universal kriging is applied with two external drift (co-variate) variables:

• Elevation (alt, m): accounts for the temperature lapse rate and orographic effects on precipitation and radiation.
• Distance to coast (dis, km): captures maritime / continental gradients in temperature, humidity, and wind.

The empirical variogram is fitted automatically using the automap package (automap::autofitVariogram()), which selects the best-fitting variogram model (spherical, exponential, Matérn, etc.) from a candidate set. Kriging is then performed on the fitted variogram to produce a complete grid for each day.

Precipitation-specific treatment

For pr, a two-step approach is used:

1. Indicator kriging: first determines the wet/dry classification at each grid point, using a binary (0/1) indicator variable.
2. Ordinary kriging on log-transformed wet-day amounts: interpolates rainfall amounts only for grid points predicted as wet.

Outputs

Grids are written to NetCDF files following CF conventions:

• <gr_dir>/<var>_pen.nc — full time series for the peninsular domain
• <gr_dir>/<var>_can.nc — full time series for the Canary Islands domain

Gridding report

The HTML report (gr_report.Rmd) shows example grid maps for the first, middle, and last day of the analysis period, with observed station values overlaid on the interpolated field, for both domains.

Gridding validation (gr_valid)

The gridding validation applies the same LOO cross-validation framework as the gap-filling validation, but at the kriging step. For each station in turn, its value is withheld from the interpolation and then predicted at that location using the kriging model fitted to the remaining stations.

The same set of statistics as gf_valid is computed in both T-mode and S-mode: MAE, ME, r, rSD, KGE, d, and (for precipitation and wind speed) PPV, NPV, TPR, TNR, F1, MCC.

The validation report additionally breaks down statistics separately for automatic and manual observatories (section 4 of gr_validation.Rmd), to assess whether the two station networks show systematic differences in interpolation accuracy. This is particularly relevant for precipitation, where automatic tipping-bucket gauges may have different systematic errors from manual observers.

Spatial maps of per-station statistics allow the identification of structurally difficult regions (mountain ranges, coastlines, island peripheries) where interpolation is consistently less accurate.

Grid comparison (gr_comp)

The grid comparison step is available for variables where reference gridded products exist. The produced LCSC grid is compared against up to six reference datasets:

Reference	Description
AEMET	Official AEMET gridded climate product
ECMWF	ERA5 reanalysis from ECMWF
E-OBS	Pan-European observation-based gridded dataset
ERA5 Land	High-resolution ERA5 land-surface reanalysis
HARMONIE	NWP-based high-resolution analysis
LCSC old	Previous version of the LCSC temperature grid

Differences are quantified in two modes:

• Spatial mode: statistics are computed grid-point by grid-point, averaged over the full time series, and displayed as spatial maps. Metrics: mean difference (bias), MAE, RMSE, rSD, and percentile differences at the 1st, 10th, 90th, and 99th percentiles.
• Temporal mode: a spatially-averaged metric is computed for each day and displayed as a time series with a fitted polynomial trend. This reveals whether the agreement between datasets drifts over time (e.g. due to changes in observation density or reference product versions).

Results for all metrics and all reference datasets are also summarised as violin plots, providing a compact comparison of the error distribution across products.

R5 Class Field Reference

Each pipeline stage is represented by an R5 reference-class object. The objects carry all data, metadata, and configuration, and can be loaded independently from their .rds files. Below is a complete field reference for each class.

pp — Pre-processed data

Field	Type	Description
sts	data.frame	Station metadata: id, name, lon, lat, alt, utm_x, utm_y, prov, sp_domain (province or island), sp_domain_big ("pen" or "can"), is_auto, is_active, cand (TRUE/FALSE/NA), complete_years, resol, source
dat	matrix	Data matrix: rows = days, columns = stations (all, not just candidates); values in native physical units
dates	Date	Vector of daily dates for the rows of dat
ovr	matrix	N × N pairwise fractional data overlap matrix
dst	matrix	N × N correlation-based distance matrix (1 −
dis	matrix	N × N Euclidean distance matrix (km)
config	list	Full configuration list for this variable, as read from config.yml
args	list	Command-line arguments used for this run
info	list	Run metadata: R version, platform, start time, hostname

gf — Gap-filled data

Inherits all pp fields, plus:

Field	Type	Description
coef	list	Per-station regression coefficients; structure depends on gf$seasonal: if TRUE, a list of 12 month-specific coefficient vectors; if FALSE, a single per-station vector
stats	data.frame	Per-station gap-filling summary: number of gaps filled, method used, fraction complete before and after

hg — Homogenised data

Inherits all gf fields, plus:

Field	Type	Description
fill	matrix	Logical matrix (same dimensions as dat) where TRUE marks gap-filled values; used to distinguish original from reconstructed observations during homogenisation
brk_pts	list	Per-candidate-station list of detected break-point dates (as Date vectors); NA if no break was detected
hg_scores	list	Per-station SNHT test statistics (full T(k) series) from the daily pass
hg_stats	data.frame	Summary table: one row per candidate station; columns include ID, Type, N_d (break points daily), Score_d (max SNHT statistic daily), N_m, Score_m, Trend_diff_d, Trend_diff_m, Wetday_diff (precipitation only), plus per-month trend diagnostics and embedded HTML diagnostic plots

cl — Climatologies

Inherits all hg fields, plus:

Field	Type	Description
cli	list	Named list of six 365 × N_cand numeric matrices, one per statistic: mean, sd, q05, q95, max, min. Row index = day of year (1–365); column index = candidate station

gr — Gridded data

The gr class does not inherit from the station-data classes. Its fields are:

Field	Type	Description
sts	data.frame	Station metadata (same structure as pp$sts)
dat	matrix	Data matrix used as input to kriging
dates	Date	Vector of daily dates
config	list	Configuration list
pts	list	Two named sf point objects ($pen, $can) holding the station coordinates and observed values per domain, used for kriging input
grd	list	Two named lists of daily grids ($pen, $can); each element is an sf or SpatRaster object with the interpolated field for one day
var	list	Two named numeric vectors ($pen, $can) of spatially-averaged daily values, derived from the grids
comp	list	Grid comparison results (populated by gr_comp): named sub-lists per reference product, each containing a spatial data frame $sp and a temporal data frame $ts with error metrics
info	list	Run metadata

Lines of Improvement

While the pipeline provides a complete and operational solution for automated climate data processing, several areas could be enhanced to improve usability, robustness, and maintainability.

Packaging and modularity. The pipeline is currently structured as a pair of scripts (main.R + functions.R) plus a YAML config. Wrapping it as a proper R package (with a DESCRIPTION, NAMESPACE, and roxygen2-generated documentation) would improve dependency management, allow installation via devtools::install_github(), and enable a pkgdown documentation site. Within functions.R, splitting the code into thematic files (one per pipeline stage) would make the codebase easier to navigate and maintain.

Unit and integration testing. There are currently no automated tests. Introducing a testthat test suite — even just for the core standalone functions such as fill_gaps, snht, error_stats, and krige_variograms — would guard against regressions and clarify expected function behaviour. An integration test that runs the full pipeline on a small synthetic dataset and checks that outputs match known values would be particularly valuable, as it would catch subtle interactions between stages that unit tests alone cannot detect.

Error handling and operational robustness. In production, the pipeline is expected to run unattended for many hours. Currently, if a single station or a single day fails inside a parallel loop, the entire stage may stop or produce a silent bad result. Adding structured error handling (e.g. tryCatch around per-station operations) and writing failures to a log file would make it much easier to diagnose problems in long runs. Similarly, pre-flight validation — checking that the QC output exists and is readable before starting — would give clearer early messages.

Parallelism improvements. The gap-filling coefficient computation and the kriging loop use snowfall for parallelism. Migrating to the future / future.apply framework (already a declared dependency) would give more flexible parallelism, better error propagation, and compatibility with cluster and cloud computing environments. The gridding stage in particular, which applies kriging independently to each of many thousands of days, is an embarrassingly parallel problem and would benefit from a more efficient parallel backend.

Database schema and versioning. In database mode, results are written to the data_flow schema without explicit versioning. Introducing a version or run-date column in the database tables, and keeping previous runs, would make it possible to compare results between pipeline runs (e.g. after changing homogenisation parameters). A lightweight provenance record — storing which config values and software versions were used for a given run — would also improve reproducibility.

Homogenisation tuning. The SNHT procedure currently uses fixed confidence levels (clevel_daily, clevel_monthly) configured per variable. An adaptive approach — for example, adjusting the confidence level based on the number of reference stations available — could reduce false-positive break-point detections in data-sparse periods. Additionally, the current implementation detects breaks sequentially (removing the most significant break and retesting) but does not yet implement a full multiple-changepoint model, which could improve accuracy for series with many closely spaced breaks.

Report quality and interactivity. The current R Markdown reports produce static HTML. Migrating to Quarto or adding plotly-based interactive plots (zooming, station selection, time series comparison) would make the diagnostic reports considerably more useful for operational review. A summary dashboard aggregating results across all variables and all stages in a single view would be particularly valuable for monitoring the health of the full pipeline.

Extensibility to new variables and spatial domains. Adding a new variable currently requires adding a new section in config.yml and, in some cases, a new aggregation branch in functions.R (e.g. a new formula in hr_tdew). Making the variable-specific logic more declarative — encoding it fully in the config where possible — would reduce the code changes needed to add a new variable. Similarly, extending the pipeline to additional spatial domains (e.g. the Spanish territories in Africa) would require adding new support grids to data_raw/ and adjusting the domain-routing logic in the gridding stage.