WSL/SLF GitLab Repository

libinterpol2D.cc 26.6 KB
Newer Older
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 <http://www.gnu.org/licenses/>.
*/
18
19
#include <cmath>

20
#include <meteoio/meteostats/libinterpol2D.h>
21
#include <meteoio/meteolaws/Atmosphere.h>
22
#include <meteoio/meteolaws/Meteoconst.h> //for math constants
23
#include <meteoio/MathOptim.h> //math optimizations
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)
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
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)
73
{
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 );
81
82
}

83
/**
84
* @brief Build the list of (distance to grid cell, stations index) ordered by their distance to a grid cell
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)
88
89
90
*/
void Interpol2D::getNeighbors(const double& x, const double& y,
                              const std::vector<StationData>& vecStations,
91
                              std::vector< std::pair<double, size_t> >& list)
92
{
93
	list.resize(vecStations.size());
94

95
	for(size_t i=0; i<vecStations.size(); i++) {
96
97
98
99
		const Coords& position = vecStations[i].position;
		const double DX = x-position.getEasting();
		const double DY = y-position.getNorthing();
		const double d2 = (DX*DX + DY*DY);
100
		const std::pair <double, size_t> tmp(d2,i);
101
		list[i] = tmp;
102
103
104
105
106
	}

	sort (list.begin(), list.end());
}

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<StationData>& vecStations, std::vector<double>& vecEastings, std::vector<double>& vecNorthings)
{
	vecEastings.resize( vecStations.size() );
	vecNorthings.resize( vecStations.size() );
	for (size_t i=0; i<vecStations.size(); i++) {
		const Coords& position = vecStations[i].position;
		vecEastings[i] = position.getEasting();
		vecNorthings[i] = position.getNorthing();
	}
}

119
//these weighting functions take the square of a distance as an argument and return a weight
120
inline double Interpol2D::weightInvDist(const double& d2)
121
{
122
	return Optim::invSqrt( d2 ); //we use the optimized approximation for 1/sqrt
123
}
124
inline double Interpol2D::weightInvDistSqrt(const double& d2)
125
{
126
	return Optim::fastSqrt_Q3( Optim::invSqrt(d2) ); //we use the optimized approximation for 1/sqrt
127
}
128
inline double Interpol2D::weightInvDist2(const double& d2)
129
130
131
{
	return 1./d2; //we use the optimized approximation for 1/sqrt
}
132
inline double Interpol2D::weightInvDistN(const double& d2)
133
{
134
	return pow( Optim::invSqrt(d2) , dist_pow); //we use the optimized approximation for 1/sqrt
135
}
136
137
138

//Filling Functions
/**
139
* @brief Grid filling function:
140
141
142
143
* This implementation builds a standard air pressure as a function of the elevation
* @param dem array of elevations (dem)
* @param grid 2D array to fill
*/
144
145
void Interpol2D::stdPressure(const DEMObject& dem, Grid2DObject& grid) {
                             grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);
146
147

	//provide each point with an altitude dependant pressure... it is worth what it is...
148
149
	for (size_t j=0; j<grid.nrows; j++) {
		for (size_t i=0; i<grid.ncols; i++) {
150
			const double& cell_altitude=dem.grid2D(i,j);
151
			if (cell_altitude!=IOUtils::nodata) {
152
				grid.grid2D(i,j) = Atmosphere::stdAirPressure(cell_altitude);
153
154
155
156
157
158
159
160
161
162
163
164
165
166
			} else {
				grid.grid2D(i,j) = IOUtils::nodata;
			}
		}
	}
}

/**
* @brief Grid filling function:
* This implementation fills the grid with a constant value
* @param value value to put in the grid
* @param dem array of elevations (dem). This is needed in order to know if a point is "nodata"
* @param grid 2D array to fill
*/
167
void Interpol2D::constant(const double& value, const DEMObject& dem, Grid2DObject& grid)
168
169
170
171
{
	grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);

	//fills a data table with constant values
172
173
	for (size_t j=0; j<grid.nrows; j++) {
		for (size_t i=0; i<grid.ncols; i++) {
174
175
176
177
178
179
180
181
182
183
			if (dem.grid2D(i,j)!=IOUtils::nodata) {
				grid.grid2D(i,j) = value;
			} else {
				grid.grid2D(i,j) = IOUtils::nodata;
			}
		}
	}
}

