WSL/SLF GitLab Repository

InterpolationAlgorithms.h 30.3 KB
Newer Older
 Mathias Bavay committed Sep 23, 2013 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 /***********************************************************************************/ /* Copyright 2010 WSL Institute for Snow and Avalanche Research SLF-DAVOS */ /***********************************************************************************/ /* This file is part of MeteoIO. MeteoIO is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. MeteoIO is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with MeteoIO. If not, see . */ #ifndef __INTERPOLATIONALGORITHMS_H__ #define __INTERPOLATIONALGORITHMS_H__ #include #include #include #include #include #include #include #include namespace mio { class IOManager; class Meteo2DInterpolator; // forward declaration, cyclic header include /** * @page interpol2d 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 \ref 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. * * @section interpol2D_section Spatial interpolations section * 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: * @code * [Interpolations2D] * TA::algorithms = IDW_LAPSE CST_LAPSE * TA::cst_lapse = -0.008 * * RH::algorithms = RH IDW_LAPSE CST_LAPSE CST * * HNW::algorithms = HNW_SNOW IDW_LAPSE CST_LAPSE CST * HNW::hnw_snow = cst_lapse * HNW::cst_lapse = 0.0005 frac * * VW::algorithms = IDW_LAPSE CST_LAPSE * * P::algorithms = STD_PRESS * @endcode * * @section interpol2D_keywords Available algorithms * The keywords defining the algorithms are the following:  77  * - NONE: returns a nodata filled grid (see NoneAlgorithm)  Mathias Bavay committed Sep 23, 2013 78 79 80 81 82  * - STD_PRESS: standard atmospheric pressure as a function of the elevation of each cell (see StandardPressureAlgorithm) * - CST: constant value in each cell (see ConstAlgorithm) * - CST_LAPSE: constant value reprojected to the elevation of the cell (see ConstLapseRateAlgorithm) * - IDW: Inverse Distance Weighting averaging (see IDWAlgorithm) * - IDW_LAPSE: Inverse Distance Weighting averaging with reprojection to the elevation of the cell (see IDWLapseAlgorithm)  83  * - LIDW_LAPSE: IDW_LAPSE restricted to a local scale (n neighbor stations, see LocalIDWLapseAlgorithm)  Mathias Bavay committed Sep 23, 2013 84 85  * - RH: the dew point temperatures are interpolated using IDW_LAPSE, then reconverted locally to relative humidity (see RHAlgorithm) * - ILWR: the incoming long wave radiation is converted to emissivity and then interpolated (see ILWRAlgorithm)  Mathias Bavay committed Jun 14, 2014 86  * - 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)  Mathias Bavay committed May 21, 2014 87  * - RYAN: the wind direction is interpolated using IDW and then altered depending on the local slope (see RyanAlgorithm)  88  * - WINSTRAL: the solid precipitation is redistributed by wind according to (Winstral, 2002) (see WinstralAlgorithm)  Mathias Bavay committed Sep 23, 2013 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145  * - HNW_SNOW: precipitation interpolation according to (Magnusson, 2011) (see SnowHNWInterpolation) * - 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) * * @section interpol2D_lapse Lapse rates * Several algorithms use elevation trends, 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. * * @section interpol2D_dev_use Developer usage * From the developer's point of view, all that has to be done is instantiate an IOManager object and call its * IOManager::interpolate method. * @code * 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.interpolate(date, dem, MeteoData::TA, param); * * @endcode * * @section interpol2D_biblio Bibliography * 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 [NOT YET IMPLEMENTED] * - "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. * * @author Mathias Bavay * @date 2010-04-12 */ /** * @class InterpolationAlgorithm * @brief A class to perform 2D spatial interpolations. For more, see \ref interpol2d * * @ingroup stats * @author Thomas Egger * @date 2010-04-01 */ class InterpolationAlgorithm { public: InterpolationAlgorithm(Meteo2DInterpolator& i_mi, const std::vector& i_vecArgs, const std::string& i_algo, IOManager& iom) :  146  algo(i_algo), mi(i_mi), date(0., 0), vecArgs(i_vecArgs), vecMeteo(), vecData(),  Mathias Bavay committed Sep 23, 2013 147 148 149 150 151 152 153 154 155 156 157 158  vecMeta(), info(), param(MeteoData::firstparam), nrOfMeasurments(0), iomanager(iom) {}; virtual ~InterpolationAlgorithm() {}; //if anything is not ok (wrong parameter for this algo, insufficient data, etc) -> return zero virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param) = 0; virtual void calculate(const DEMObject& dem, Grid2DObject& grid) = 0; std::string getInfo() const; const std::string algo; protected: size_t getData(const Date& i_date, const MeteoData::Parameters& i_param, std::vector& o_vecData); size_t getData(const Date& i_date, const MeteoData::Parameters& i_param, std::vector& o_vecData, std::vector& o_vecMeta);  Mathias Bavay committed Feb 18, 2014 159  static size_t getStationAltitudes(const std::vector& i_vecMeta, std::vector& o_vecData);  Mathias Bavay committed Sep 23, 2013 160  void getTrend(const std::vector& vecAltitudes, const std::vector& vecDat, Fit1D &trend) const;  Mathias Bavay committed Feb 18, 2014 161 162  static void detrend(const Fit1D& trend, const std::vector& vecAltitudes, std::vector &vecDat, const double& min_alt=-1e4, const double& max_alt=1e4); static void retrend(const DEMObject& dem, const Fit1D& trend, Grid2DObject &grid, const double& min_alt=-1e4, const double& max_alt=1e4);  163  void simpleWindInterpolate(const DEMObject& dem, const std::vector& vecDataVW, const std::vector& vecDataDW, Grid2DObject &VW, Grid2DObject &DW);  Mathias Bavay committed Sep 23, 2013 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184  Meteo2DInterpolator& mi; Date date; const std::vector vecArgs; //we must keep our own copy, it is different for each algorithm! std::vector vecMeteo; std::vector vecData; /// vecMeta; ///& i_vecArgs, IOManager& iom); };  185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 /** * @class NoneAlgorithm * @brief Returns a nodata filled grid * This allows to tolerate missing data, which can be usefull if an alternate strategy could * later be used to generate the data (ie. a parametrization). This algorithm will only run * after all others failed. */ class NoneAlgorithm : public InterpolationAlgorithm { public: NoneAlgorithm(Meteo2DInterpolator& i_mi, const std::vector& i_vecArgs, const std::string& i_algo, IOManager& iom) : InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom) {} virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param); virtual void calculate(const DEMObject& dem, Grid2DObject& grid); };  Mathias Bavay committed Sep 23, 2013 202 203 204 205 /** * @class ConstAlgorithm * @brief Constant filling interpolation algorithm. * Fill the grid with the average of the inputs for this parameter.  206 207  * Optionally, it is also possible to provide the constant that should be used if no measurements * are avaiblable.  Mathias Bavay committed Sep 23, 2013 208 209 210 211 212 213  */ class ConstAlgorithm : public InterpolationAlgorithm { public: ConstAlgorithm(Meteo2DInterpolator& i_mi, const std::vector& i_vecArgs, const std::string& i_algo, IOManager& iom)  214  : InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom), user_cst(0.), user_provided(false) {}  Mathias Bavay committed Sep 23, 2013 215 216  virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param); virtual void calculate(const DEMObject& dem, Grid2DObject& grid);  217 218 219  private: double user_cst; bool user_provided;  Mathias Bavay committed Sep 23, 2013 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 }; /** * @class StandardPressureAlgorithm * @brief Standard atmospheric pressure interpolation algorithm. * Fill the grid with the standard atmosphere's pressure, depending on the local elevation. */ class StandardPressureAlgorithm : public InterpolationAlgorithm { public: StandardPressureAlgorithm(Meteo2DInterpolator& i_mi, const std::vector& i_vecArgs, const std::string& i_algo, IOManager& iom) : InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom) {} virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param); virtual void calculate(const DEMObject& dem, Grid2DObject& grid); }; /** * @class ConstLapseRateAlgorithm * @brief Constant filling with elevation lapse rate interpolation algorithm. * Assuming that average values occured at the average of the elevations, the grid is filled with average values * reprojected to real grid elevation according to a lapse rate. The lapse rate is either calculated from the data * (if no extra argument is provided), or given by the user-provided the optional argument "cst_lapse". * If followed by "soft", then an attempt to calculate the lapse rate from the data is made, any only if * unsuccessful, then user provided lapse rate is used as a fallback. If the optional user given lapse rate is * followed by "frac", then the lapse rate is understood as a fractional lapse rate, that is a relative change * of the value as a function of the elevation (for example, +0.05% per meters given as 0.0005). In this case, no attempt to calculate * the fractional lapse from the data is made. */ class ConstLapseRateAlgorithm : public InterpolationAlgorithm { public: ConstLapseRateAlgorithm(Meteo2DInterpolator& i_mi, const std::vector& i_vecArgs, const std::string& i_algo, IOManager& iom) : InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom) {} virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param); virtual void calculate(const DEMObject& dem, Grid2DObject& grid); }; /** * @class IDWAlgorithm * @brief Inverse Distance Weighting interpolation algorithm. * Each cell receives the weighted average of the whole data set with weights being 1/r² * (r being the distance of the current cell to the contributing station) and renormalized * (so that the sum of the weights is equal to 1.0). */ class IDWAlgorithm : public InterpolationAlgorithm { public: IDWAlgorithm(Meteo2DInterpolator& i_mi, const std::vector& i_vecArgs, const std::string& i_algo, IOManager& iom) : InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom) {} virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param); virtual void calculate(const DEMObject& dem, Grid2DObject& grid); }; /** * @class IDWLapseAlgorithm * @brief Inverse Distance Weighting interpolation algorithm with elevation detrending/reprojection. * The input data is projected to a reference elevation and spatially interpolated using an Inverse Distance * Weighting interpolation algorithm (see IDWAlgorithm). Then, each value is reprojected to the real * elevation of the relative cell. The lapse rate is either calculated from the data * (if no extra argument is provided), or given by the user-provided the optional argument "idw_lapse". * If followed by "soft", then an attempt to calculate the lapse rate from the data is made, any only if * unsuccessful or too bad (r^2<0.6), then the user provided lapse rate is used as a fallback. * If the optional user given lapse rate is * followed by "frac", then the lapse rate is understood as a fractional lapse rate, that is a relative change * of the value as a function of the elevation (for example, +0.05% per meters given as 0.0005). In this case, no attempt to calculate * the fractional lapse from the data is made. */ class IDWLapseAlgorithm : public InterpolationAlgorithm { public: IDWLapseAlgorithm(Meteo2DInterpolator& i_mi, const std::vector& i_vecArgs, const std::string& i_algo, IOManager& iom) : InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom) {} virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param); virtual void calculate(const DEMObject& dem, Grid2DObject& grid); }; /** * @class LocalIDWLapseAlgorithm * @brief Inverse Distance Weighting interpolation algorithm with elevation detrending/reprojection.  304  * The closest n stations (n being given as an extra argument of "lidw_lapse") to each pixel are  Mathias Bavay committed Sep 23, 2013 305  * used to compute the local lapse rate, allowing to project the contributions of these n stations to the  306 307  * local pixel with an inverse distance weight. Beware, this method sometimes produces very sharp transitions * as it spatially moves from one station's area of influence to another one!  Mathias Bavay committed Sep 23, 2013 308 309 310 311 312  */ class LocalIDWLapseAlgorithm : public InterpolationAlgorithm { public: LocalIDWLapseAlgorithm(Meteo2DInterpolator& i_mi, const std::vector& i_vecArgs,  Mathias Bavay committed Jan 02, 2014 313  const std::string& i_algo, IOManager& iom);  Mathias Bavay committed Sep 23, 2013 314 315  virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param); virtual void calculate(const DEMObject& dem, Grid2DObject& grid);  Mathias Bavay committed Jan 02, 2014 316 317  private: size_t nrOfNeighbors;  Mathias Bavay committed Sep 23, 2013 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 }; /** * @class RHAlgorithm * @brief Relative humidity interpolation algorithm. * This is an implementation of the method described in (Liston & Elder, 2006): for each input point, the dew * point temperature is calculated. Then, the dew point temperatures are spatially interpolated using IDWLapseAlgorithm. * Finally, each local dew point temperature is converted back to a local relative humidity. * * As a side effect, the user must have defined algorithms to be used for air temperature (since this is needed for dew * point to RH conversion) */ class RHAlgorithm : public InterpolationAlgorithm { public: RHAlgorithm(Meteo2DInterpolator& i_mi, const std::vector& i_vecArgs, const std::string& i_algo, IOManager& iom) : InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom), vecDataTA(), vecDataRH() {} virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param); virtual void calculate(const DEMObject& dem, Grid2DObject& grid); private: std::vector vecDataTA, vecDataRH; ///& i_vecArgs, const std::string& i_algo, IOManager& iom) : InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom), vecDataEA() {} virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param); virtual void calculate(const DEMObject& dem, Grid2DObject& grid); private: std::vector vecDataEA; ///"A meteorological distribution system for high-resolution terrestrial modeling (MicroMet)", Journal of Hydrometeorology, 7.2, 2006. * The wind speed and direction are spatially interpolated using IDWLapseAlgorithm. Then, the wind speed and * direction fields are altered by wind weighting factors and wind diverting factors (respectively) calculated * from the local curvature and slope (as taken from the DEM, see DEMObject). The wind diverting factor is * actually the same as in RyanAlgorithm.  Mathias Bavay committed Sep 23, 2013 372  */  Mathias Bavay committed Jun 14, 2014 373 class ListonWindAlgorithm : public InterpolationAlgorithm {  Mathias Bavay committed Sep 23, 2013 374  public:  Mathias Bavay committed Jun 14, 2014 375  ListonWindAlgorithm(Meteo2DInterpolator& i_mi,  Mathias Bavay committed Sep 23, 2013 376 377 378 379 380 381 382 383 384  const std::vector& i_vecArgs, const std::string& i_algo, IOManager& iom) : InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom), vecDataVW(), vecDataDW() {} virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param); virtual void calculate(const DEMObject& dem, Grid2DObject& grid); private: std::vector vecDataVW, vecDataDW; ///"a mathematical model for diagnosis and prediction of surface winds in mountainous terrain", * 1977, journal of applied meteorology, 16, 6. * The DEM is used to compute wind drection changes that are used to alter the wind direction fields. * @code * DW::algorithms = RYAN * @endcode */ class RyanAlgorithm : public InterpolationAlgorithm { public: RyanAlgorithm(Meteo2DInterpolator& i_mi, const std::vector& i_vecArgs, const std::string& i_algo, IOManager& iom)  401  : InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom), vecDataVW(), vecDataDW() {}  Mathias Bavay committed May 21, 2014 402 403  virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param); virtual void calculate(const DEMObject& dem, Grid2DObject& grid);  404 405  private: std::vector vecDataVW, vecDataDW; ///"Spatial snow modeling of wind-redistributed snow using terrain-based parameters", 2002, * Journal of Hydrometeorology, 3(5), 524-538. * The DEM is used to compute wind exposure factors that are used to alter the precipitation fields. * It is usually a good idea to provide a DEM that also contain the accumulated snow height in order  416  * to get a progressive softening of the terrain features.  417 418 419 420 421 422 423 424 425 426 427 428 429 430  * * This method must therefore first use another algorithm to generate an initial precipitation field, * and then modifies this field accordingly. This base method is "idw_lapse" by default. * Then it requires a synoptic wind direction that can be provided by different means: * - without any extra argument, the stations are located in the DEM and their wind shading (or exposure) * is computed. If at least one station is found that is not sheltered from the wind (in every direction), it * provides the synoptic wind (in case of multiple stations, the vector average is used). Please note that * the stations that are not included in the DEM are considered to be sheltered. If no such station * is found, the vector average of all the available stations is used. * - by providing a fixed synoptic wind bearing that is used for all time steps * - by providing the station_id of the station to get the wind direction from. In this case, the base algorithm * for generating the initial wind field must be specified in the first position. * * @remarks Only cells with an air temperature below freezing participate in the redistribution  Mathias Bavay committed May 07, 2014 431 432 433 434  * @code * HNW::algorithms = WINSTRAL * HNW::winstral = idw_lapse 180 * @endcode  Mathias Bavay committed Apr 11, 2014 435 436 437 438  */ class WinstralAlgorithm : public InterpolationAlgorithm { public: WinstralAlgorithm(Meteo2DInterpolator& i_mi,  439 440  const std::vector& i_vecArgs, const std::string& i_algo, IOManager& iom);  Mathias Bavay committed Apr 11, 2014 441 442 443  virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param); virtual void calculate(const DEMObject& dem, Grid2DObject& grid); private:  444 445  void initGrid(const DEMObject& dem, Grid2DObject& grid); static bool isExposed(const DEMObject& dem, Coords location);  446  static double getSynopticBearing(const std::vector& vecMeteo, const std::string& ref_station, const std::string& algo);  Mathias Bavay committed May 07, 2014 447  static double getSynopticBearing(const std::vector& vecMeteo);  448  static double getSynopticBearing(const DEMObject& dem, const std::vector& vecMeteo);  Mathias Bavay committed May 21, 2014 449   450 451  std::string base_algo, ref_station; double user_synoptic_bearing;  Mathias Bavay committed May 21, 2014 452  static const double dmax;  Mathias Bavay committed Apr 11, 2014 453 454 };  Mathias Bavay committed Sep 23, 2013 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 /** * @class USERInterpolation * @brief Reads user provided gridded data on the disk. * The grids are all in a directory that is given as the algorithm's argument. The files must be named * according to the following schema: * - {numeric date}_{capitalized meteo parameter}.asc, for example 200812011500_TA.asc * - Default_{capitalized meteo parameter}.asc for the grid to use when no measurements exist (which prevents * retrieving the date for the interpolation) * The meteo parameters can be found in \ref meteoparam "MeteoData". Example of use: * @code * TA::algorithms = USER * TA::user = ./meteo_grids * @endcode * */ class USERInterpolation : public InterpolationAlgorithm { public: USERInterpolation(Meteo2DInterpolator& i_mi, const std::vector& i_vecArgs, const std::string& i_algo, IOManager& iom) : InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom), filename() {nrOfMeasurments=0;} virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param); virtual void calculate(const DEMObject& dem, Grid2DObject& grid); private: std::string getGridFileName() const; std::string filename; }; /** * @class SnowHNWInterpolation * @brief Precipitation distribution according to the local slope and curvature. * The precipitation distribution is initialized using a specified algorithm (IDW_LAPSE by default, see IDWLapseAlgorithm). * An optional parameter can be given to specify which algorithm has to be used for initializing the grid. * Please do not forget to provide the arguments of the chosen algorithm itself if necessary! * * After this initialization, the pixels whose air temperatures are below or at freezing are modified according * to the method described in "Quantitative evaluation of different hydrological modelling approaches * in a partly glacierized Swiss watershed", Magnusson et Al., Hydrological Processes, 25, 2071-2084, 2011 and * "Modelling runoff from highly glacierized alpine catchments in a changing climate", Huss et All., Hydrological Processes, 22, 3888-3902, 2008. * * An example using this algorithm, initializing the grid with a constant lapse rate fill using +0.05% precipitation increase per meter of elevation, is given below: * @code * HNW::algorithms = HNW_SNOW * HNW::hnw_snow = cst_lapse * HNW::cst_lapse = 0.0005 frac * @endcode * * @author Florian Kobierska, Jan Magnusson and Mathias Bavay */ class SnowHNWInterpolation : public InterpolationAlgorithm { public: SnowHNWInterpolation(Meteo2DInterpolator& i_mi, const std::vector& i_vecArgs, const std::string& i_algo, IOManager& iom)  509  : InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom) {}  Mathias Bavay committed Sep 23, 2013 510  virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);  511  virtual void calculate(const DEMObject& dem, Grid2DObject& grid);  Mathias Bavay committed Sep 23, 2013 512 513 514 515 516 517 518 }; /** * @class OrdinaryKrigingAlgorithm * @brief Ordinary kriging. * This implements ordinary krigging (see https://secure.wikimedia.org/wikipedia/en/wiki/Kriging) * with user-selectable variogram model (see https://secure.wikimedia.org/wikipedia/en/wiki/Variogram).  Mathias Bavay committed Feb 25, 2014 519  * More details about the specific computation steps of kriging are provided in Interpol2D::ODKriging.  Mathias Bavay committed Sep 23, 2013 520 521 522 523 524  * * The variogram is currently computed with the current data (as 1/2*(X1-X2)^2), which makes it quite * uninteresting... The next improvement will consist in calculating the covariances (used to build the * variogram) from time series (thus reflecting the time-correlation between stations). *  Mathias Bavay committed Feb 25, 2014 525  * Please note that the variogram and krigging coefficients are re-computed fresh for each new grid (or time step).  Mathias Bavay committed Sep 23, 2013 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544  * The available variogram models are found in Fit1D::regression and given as optional arguments * (by default, LINVARIO is used). Several models can be given, the first that can fit the data will be used * for the current timestep: * @code * TA::algorithms = ODKRIG * TA::odkrig = SPHERICVARIO linvario * @endcode * * @author Mathias Bavay */ class OrdinaryKrigingAlgorithm : public InterpolationAlgorithm { public: OrdinaryKrigingAlgorithm(Meteo2DInterpolator& i_mi, const std::vector& i_vecArgs, const std::string& i_algo, IOManager& iom) : InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom), variogram() {} virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param); virtual void calculate(const DEMObject& dem, Grid2DObject& grid); protected:  545 546 547  size_t getTimeSeries(const bool& detrend_data, std::vector< std::vector > &vecVecData) const; void getDataForEmpiricalVariogram(std::vector &distData, std::vector &variData) const; void getDataForVariogram(std::vector &distData, std::vector &variData, const bool& detrend_data=false) const;  Mathias Bavay committed Jun 17, 2014 548  bool computeVariogram(const bool& detrend_data=false);  Mathias Bavay committed Sep 23, 2013 549 550 551 552 553 554 555 556 557 558  Fit1D variogram; }; /** * @class LapseOrdinaryKrigingAlgorithm * @brief Ordinary kriging with detrending. * This is very similar to OrdinaryKrigingAlgorithm but performs detrending on the data. * @code * TA::algorithms = ODKRIG_LAPSE  Mathias Bavay committed Jun 16, 2014 559  * TA::odkrig_lapse = SPHERICVARIO  Mathias Bavay committed Sep 23, 2013 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575  * @endcode * * @author Mathias Bavay */ class LapseOrdinaryKrigingAlgorithm : public OrdinaryKrigingAlgorithm { public: LapseOrdinaryKrigingAlgorithm(Meteo2DInterpolator& i_mi, const std::vector& i_vecArgs, const std::string& i_algo, IOManager& iom) : OrdinaryKrigingAlgorithm(i_mi, i_vecArgs, i_algo, iom) {} virtual void calculate(const DEMObject& dem, Grid2DObject& grid); }; } //end namespace mio #endif