Spatial interpolations

Using the vectors of MeteoData and StationData as filled by the IOInterface::readMeteoData call as well as a grid of elevations (DEM, stored as a DEMObject), it is possible to get spatially interpolated parameters.

First, an interpolation method has to be selected for each variable which needs interpolation. Then the class computes the interpolation for each 2D grid point, combining the inputs provided by the available data sources. Any parameter of MeteoData can be interpolated, using the names given by meteoparam. One has to keep in mind that the interpolations are time-independent: each interpolation is done at a given time step and no memory of (eventual) previous time steps is kept. This means that all parameters and variables that are automatically calculated get recalculated anew for each time step.

- Note
- Please keep in mind that you need to specify a proper cartesian coordinate system in your [Input] section in order to be able to perform most of the spatial interpolations (Lat/lon coordinates are
**not**cartesian, they are spherical!).

Practically, the user has to specify in his configuration file (typically io.ini), for each parameter to be interpolated, which spatial interpolations algorithms should be considered, in the [Interpolations2D] section. This is provided as a space separated list of keywords (one per interpolation algorithm). Please notice that some algorithms may require extra arguments. Then, each algorithm will be evaluated (through the use of its rating method) and receive a grade (that might depend on the number of available data, the quality of the data, etc). The algorithm that receives the higher score within the user list, will be used for interpolating the selected variable at the given timestep. This means that at another timestep, the same parameter might get interpolated by a different algorithm. An example of such section is given below:

[Interpolations2D]

TA::algorithms = IDW_LAPSE AVG_LAPSE

TA::avg_lapse::rate = -0.008

RH::algorithms = LISTON_RH IDW_LAPSE AVG_LAPSE AVG

PSUM::algorithms = IDW_LAPSE AVG_LAPSE AVG CST

PSUM::avg_lapse::rate = 0.0005

PSUM::avg_lapse::frac = true

PSUM::cst::value = 0

VW::algorithms = IDW_LAPSE AVG_LAPSE

P::algorithms = STD_PRESS

P::std_press::USE_RESIDUALS = true

ILWR::algorithms = AVG_LAPSE

ILWR::avg_lapse::rate = -0.03125

RSWR::algorithms = IDW AVG

ISWR::algorithms = SWRAD

The keywords defining the algorithms are the following:

- NONE: returns a nodata filled grid (see NoneAlgorithm)
- STD_PRESS: standard atmospheric pressure as a function of the elevation of each cell (see StandardPressureAlgorithm)
- CST: constant value in each cell (see ConstAlgorithm)
- NEAREST: the value at the closest station is taken for each cell (see NearestNeighbourAlgorithm)
- AVG: average of the measurements in each cell (see AvgAlgorithm)
- AVG_LAPSE: constant value reprojected to the elevation of the cell (see AvgLapseRateAlgorithm)
- IDW: Inverse Distance Weighting averaging (see IDWAlgorithm)
- IDW_LAPSE: Inverse Distance Weighting averaging with reprojection to the elevation of the cell (see IDWLapseAlgorithm)
- IDW_SLOPES: IDW_LAPSE with separate processing for each of the 4 aspects+flat before merging with weighted average (see IDWSlopesAlgorithm)
- LIDW_LAPSE: IDW_LAPSE restricted to a local scale (n neighbor stations, see LocalIDWLapseAlgorithm)
- LISTON_RH: the dew point temperatures are interpolated using IDW_LAPSE, then reconverted locally to relative humidity (see RHListonAlgorithm)
- ILWR_EPS: the incoming long wave radiation is converted to emissivity and then interpolated (see ILWREpsAlgorithm)
- SWRAD: The atmospheric attenuation and splitting coefficients are evaluated and used to compute the short wave radiation with topographic shading (see SWRadInterpolation)
- LISTON_WIND: the wind field (VW and DW) is interpolated using IDW_LAPSE and then altered depending on the local curvature and slope (taken from the DEM, see ListonWindAlgorithm)
- RYAN: the wind direction is interpolated using IDW and then altered depending on the local slope (see RyanAlgorithm)
- WINSTRAL: the solid precipitation is redistributed by wind according to (Winstral, 2002) (see WinstralAlgorithm)
- PSUM_SNOW: precipitation interpolation according to (Magnusson, 2011) (see SnowPSUMInterpolation)
- PPHASE: precipitation phase parametrization performed at each cell (see PPHASEInterpolation)
- ODKRIG: ordinary kriging (see OrdinaryKrigingAlgorithm)
- ODKRIG_LAPSE: ordinary kriging with lapse rate (see LapseOrdinaryKrigingAlgorithm)
- USER: user provided grids to be read from disk (if available, see USERInterpolation)
- ALS_SCALING: scaling from Airborn Laser Scan data (see ALS_Interpolation)
- SNOWLINE: assimilation of snowline elevation information from external data sources (see SnowlineAlgorithm)

Several algorithms use elevation trends, all of them relying on the same principles: the lapse rates are recomputed at each time steps (see section Lapse rates), all stations' data are detrended with this lapse rate, the residuals are spatially interpolated with the algorithm as configured by the user and finally, the values at each cell are retrended (ie the lapse rates are re-applied using the cell's elevation).

The altitudinal trends are currently modelled as a linear relation. The slope of this linear relation can sometimes be provided by the end user (through his io.ini configuration file), otherwise it is computed from the data. In order to bring slightly more robustness, if the correlation between the input data and the computed linear regression is not good enought (below 0.7, as defined in Interpol2D::LinRegression), the same regression will get re-calculated with one point less (cycling throught all the points). The best result (ie: highest correlation coefficient) will be kept. If the final correlation coefficient is less than 0.7, a warning is displayed.

A list of supported options controlling the lapse rates is given in Trend::Trend().

From the users's point of view, all that has to be done is instantiate an IOManager object and call its IOManager::getMeteoData method with an elevation model and a grid. The spatial interpolation will then be done according to the settings in the INI file. The following lines show an example of the workflow:

Config cfg("io.ini");

IOManager io(cfg);

//reading the dem (necessary for several spatial interpolations algoritms)

DEMObject dem;

io.readDEM(dem);

//performing spatial interpolations

Grid2DObject param;

io.getMeteoData(date, dem, MeteoData::TA, param);

The interpolation algorithms have been inspired by the following papers:

*"A Meteorological Distribution System for High-Resolution Terrestrial Modeling (MicroMet)"*, Liston and Elder, Journal of Hydrometeorology**7**(2006), 217-234.*"Simulating wind ﬁelds and snow redistribution using terrain-based parameters to model snow accumulation and melt over a semi-arid mountain catchment"*, Adam Winstral and Danny Marks, Hydrological Processes**16**(2002), 3585– 3603. DOI: 10.1002/hyp.1238*"Quantitative evaluation of different hydrological modelling approaches in a partly glacierized Swiss watershed"*, Jan Magnusson, Daniel Farinotti, Tobias Jonas and Mathias Bavay, Hydrological Processes, 2010, under review.*"Modelling runoff from highly glacierized alpine catchments in a changing climate"*, Matthias Huss, Daniel Farinotti, Andreas Bauder and Martin Funk, Hydrological Processes,**22**, 3888-3902, 2008.*"Geostatistics for Natural Resources Evaluation"*, Pierre Goovaerts, Oxford University Press, Applied Geostatistics Series, 1997, 483 p., ISBN 0-19-511538-4*"Statistics for spatial data"*, Noel A. C. Cressie, John Wiley & Sons, revised edition, 1993, 900 p.