double Interpol2D::IDWCore(const double& x, const double& y, const std::vector<double>& vecData_in,
184
                           const std::vector<double>& vecEastings, const std::vector<double>& vecNorthings)
185
186
{
	//The value at any given cell is the sum of the weighted contribution from each source
187
	const size_t n_stations=vecEastings.size();
188
	double parameter=0., norm=0.;
189
	const double scale = 1.e6;
190
191

	for (size_t i=0; i<n_stations; i++) {
192
193
		const double DX=x-vecEastings[i];
		const double DY=y-vecNorthings[i];
194
		const double weight = Optim::invSqrt( DX*DX + DY*DY + scale ); //use the optimized 1/sqrt approximation
195
196
197
198
199
200
		parameter += weight*vecData_in[i];
		norm += weight;
	}
	return (parameter/norm); //normalization
}

201
/*
202
* @brief Grid filling function:
203
* Similar to Interpol2D::LapseIDW but using a limited number of stations for each cell. We also assume a two segments regression for altitude detrending with
204
* a fixed 1200m above sea level inflection point.
205
206
207
208
209
210
211
* @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)
* @param nrOfNeighbors number of neighboring stations to use for each pixel
* @param grid 2D array to fill
* @param r2 average r² coefficient of the lapse rate regressions
*/
212
/*void Interpol2D::LocalLapseIDW(const std::vector<double>& vecData_in, const std::vector<StationData>& vecStations_in,
213
                               const DEMObject& dem, const size_t& nrOfNeighbors,
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
221
222
	for (size_t j=0; j<grid.nrows; j++) {
		for (size_t i=0; i<grid.ncols; i++) {
223
224
			//LL_IDW_pixel returns nodata when appropriate
			double r;
225
			const double value = LLIDW_pixel(i,j,vecData_in, vecStations_in, dem, nrOfNeighbors, r); //TODO: precompute in vectors
226
227
228
229
230
231
232
233
234
235
236
237
238
239
			grid.grid2D(i,j) = value;
			if(value!=IOUtils::nodata) {
				sum += fabs(r);
				count++;
			}
		}
	}
	if(count>0)
		r2 = sum/(double)count;
	else
		r2 = 0.;
}

//calculate a local pixel for LocalLapseIDW
240
double Interpol2D::LLIDW_pixel(const size_t& i, const size_t& j,
241
                                const std::vector<double>& vecData_in, const std::vector<StationData>& vecStations_in,
242
                                const DEMObject& dem, const size_t& nrOfNeighbors, double& r2)
243
244
245
246
247
{
	const double& cell_altitude=dem.grid2D(i,j);
	if(cell_altitude==IOUtils::nodata)
		return IOUtils::nodata;

248
	std::vector< std::pair<double, size_t> > list;
249
250
251
252
253
254
	std::vector<double> 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);
255
256
257
	const size_t max_stations = std::min(list.size(), nrOfNeighbors);
	for(size_t st=0; st<max_stations; st++) {
		const size_t st_index=list[st].second;
258
259
260
261
262
263
264
265
266
		const double value = vecData_in[st_index];
		const double alt = vecStations_in[st_index].position.getAltitude();
		if ((value != IOUtils::nodata) && (alt != IOUtils::nodata)) {
			X.push_back( alt );
			Y.push_back( value );
		}
	}

	//compute lapse rate
267
	if(X.empty())
268
		return IOUtils::nodata;
269
	coeffs.resize(7,0.);
270
	//BiLinRegression(X, Y, coeffs);
271
	r2=coeffs[3]*coeffs[6]; //Is it correct?
272
273
274
275
276

	//compute local pixel value
	unsigned int count=0;
	double pixel_value=0., norm=0.;
	const double scale=0.;
277
278
	for(size_t st=0; st<max_stations; st++) {
		const size_t st_index=list[st].second;
279
280
281
		const double value = vecData_in[st_index];
		const double alt = vecStations_in[st_index].position.getAltitude();
		if ((value != IOUtils::nodata) && (alt != IOUtils::nodata)) {
282
			//const double contrib = LinProject(value, alt, cell_altitude, coeffs);
Mathias Bavay's avatar
Mathias Bavay committed
283
			//const double contrib = BiLinProject(value, alt, cell_altitude, coeffs);
284
			const double weight = Optim::invSqrt( list[st].first + scale + 1.e-6 );
285
			//pixel_value += weight*contrib;
286
287
288
289
290
291
292
293
294
			norm += weight;
			count++;
		}
	}

	if(count>0)
		return (pixel_value/norm);
	else
		return IOUtils::nodata;
295
}*/
296
297

