WSL/SLF GitLab Repository

libinterpol2D.cc 26.6 KB
 Thomas Egger committed Jun 21, 2010 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 ``````/***********************************************************************************/ /* Copyright 2009 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 . */ `````` Mathias Bavay committed Mar 09, 2011 18 19 ``````#include `````` Mathias Bavay committed Jul 19, 2011 20 ``````#include `````` Mathias Bavay committed Jan 26, 2011 21 ``````#include `````` Mathias Bavay committed Jan 19, 2012 22 ``````#include //for math constants `````` 23 ``````#include //math optimizations `````` Thomas Egger committed Jun 21, 2010 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 `````` using namespace std; namespace mio { const double Interpol2D::wind_ys = 0.58; const double Interpol2D::wind_yc = 0.42; //Usefull functions /** * @brief Computes the horizontal distance between points, given by coordinates in a geographic grid * @param X1 (const double) first point's X coordinate * @param Y1 (const double) first point's Y coordinate * @param X2 (const double) second point's X coordinate * @param Y2 (const double) second point's Y coordinate * @return (double) distance in m */ `````` 40 ``````inline double Interpol2D::HorizontalDistance(const double& X1, const double& Y1, const double& X2, const double& Y2) `````` Thomas Egger committed Jun 21, 2010 41 42 43 ``````{ //This function computes the horizontaldistance between two points //coordinates are given in a square, metric grid system `````` 44 45 46 47 48 49 50 51 52 53 54 55 `````` const double DX=(X1-X2), DY=(Y1-Y2); return sqrt( DX*DX + DY*DY ); } /** * @brief Computes the 1/horizontal distance between points, given by coordinates in a geographic grid * @param X1 (const double) first point's X coordinate * @param Y1 (const double) first point's Y coordinate * @param X2 (const double) second point's X coordinate * @param Y2 (const double) second point's Y coordinate * @return (double) 1/distance in m */ `````` 56 ``````inline double Interpol2D::InvHorizontalDistance(const double& X1, const double& Y1, const double& X2, const double& Y2) `````` 57 58 59 60 ``````{ //This function computes 1/horizontaldistance between two points //coordinates are given in a square, metric grid system const double DX=(X1-X2), DY=(Y1-Y2); `````` 61 `````` return Optim::invSqrt( DX*DX + DY*DY ); //we use the optimized approximation for 1/sqrt `````` Thomas Egger committed Jun 21, 2010 62 63 64 65 66 67 68 69 70 71 ``````} /** * @brief Computes the horizontal distance between points, given by their cells indexes * @param X1 (const double) first point's i index * @param Y1 (const double) first point's j index * @param X2 (const double) second point's X coordinate * @param Y2 (const double) second point's Y coordinate * @return (double) distance in m */ `````` 72 ``````inline double Interpol2D::HorizontalDistance(const DEMObject& dem, const int& i, const int& j, const double& X2, const double& Y2) `````` Mathias Bavay committed Sep 16, 2013 73 ``````{ `````` Thomas Egger committed Jun 21, 2010 74 75 76 77 78 `````` //This function computes the horizontal distance between two points //coordinates are given in a square, metric grid system //for grid points toward real coordinates const double X1 = (dem.llcorner.getEasting()+i*dem.cellsize); const double Y1 = (dem.llcorner.getNorthing()+j*dem.cellsize); `````` 79 80 `````` const double DX=(X1-X2), DY=(Y1-Y2); return sqrt( DX*DX + DY*DY ); `````` Thomas Egger committed Jun 21, 2010 81 82 ``````} `````` Mathias Bavay committed Oct 27, 2010 83 ``````/** `````` 84 ``````* @brief Build the list of (distance to grid cell, stations index) ordered by their distance to a grid cell `````` Mathias Bavay committed Oct 27, 2010 85 86 ``````* @param x x coordinate of cell * @param y y coordinate of cell `````` 87 ``````* @param list list of pairs (distance to grid cell, stations index) `````` Mathias Bavay committed Oct 27, 2010 88 89 90 ``````*/ void Interpol2D::getNeighbors(const double& x, const double& y, const std::vector& vecStations, `````` Mathias Bavay committed Nov 10, 2011 91 `````` std::vector< std::pair >& list) `````` Mathias Bavay committed Oct 27, 2010 92 ``````{ `````` Mathias Bavay committed Mar 10, 2013 93 `````` list.resize(vecStations.size()); `````` Mathias Bavay committed May 17, 2011 94 `````` `````` Mathias Bavay committed Nov 10, 2011 95 `````` for(size_t i=0; i tmp(d2,i); `````` Mathias Bavay committed Mar 10, 2013 101 `````` list[i] = tmp; `````` Mathias Bavay committed Oct 27, 2010 102 103 104 105 106 `````` } sort (list.begin(), list.end()); } `````` Mathias Bavay committed Jul 25, 2013 107 108 109 110 111 112 113 114 115 116 117 118 ``````//convert a vector of stations into two vectors of eastings and northings void Interpol2D::buildPositionsVectors(const std::vector& vecStations, std::vector& vecEastings, std::vector& vecNorthings) { vecEastings.resize( vecStations.size() ); vecNorthings.resize( vecStations.size() ); for (size_t i=0; i& vecData_in, `````` Mathias Bavay committed Jan 12, 2012 184 `````` const std::vector& vecEastings, const std::vector& vecNorthings) `````` Thomas Egger committed Jun 21, 2010 185 186 ``````{ //The value at any given cell is the sum of the weighted contribution from each source `````` Mathias Bavay committed Jan 12, 2012 187 `````` const size_t n_stations=vecEastings.size(); `````` 188 `````` double parameter=0., norm=0.; `````` Mathias Bavay committed Aug 19, 2010 189 `````` const double scale = 1.e6; `````` Mathias Bavay committed May 17, 2011 190 191 `````` for (size_t i=0; i& vecData_in, const std::vector& vecStations_in, `````` Mathias Bavay committed Nov 10, 2011 213 `````` const DEMObject& dem, const size_t& nrOfNeighbors, `````` Mathias Bavay committed Oct 27, 2010 214 215 216 217 218 219 220 `````` Grid2DObject& grid, double& r2) { unsigned int count=0; double sum=0; grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner); //run algorithm `````` Mathias Bavay committed Sep 18, 2013 221 222 `````` for (size_t j=0; j0) r2 = sum/(double)count; else r2 = 0.; } //calculate a local pixel for LocalLapseIDW `````` Mathias Bavay committed Sep 18, 2013 240 ``````double Interpol2D::LLIDW_pixel(const size_t& i, const size_t& j, `````` Mathias Bavay committed Oct 27, 2010 241 `````` const std::vector& vecData_in, const std::vector& vecStations_in, `````` Mathias Bavay committed Nov 10, 2011 242 `````` const DEMObject& dem, const size_t& nrOfNeighbors, double& r2) `````` Mathias Bavay committed Oct 27, 2010 243 244 245 246 247 ``````{ const double& cell_altitude=dem.grid2D(i,j); if(cell_altitude==IOUtils::nodata) return IOUtils::nodata; `````` Mathias Bavay committed Nov 10, 2011 248 `````` std::vector< std::pair > list; `````` Mathias Bavay committed Oct 27, 2010 249 250 251 252 253 254 `````` std::vector X, Y, coeffs; //fill vectors with appropriate neighbors const double x = dem.llcorner.getEasting()+i*dem.cellsize; const double y = dem.llcorner.getNorthing()+j*dem.cellsize; getNeighbors(x, y, vecStations_in, list); `````` Mathias Bavay committed Nov 10, 2011 255 256 257 `````` const size_t max_stations = std::min(list.size(), nrOfNeighbors); for(size_t st=0; st0) return (pixel_value/norm); else return IOUtils::nodata; `````` Mathias Bavay committed Jul 25, 2013 295 ``````}*/ `````` Mathias Bavay committed Oct 27, 2010 296 297 `````` /** `````` Mathias Bavay committed May 17, 2011 298 ``````* @brief Grid filling function: `````` Thomas Egger committed Jun 21, 2010 299 300 301 302 303 304 305 306 307 308 309 310 ``````* This implementation fills a grid using Inverse Distance Weighting. * for example, the air temperatures measured at several stations would be given as values, the stations positions * as positions and projected to a grid. No elevation detrending is performed, the DEM is only used for checking if a grid point is "nodata". * @param vecData_in input values to use for the IDW * @param vecStations_in position of the "values" (altitude and coordinates) * @param dem array of elevations (dem). This is needed in order to know if a point is "nodata" * @param grid 2D array to fill */ void Interpol2D::IDW(const std::vector& vecData_in, const std::vector& vecStations_in, const DEMObject& dem, Grid2DObject& grid) { grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner); `````` Mathias Bavay committed Jan 12, 2012 311 312 `````` std::vector vecEastings, vecNorthings; buildPositionsVectors(vecStations_in, vecEastings, vecNorthings); `````` Thomas Egger committed Jun 21, 2010 313 314 315 316 317 `````` //multiple source stations: simple IDW Krieging const double xllcorner = dem.llcorner.getEasting(); const double yllcorner = dem.llcorner.getNorthing(); const double cellsize = dem.cellsize; `````` Mathias Bavay committed Sep 18, 2013 318 319 320 321 `````` for (size_t jj=0; jjsetUpdatePpt((DEMObject::update_type)(DEMObject::SLOPE|DEMObject::CURVATURE)); intern_dem->update(); } const DEMObject *dem = (recomputeDEM)? intern_dem : &i_dem; `````` Thomas Egger committed Jun 21, 2010 354 355 356 357 358 359 360 361 362 363 364 `````` //This method computes the speed of the wind and returns a table in 2D with this values double speed; // Wind speed (m s-1) double dir; // Wind direction double u; // Zonal component u (m s-1) double v; // Meridional component v (m s-1) double beta; // Terrain slope double azi; // Topographic slope azimuth double curvature; // Topographic curvature double slopeDir; // Slope in the direction of the wind double Ww; // Wind weighting double Od; // Diverting factor `````` 365 `````` `````` Mathias Bavay committed Sep 01, 2013 366 367 368 369 `````` const double dem_min_slope=dem->min_slope*Cst::to_rad; const double dem_min_curvature=dem->min_curvature; double dem_range_slope=(dem->max_slope-dem_min_slope)*Cst::to_rad; double dem_range_curvature=(dem->max_curvature-dem_min_curvature); `````` 370 371 `````` if(dem_range_slope==0.) dem_range_slope = 1.; //to avoid division by zero below if(dem_range_curvature==0.) dem_range_curvature = 1.; //to avoid division by zero below `````` Mathias Bavay committed May 17, 2011 372 `````` `````` Mathias Bavay committed Sep 18, 2013 373 374 `````` for (size_t j=0;jslope(i, j)*Cst::to_rad; azi = dem->azi(i, j)*Cst::to_rad; curvature = dem->curvature(i, j); `````` Thomas Egger committed Jun 21, 2010 381 382 383 384 385 386 `````` if(speed==IOUtils::nodata || dir==IOUtils::nodata || beta==IOUtils::nodata || azi==IOUtils::nodata || curvature==IOUtils::nodata) { VW.grid2D(i, j) = IOUtils::nodata; DW.grid2D(i, j) = IOUtils::nodata; } else { //convert direction to rad `````` 387 `````` dir *= Cst::to_rad; `````` Thomas Egger committed Jun 21, 2010 388 `````` //Speed and direction converted to zonal et meridional `````` Mathias Bavay committed May 17, 2011 389 `````` //components `````` Thomas Egger committed Jun 21, 2010 390 391 392 393 394 `````` u = (-1.) * (speed * sin(dir)); v = (-1.) * (speed * cos(dir)); // Converted back to speed and direction speed = sqrt(u*u + v*v); `````` Mathias Bavay committed Jan 19, 2012 395 `````` dir = (1.5 * Cst::PI) - atan(v/u); `````` Thomas Egger committed Jun 21, 2010 396 `````` `````` Mathias Bavay committed May 17, 2011 397 `````` //normalize curvature and beta. `````` Thomas Egger committed Jun 21, 2010 398 399 `````` //Note: it should be slopeDir instead of beta, but beta is more efficient //to compute (only once for each dem) and it should not be that different... `````` 400 401 `````` beta = (beta - dem_min_slope)/dem_range_slope - 0.5; curvature = (curvature - dem_min_curvature)/dem_range_curvature - 0.5; `````` Thomas Egger committed Jun 21, 2010 402 403 404 `````` // Calculate the slope in the direction of the wind slopeDir = beta * cos(dir - azi); `````` Mathias Bavay committed May 17, 2011 405 `````` `````` Thomas Egger committed Jun 21, 2010 406 407 408 409 410 411 412 413 `````` // Calculate the wind weighting factor Ww = 1. + wind_ys * slopeDir + wind_yc * curvature; // Modify the wind direction by a diverting factor Od = -0.5 * slopeDir * sin(2.*(azi - dir)); // Calculate the terrain-modified wind speed VW.grid2D(i, j) = Ww * speed; `````` Mathias Bavay committed Mar 29, 2012 414 `````` `````` Thomas Egger committed Jun 21, 2010 415 `````` // Add the diverting factor to the wind direction and convert to degrees `````` 416 `````` DW.grid2D(i, j) = (dir + Od) * Cst::to_deg; `````` Thomas Egger committed Jun 21, 2010 417 418 419 420 421 422 `````` if( DW.grid2D(i, j)>360. ) { DW.grid2D(i, j) -= 360.; } } } } `````` Mathias Bavay committed Sep 01, 2013 423 424 `````` if(intern_dem!=NULL) delete (intern_dem); `````` 425 ``````} `````` Thomas Egger committed Jun 21, 2010 426 `````` `````` Mathias Bavay committed Sep 30, 2010 427 428 429 430 ``````/** * @brief Distribute precipitation in a way that reflects snow redistribution on the ground, according to (Huss, 2008) * This method modifies the solid precipitation distribution according to the local slope and curvature. See * "Quantitative evaluation of different hydrological modelling approaches in a partly glacierized Swiss watershed", Magnusson et All., Hydrological Processes, 2010, under review. `````` Mathias Bavay committed May 17, 2011 431 ``````* and `````` Mathias Bavay committed Sep 30, 2010 432 433 434 435 ``````* "Modelling runoff from highly glacierized alpine catchments in a changing climate", Huss et All., Hydrological Processes, 22, 3888-3902, 2008. * @param dem array of elevations (dem). The slope must have been updated as it is required for the DEM analysis. * @param ta array of air temperatures used to determine if precipitation is rain or snow * @param grid 2D array of precipitation to fill `````` 436 ``````* @author Florian Kobierska, Jan Magnusson, Rob Spence and Mathias Bavay `````` Mathias Bavay committed Sep 30, 2010 437 ``````*/ `````` 438 ``````void Interpol2D::CurvatureCorrection(DEMObject& dem, const Grid2DObject& ta, Grid2DObject& grid) `````` Mathias Bavay committed Sep 30, 2010 439 440 441 442 ``````{ if(!grid.isSameGeolocalization(dem)) { throw IOException("Requested grid does not match the geolocalization of the DEM", AT); } `````` 443 444 445 446 447 448 449 450 451 452 `````` const double dem_max_curvature = dem.max_curvature, dem_range_curvature=(dem.max_curvature-dem.min_curvature); const double orig_mean = grid.grid2D.getMean(); for (size_t j=0;jmax_slope) && (height_ratio>1.)){ max_slope = new_slope; d_i_dest = d_i; d_j_dest = d_j; } } } } double Interpol2D::depositAroundCell(const DEMObject& dem, const size_t& ii, const size_t& jj, const double& precip, Grid2DObject &grid) { //else add precip to the cell and remove the same amount from the precip variable grid.grid2D(ii, jj) += precip; double distributed_precip = precip; for(short d_i=-1;d_i<=1;d_i++){ for(short d_j=-1;d_j<=1;d_j++){ const double elev_pt1 = dem.grid2D(ii, jj); const double elev_pt2 = dem.grid2D(ii + d_i, jj + d_j); const double precip_1 = grid.grid2D(ii, jj); const double precip_2 = grid.grid2D(ii + d_i, jj + d_j); const double height_ratio = (elev_pt1+precip_1) / (elev_pt2+precip_2); if ((d_i!=0)||(d_j!=0)){ if (height_ratio>1.){ grid.grid2D(ii + d_i, jj + d_j) += precip; distributed_precip += precip; } } } } return distributed_precip; } /** *@brief redistribute precip from steeper slopes to gentler slopes by following the steepest path from top to bottom *and gradually depositing precip during descent *@param precip variable to keep track of removed and deposited precip in order to maintain a (roughly) constant global sum *@param counter variable that is increased during descent in order to deposit more precip towards the foot *of the slope than at the top *@param height_ratio ratio of the height of the dem plus accumulated precip height between one cell and the next *in order to determine a route down the slope *@author Rob Spence and Mathias Bavay */ void Interpol2D::SteepSlopeRedistribution(const DEMObject& dem, const Grid2DObject& ta, Grid2DObject& grid) { for (size_t jj=1; jj<(grid.nrows-1); jj++) { for (size_t ii=1; ii<(grid.ncols-1); ii++) { if(grid.grid2D(ii,jj)==IOUtils::nodata) continue; //we only modify the grid of precipitations if air temperature //at this point is below or at freezing if(ta.grid2D(ii, jj)<=273.15) { const double slope = dem.slope(ii, jj); const double curvature = dem.curvature(ii, jj); if(slope==IOUtils::nodata || curvature==IOUtils::nodata) continue; if (slope>40.) { //redistribution only above 40 degrees //remove all precip above 60 deg or linearly decrease it double precip = (slope>60.)? grid.grid2D(ii, jj) : grid.grid2D(ii, jj) * ((40.-slope)/-30.); grid.grid2D(ii, jj) -= precip; //we will redistribute the precipitation in a different way const double increment = precip / 50.; //break removed precip into smaller amounts to be redistributed double counter = 0.5; //counter will determine amount of precip deposited size_t ii_dest = ii, jj_dest = jj; while(precip>0.) { short d_i, d_j; steepestDescentDisplacement(dem, grid, ii_dest, jj_dest, d_i, d_j); //move to the destination cell ii_dest += d_i; jj_dest += d_j; //if (((ii_dest+d_i)<0) || ((jj_dest+d_j)<0) || ((ii_dest+d_i)>=grid.ncols)|| ((jj_dest+d_j)>=grid.nrows)) { if ((ii_dest==0) || (jj_dest==0) || (ii_dest==(grid.ncols-1))|| (jj_dest==(grid.nrows-1))){ //we are getting out of the domain: deposit local contribution grid.grid2D(ii_dest, jj_dest) += counter*increment; break; } if(d_i==0 && d_j==0) { //local minimum, everything stays here... grid.grid2D(ii_dest, jj_dest) += precip; break; } else {precip -= depositAroundCell(dem, ii_dest, jj_dest, counter*increment, grid); counter += 0.25; //greater amount of precip is deposited as we move down the slope } } } } } } } `````` 586 587 ``````/** * @brief Ordinary Kriging matrix formulation `````` Mathias Bavay committed May 17, 2011 588 ``````* This implements the matrix formulation of Ordinary Kriging, as shown (for example) in `````` 589 590 591 592 ``````* "Statistics for spatial data", Noel A. C. Cressie, John Wiley & Sons, revised edition, 1993, pp122. * @param vecData vector containing the values as measured at the stations * @param vecStations vector of stations * @param dem digital elevation model `````` Mathias Bavay committed Aug 16, 2011 593 ``````* @param variogram variogram regression model `````` 594 595 596 ``````* @param grid 2D array of precipitation to fill * @author Mathias Bavay */ `````` Mathias Bavay committed Aug 16, 2011 597 ``````void Interpol2D::ODKriging(const std::vector& vecData, const std::vector& vecStations, const DEMObject& dem, const Fit1D& variogram, Grid2DObject& grid) `````` 598 599 ``````{ grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner); `````` 600 `````` size_t nrOfMeasurments = vecStations.size(); `````` Mathias Bavay committed Aug 16, 2011 601 602 603 604 605 `````` //precompute various coordinates in the grid const double llcorner_x = grid.llcorner.getEasting(); const double llcorner_y = grid.llcorner.getNorthing(); const double cellsize = grid.cellsize; `````` Mathias Bavay committed Sep 13, 2013 606 `````` Matrix Ginv(nrOfMeasurments+1, nrOfMeasurments+1); `````` Mathias Bavay committed Nov 10, 2011 607 608 `````` //fill the Ginv matrix `````` 609 `````` for(size_t j=1; j<=nrOfMeasurments; j++) { `````` Mathias Bavay committed Aug 16, 2011 610 `````` const Coords& st1 = vecStations[j-1].position; `````` 611 612 613 `````` const double x1 = st1.getEasting(); const double y1 = st1.getNorthing(); `````` 614 `````` for(size_t i=1; i<=j; i++) { `````` 615 `````` //compute distance between stations `````` Mathias Bavay committed Aug 16, 2011 616 `````` const Coords& st2 = vecStations[i-1].position; `````` 617 618 `````` const double DX = x1-st2.getEasting(); const double DY = y1-st2.getNorthing(); `````` 619 `````` const double distance = Optim::fastSqrt_Q3(DX*DX + DY*DY); `````` Mathias Bavay committed Nov 10, 2011 620 `````` Ginv(i,j) = variogram.f(distance); `````` Mathias Bavay committed Aug 16, 2011 621 `````` } `````` Mathias Bavay committed Nov 10, 2011 622 623 `````` Ginv(j,j)=1.; //HACK diagonal should contain the nugget... Ginv(nrOfMeasurments+1,j) = 1.; //last line filled with 1s `````` Mathias Bavay committed Aug 16, 2011 624 625 `````` } //fill the upper half (an exact copy of the lower half) `````` 626 627 `````` for(size_t j=1; j<=nrOfMeasurments; j++) { for(size_t i=j+1; i<=nrOfMeasurments; i++) { `````` Mathias Bavay committed Nov 10, 2011 628 `````` Ginv(i,j) = Ginv(j,i); `````` 629 630 `````` } } `````` Mathias Bavay committed Nov 10, 2011 631 632 633 634 635 `````` //add last column of 1's and a zero for(size_t i=1; i<=nrOfMeasurments; i++) Ginv(i,nrOfMeasurments+1) = 1.; Ginv(nrOfMeasurments+1,nrOfMeasurments+1) = 0.; //invert the matrix Ginv.inv(); `````` 636 `````` `````` Mathias Bavay committed Sep 13, 2013 637 `````` Matrix G0(nrOfMeasurments+1, (size_t)1); `````` 638 `````` //now, calculate each point `````` 639 640 `````` for(size_t j=0; j(i)*cellsize; const double y = llcorner_y+static_cast(j)*cellsize; `````` 643 644 `````` //fill gamma `````` 645 `````` for(size_t st=0; st