Every hydrological model starts from a map of precipitation and temperature that nobody measured. We measured points, at rain gauges and thermometers, and we filled the space between them with an interpolation. How well we did that filling sets the quality of everything downstream: snow, evapotranspiration, discharge, the water budget itself.
This post collects what I have learned, from the literature and from building a 1 km daily dataset for the Po River District (Salehi et al., preprint), about doing this filling well. It is written as a set of practical guidelines, with the papers behind each one.
1. A fine grid is not fine
information
The starting point is still Daly (2006): do not equate resolution with realism. A 1 km grid is a rendering choice. The information it contains is set by the stations and by the spatial structure of the field, not by the grid spacing.
Daly also tells us where the difficulty lies. Terrain and coastal effects matter least at scales above 100 km and most below 10 km. With a typical station spacing of 100 km they cannot be resolved, “except in densely populated regions of developed countries”. Northern Italy is one of those regions, so the question is not whether to interpolate below 10 km, but how to say where the result is informative and where it is not.
Two things are often read into Daly’s paper that it does not say:
- it does not reject interpolation below 10 km; its conclusions depend on the setting (flat land, mountains without rain shadows, complex terrain) and on station density;
- it does not reject cross-validation; it warns against comparing cross-validation errors computed on different networks, which is a different and correct point.
2. Know your network before choosing a method
Most interpolation problems are network problems. Before any method, describe the network:
- Stations active each day, not in total. A network of 1,500 stations with 50% missing data can mean 500 stations active in the 1990s and 1,100 after 2000. In the Po dataset, the active precipitation network grew from 472 stations (1991–2000) to 1,108 (2011–2020).
- Distances, not densities. For each cell and day, compute the distance to the nearest station and to the k-th nearest, where k is the number of stations your interpolation uses. These maps show where the field is supported and where it is extrapolated. In the Po dataset, the median distance to the seventh-nearest station fell from about 28 km to 15 km over the three decades.
- Hypsometry. Compare the elevation distribution of stations with that of the DEM, as Manara et al. (2026) and Pavan et al. (2026) do. High elevations are almost always under-sampled.
- Time of observation. Many manual gauges in Italy recorded 9-to-9 local time; ARCIS corrects for it (Pavan et al., 2019). Mixing 9-to-9 and 00–24 days shifts daily values and misaligns them with discharge.
- Record length and homogeneity. Short series are useful for spatial detail. Long homogenized series are needed for trends. You rarely get both from the same selection.
3. Temperature and precipitation need different treatment
An elevation trend works at the daily scale for temperature, and only at the climatological scale for precipitation. Every dataset I reviewed is consistent with this.
Temperature. The lapse rate is physically constrained and changes smoothly from day to day, so a daily regression on elevation, with kriging or weighting of the residuals, works well. The main difficulty is inversions, frequent in the Po Valley and in Alpine valleys in winter. A single lapse rate for a large domain misses them. Pavan et al. (2026) show that fitting a piecewise (two-slope) lapse rate in each of 28 sub-areas gives lower cross-validation errors, every year, than one regression for the whole domain. In our own dataset the lapse rate is domain-wide, re-estimated daily, and temperature errors grow from 1.4 °C below 500 m to 2.1 °C above 1,500 m. That is the price.
Precipitation. On a single day the relation with elevation is weak and changes with the weather. Kirshbaum et al. (2018) explain why:
- in unblocked flow, convection starts on windward slopes;
- in blocked flow, it starts upstream or downstream of the barrier;
- in summer, thermally driven flows make crests hot spots of convective initiation.
The same mountain can be wetter or drier than the plain depending on the day. When we tried a daily elevation trend for precipitation, the slopes were unstable and, multiplied by the elevation of 1 km cells, produced an obviously artificial topographic imprint. We dropped it.
At the climatological scale the relation is robust, mostly because it rains more often at higher elevations. This is why the best precipitation products use elevation only for the normals:
- the anomaly method: interpolate monthly normals with an elevation model, then interpolate the anomalies (or daily fractions) separately. APGD (Isotta et al., 2014), Crespi et al. (2021) and IH-GAR (Manara et al., 2026) work this way;
- PRISM-like local regression for the normals (Daly et al., 2008; Crespi et al., 2018), with weights that account for distance, elevation, slope, aspect and distance from the sea.
The anomaly method has a further advantage that a reviewer of our paper rightly stressed: anomalies are spatially smoother than absolute values, so they are much less sensitive to missing stations.
It also has a cost that is easy to miss when the goal is a daily field. A monthly normal is the number of wet days times their mean intensity. Scaling a station’s daily value by the ratio between the cell’s normal and the station’s normal multiplies each storm by the ratio of wet-day frequencies as well as of intensities. If a ridge is wetter than the valley because it rains there more often, the anomaly method gives the ridge the valley’s wet days with heavier storms. Monthly totals come out right by construction, so monthly validation cannot see the problem; wet-day frequency and intensity can. For monthly climatology this does not matter. For daily forcing of a hydrological model, where the intensity of each storm controls runoff and the split between rain and snow, it does. There is a further reason: the elevation gradient of a monthly normal is an average over the weather types of that month, each with its own gradient (Kirshbaum et al., 2018), so it describes no single storm.
The remedy is to ask the data where elevation acts. Group wet stations into events, and estimate the elevation effect on the probability of rain and on storm intensity from contrasts between stations in the same storm (a mixed model with one random effect per event does this). Then put elevation where it is found: if it acts on occurrence, interpolate occurrence with an elevation drift and intensity without, and combine them. To claim that intensity does not depend on elevation, use an equivalence test against a margin fixed in advance; a non-significant slope is not enough. One practical detail: the product of an occurrence probability and an intensity is an expected value, which smears the wet/dry boundary again. For a daily map, decide which cells are wet by thresholding the probability, or produce conditional simulations.
4. The main families of methods
| Method | Weights from | Elevation | Uncertainty per cell | Typical use |
|---|---|---|---|---|
| Inverse distance, Shepard | Prescribed distance functions | Only if added as a trend | No | ARCIS precipitation (Pavan et al., 2019) |
| Angular distance weighting (SPHEREMAP) | Prescribed, with angular correction | No | Yes, with the Yamamoto (2000) method | EMO-5 and EMO-1 (Thiemig et al., 2022) |
| Local weighted regression (PRISM-like) | Prescribed kernels on distance and terrain similarity | Explicit | No | Precipitation normals (Crespi et al., 2018) |
| Anomaly method | Regression for normals, distance weighting for anomalies | Through the normals | No | APGD, IH-GAR |
| Ordinary kriging | Variogram fitted to the data | No | Kriging variance | Precipitation, when no drift is reliable |
| Detrended kriging or kriging with external drift | Variogram of residuals | Through the trend | Kriging variance | Daily temperature |
| Occurrence–intensity (hurdle) kriging | Variograms of indicator and of log-intensity | In occurrence, intensity or both, as the data show | Kriging variance of each part | Daily precipitation |
Two remarks on kriging, since I am a long-time user of it in GEOframe (Bancheri et al., 2018):
- It is often lumped together with inverse distance weighting, as in Daly’s tables. That throws away its most useful part: the variogram, which is a description of the data. In kriging the spatial structure is estimated from the observations, not assumed through prescribed weights or tuned kernels; the estimate is the best linear unbiased one for that structure; and stations close to each other share their weight instead of doubling it, which no distance-weighting scheme does.
- Its variance per cell is a real advantage, but not a complete one. For a given variogram it depends only on station geometry, not on how different the neighbouring values are. Thiemig et al. (2022) chose the Yamamoto (2000) estimator for this reason. Report the kriging variance together with cross-validation errors, not instead of them.
5. Read the variogram
The variogram is not just a step in the kriging machinery; it is information about the field. Its three parameters have a direct reading:
- Nugget: variability at scales below the station spacing, plus measurement error. A large nugget says the network does not resolve the field; denser sampling would be needed.
- Sill: the total variance of the field. A sill that keeps growing with distance signals a trend that should be modelled.
- Range: the distance over which values stay correlated. If it is longer than the typical distance to the stations used, the network resolves the dominant structure that day. If it is shorter, kriging returns towards the mean, and its variance says so.
Fitting a variogram every day is possible, and in our case it is done automatically with particle swarm optimization over five models. Two cautions:
- automatic daily fitting can be unstable; Thiemig et al. (2022) excluded kriging from EMO-5 partly for this reason. Keep fit diagnostics and flag bad days;
- a variogram fitted over a whole domain assumes the same structure everywhere. In large and varied domains, regional variograms may be needed.
6. Evaluate the errors honestly
Leave-one-out cross-validation (LOO) is the most useful error estimate we have: remove a station, estimate it from the others, repeat for all stations. Used carefully, it is sound. Its limits are well known (Daly, 2006):
- It sees only places with stations. Terrain without stations, often the highest, produces no error estimate at all.
- It is optimistic with clustered stations, because a twin station is always there to predict the one removed.
- Do not use it to tune and then to evaluate. If bandwidths or smoothing parameters are chosen by minimizing LOO error, the same LOO error is optimistic, and the selection favours smooth fields. Removing stations that the method predicts poorly lowers the reported error further.
- Compare methods only on the same network. Station-by-station comparison, with the same neighbours, is a valid paired test. MAE values from different networks, periods or wet-day thresholds are not comparable: in Daly’s example, removing the hard high-elevation stations lowered the MAE from 49.7 to 17.2 mm.
Complements that make an error assessment convincing:
- Stratify LOO by station support and elevation. In the Po dataset, below 500 m, precipitation RMSE rises from 5.4 mm where the seventh-nearest station is within 10 km to 13.9 mm where it is 50–75 km away.
- Check that the kriging variance is calibrated: standardized LOO errors, (z − ẑ)/σ, should have mean near 0 and variance near 1.
- Withhold a block of stations, for example all those above 1,500 m, to test extrapolation (Daly’s Table V).
- Thin the network in a data-rich period, deleting half the stations at random many times, to measure sensitivity to missing data.
- Read residuals against predicted values, not observed ones. Kriging always flattens peaks, so residuals against observations always tilt and prove nothing.
- Validate at the time step you distribute, and score wet-day frequency and intensity separately. Monthly or annual statistics let errors of opposite sign cancel; a daily product needs daily validation.
- Remember what no gauge comparison can show. Snow undercatch is shared by all gauges, so the reference carries the same bias as the estimate. Only independent data, such as discharge, snow or radar, can reveal it, and each is confounded with something else.
7. Report what the grid contains
The last guideline is about communication. Daly’s advice to dataset developers is to state strengths and limits frankly. In practice:
- say “provided on a 1 km grid”, not “1 km resolution”;
- release daily station counts and support maps (distance to the nearest stations);
- release the per-cell uncertainty and the station-wise cross-validation errors, as separate products;
- release the daily variogram parameters and lapse rates;
- release station metadata and, where licences allow, the quality-controlled series;
- state the known biases: undercatch, high-elevation underestimation, valleys where the field is smoother than reality.
A user who knows, for every cell and every day, how much local information the value contains can decide whether it is good enough. That is more useful than any single number for “effective resolution”.
A short checklist
References
- Bancheri, M., Serafin, F., Bottazzi, M., Abera, W., Formetta, G., Rigon, R. (2018). The design, deployment, and testing of kriging models in GEOframe with SIK-0.9.8. Geosci. Model Dev., 11, 2189–2207. doi:10.5194/gmd-11-2189-2018
- Crespi, A., Brunetti, M., Lentini, G., Maugeri, M. (2018). 1961–1990 high-resolution monthly precipitation climatologies for Italy. Int. J. Climatol., 38, 878–895. doi:10.1002/joc.5217
- Crespi, A., Matiu, M., Bertoldi, G., Petitta, M., Zebisch, M. (2021). A high-resolution gridded dataset of daily temperature and precipitation records (1980–2018) for Trentino-South Tyrol. Earth Syst. Sci. Data, 13, 2801–2818. doi:10.5194/essd-13-2801-2021
- Daly, C. (2006). Guidelines for assessing the suitability of spatial climate data sets. Int. J. Climatol., 26, 707–721. doi:10.1002/joc.1322
- Daly, C., et al. (2008). Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous United States. Int. J. Climatol., 28, 2031–2064. doi:10.1002/joc.1688
- Isotta, F. A., et al. (2014). The climate of daily precipitation in the Alps: development and analysis of a high-resolution grid dataset from pan-Alpine rain-gauge data. Int. J. Climatol., 34, 1657–1675. doi:10.1002/joc.3794
- Kirshbaum, D. J., Adler, B., Kalthoff, N., Barthlott, C., Serafin, S. (2018). Moist orographic convection: physical mechanisms and links to surface-exchange processes. Atmosphere, 9, 80. doi:10.3390/atmos9030080
- Manara, V., et al. (2026). A new daily high-resolution gridded precipitation dataset for the Italian Greater Alpine Region (1951–2023). J. Hydrol.: Reg. Stud., 67, 103776. doi:10.1016/j.ejrh.2026.103776
- Pavan, V., et al. (2019). High resolution climate precipitation analysis for north-central Italy, 1961–2015. Clim. Dyn., 52, 3435–3453. doi:10.1007/s00382-018-4337-6
- Pavan, V., et al. (2026). A new operational dataset of gridded minimum and maximum temperature over north-central Italy 1991 to present. Climate Services, 43, 100689. doi:10.1016/j.cliser.2026.100689
- Salehi, H., et al. (2026). A 30-year 1-km daily precipitation and air temperature dataset for the Po River District (Italy). Earth Syst. Sci. Data Discuss. doi:10.5194/essd-2026-467
- Thiemig, V., et al. (2022). EMO-5: a high-resolution multi-variable gridded meteorological dataset for Europe. Earth Syst. Sci. Data, 14, 3249–3272. doi:10.5194/essd-14-3249-2022
- Yamamoto, J. K. (2000). An alternative measure of the reliability of ordinary kriging estimates. Math. Geol., 32, 489–509. doi:10.1023/A:1007577916868