/**
298
* @brief Grid filling function:
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<double>& vecData_in, const std::vector<StationData>& vecStations_in,
                     const DEMObject& dem, Grid2DObject& grid)
{
	grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);
311
312
	std::vector<double> vecEastings, vecNorthings;
	buildPositionsVectors(vecStations_in, vecEastings, vecNorthings);
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's avatar
Mathias Bavay committed
318
319
320
321
	for (size_t jj=0; jj<grid.nrows; jj++) {
		for (size_t ii=0; ii<grid.ncols; ii++) {
			if (dem.grid2D(ii,jj)!=IOUtils::nodata) {
				grid.grid2D(ii,jj) = IDWCore((xllcorner+double(ii)*cellsize), (yllcorner+double(jj)*cellsize),
322
				                           vecData_in, vecEastings, vecNorthings);
323
			} else {
Mathias Bavay's avatar
Mathias Bavay committed
324
				grid.grid2D(ii,jj) = IOUtils::nodata;
325
326
327
328
329
330
			}
		}
	}
}

/**
331
* @brief Grid filling function:
332
333
* This implementation fills a grid using a curvature and slope algorithm, as described in "A Meteorological
* Distribution System for High-Resolution Terrestrial Modeling (MicroMet)", Liston and Elder, 2006.
334
* @param i_dem array of elevations (dem). The slope must have been updated as it is required for the DEM analysis.
335
336
337
* @param VW 2D array of Wind Velocity to fill
* @param DW 2D array of Wind Direction to fill
*/
338
void Interpol2D::SimpleDEMWindInterpolate(const DEMObject& i_dem, Grid2DObject& VW, Grid2DObject& DW)
339
{
340
	if ((!VW.isSameGeolocalization(DW)) || (!VW.isSameGeolocalization(i_dem))){
341
342
343
		throw IOException("Requested grid VW and grid DW don't match the geolocalization of the DEM", AT);
	}

344
345
346
347
348
349
350
351
352
353
	const bool recomputeDEM = i_dem.curvature.isEmpty();
	DEMObject *intern_dem = NULL;
	if(recomputeDEM) {
		std::cerr << "[W] WIND_CURV spatial interpolations algorithm selected but no dem curvature available! Computing it...\n";
		intern_dem = new DEMObject(i_dem);
		intern_dem->setUpdatePpt((DEMObject::update_type)(DEMObject::SLOPE|DEMObject::CURVATURE));
		intern_dem->update();
	}
	const DEMObject *dem = (recomputeDEM)? intern_dem : &i_dem;

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

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
372

373
374
	for (size_t j=0;j<VW.nrows;j++) {
		for (size_t i=0;i<VW.ncols;i++){
375
			speed = VW.grid2D(i,j);
376
			if(speed==0.) continue; //we can not apply any correction factor!
377
			dir = DW.grid2D(i,j);
378
379
380
			beta = dem->slope(i, j)*Cst::to_rad;
			azi = dem->azi(i, j)*Cst::to_rad;
			curvature = dem->curvature(i, j);
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;
388
				//Speed and direction converted to zonal et meridional
389
				//components
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);
395
				dir = (1.5 * Cst::PI) - atan(v/u);
396

397
				//normalize curvature and beta.
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;
402
403
404

				// Calculate the slope in the direction of the wind
				slopeDir = beta * cos(dir - azi);
405

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;
414

415
				// Add the diverting factor to the wind direction and convert to degrees
416
				DW.grid2D(i, j) = (dir + Od) * Cst::to_deg;
417
418
419
420
421
422
				if( DW.grid2D(i, j)>360. ) {
					DW.grid2D(i, j) -= 360.;
				}
			}
		}
	}
423
424

	if(intern_dem!=NULL) delete (intern_dem);
425
}
426

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
* <i>"Quantitative evaluation of different hydrological modelling approaches in a partly glacierized Swiss watershed"</i>, Magnusson et All., Hydrological Processes, 2010, under review.
431
* and
432
433
434
435
* <i>"Modelling runoff from highly glacierized alpine catchments in a changing climate"</i>, Huss et All., Hydrological Processes, <b>22</b>, 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
437
*/
438
void Interpol2D::CurvatureCorrection(DEMObject& dem, const Grid2DObject& ta, Grid2DObject& grid)
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;j<grid.nrows;j++) {
        for (size_t i=0;i<grid.ncols;i++) {
            const double slope = dem.slope(i, j);
            const double curvature = dem.curvature(i, j);
            double val = grid.grid2D(i, j);

            if(ta.grid2D(i,j)<=273.15) {
453
				//we only modify the grid of precipitations if air temperature
454
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
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
				//at this point is below or at freezing

                if(slope==IOUtils::nodata || curvature==IOUtils::nodata) continue;

                if(val!=IOUtils::nodata && dem_range_curvature!=0.) {
                    //cf Huss
                    grid.grid2D(i, j) = val*(0.5-(curvature-dem_max_curvature)/dem_range_curvature);
                }
            }
        }
    }

    //HACK: correction for precipitation sum over the whole domain
	//this is a cheap/crappy way of compensating for the spatial redistribution of snow on the slopes
	const double new_mean = grid.grid2D.getMean();
	if(new_mean!=0.) grid.grid2D *= orig_mean/new_mean;

}

void Interpol2D::steepestDescentDisplacement(const DEMObject& dem, const Grid2DObject& grid, const size_t& ii, const size_t& jj, short &d_i_dest, short &d_j_dest)
{
    double max_slope = 0.;
	d_i_dest = 0, d_j_dest = 0;

    //loop around all adjacent cells to find the cell with the steepest downhill slope
    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);
            const double new_slope = dem.slope(ii + d_i, jj + d_j);

            if ((new_slope>max_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
588
* This implements the matrix formulation of Ordinary Kriging, as shown (for example) in
589
590
591
592
* <i>"Statistics for spatial data"</i>, 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
593
* @param variogram variogram regression model
594
595
596
* @param grid 2D array of precipitation to fill
* @author Mathias Bavay
*/
597
void Interpol2D::ODKriging(const std::vector<double>& vecData, const std::vector<StationData>& 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();
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's avatar
Mathias Bavay committed
606
	Matrix Ginv(nrOfMeasurments+1, nrOfMeasurments+1);
607
608

	//fill the Ginv matrix
609
	for(size_t j=1; j<=nrOfMeasurments; j++) {
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
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);
620
			Ginv(i,j) = variogram.f(distance);
621
		}
622
623
		Ginv(j,j)=1.; //HACK diagonal should contain the nugget...
		Ginv(nrOfMeasurments+1,j) = 1.; //last line filled with 1s
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++) {
628
			Ginv(i,j) = Ginv(j,i);
629
630
		}
	}
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's avatar
Mathias Bavay committed
637
	Matrix G0(nrOfMeasurments+1, (size_t)1);
638
	//now, calculate each point
639
640
	for(size_t j=0; j<grid.nrows; j++) {
		for(size_t i=0; i<grid.ncols; i++) {
641
642
			const double x = llcorner_x+static_cast<double>(i)*cellsize;
			const double y = llcorner_y+static_cast<double>(j)*cellsize;
643
644

			//fill gamma
645
			for(size_t st=0; st<nrOfMeasurments; st++) {
646
647
648
649
				//compute distance between cell and each station
				const Coords& position = vecStations[st].position;
				const double DX = x-position.getEasting();
				const double DY = y-position.getNorthing();
650
				const double distance = Optim::fastSqrt_Q3(DX*DX + DY*DY);
651

652
				G0(st+1,1) = variogram.f(distance); //matrix starts at 1, not 0
653
			}
654
			G0(nrOfMeasurments+1,1) = 1.; //last value is always 1
655

656
			const Matrix lambda = Ginv*G0;
657
658
659

			//calculate local parameter interpolation
			double p = 0.;
660
			for(size_t st=0; st<nrOfMeasurments; st++) {
661
				p += lambda(st+1,1) * vecData[st]; //matrix starts at 1, not 0
662
663
664
665
			}
			grid.grid2D(i,j) = p;
		}
	}
666
}
667

668
} //namespace