1 /**********************************************************
2 * Version $Id$
3 *********************************************************/
4
5 ///////////////////////////////////////////////////////////
6 // //
7 // SAGA //
8 // //
9 // System for Automated Geoscientific Analyses //
10 // //
11 // Tool Library //
12 // ta_slope_stability //
13 // //
14 //-------------------------------------------------------//
15 // //
16 // WETNESS_01.cpp //
17 // //
18 // Copyright (C) 2012 by //
19 // Andreas G�nther //
20 // //
21 //-------------------------------------------------------//
22 // //
23 // This file is part of 'SAGA - System for Automated //
24 // Geoscientific Analyses'. SAGA is free software; you //
25 // can redistribute it and/or modify it under the terms //
26 // of the GNU General Public License as published by the //
27 // Free Software Foundation, either version 2 of the //
28 // License, or (at your option) any later version. //
29 // //
30 // SAGA is distributed in the hope that it will be //
31 // useful, but WITHOUT ANY WARRANTY; without even the //
32 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
33 // PARTICULAR PURPOSE. See the GNU General Public //
34 // License for more details. //
35 // //
36 // You should have received a copy of the GNU General //
37 // Public License along with this program; if not, see //
38 // <http://www.gnu.org/licenses/>. //
39 // //
40 //-------------------------------------------------------//
41 // //
42 // e-mail: a.guenther@bgr.de //
43 // //
44 // contact: Andreas G�nther //
45 // B2.2 //
46 // BGR //
47 // Germany //
48 // //
49 ///////////////////////////////////////////////////////////
50
51 //---------------------------------------------------------
52
53
54 ///////////////////////////////////////////////////////////
55 // //
56 // //
57 // //
58 ///////////////////////////////////////////////////////////
59
60 //---------------------------------------------------------
61 #include "WETNESS_01.h"
62 #include <math.h>
63
64 //---------------------------------------------------------
65 #define RUN_TOOL(LIBRARY, TOOL, CONDITION) {\
66 bool bResult;\
67 SG_RUN_TOOL(bResult, LIBRARY, TOOL, CONDITION)\
68 if( !bResult ) return( false );\
69 }
70
71 #define SET_PARAMETER(IDENTIFIER, VALUE) pTool->Get_Parameters()->Set_Parameter(SG_T(IDENTIFIER), VALUE)
72
73 //---------------------------------------------------------
74
75
76 ///////////////////////////////////////////////////////////
77 // //
78 // //
79 // //
80 ///////////////////////////////////////////////////////////
81
82 //---------------------------------------------------------
CWETNESS(void)83 CWETNESS::CWETNESS(void)
84 {
85 CSG_Parameters P;
86
87 //-----------------------------------------------------
88
89 Set_Name (_TL("WETNESS"));
90
91 Set_Author (SG_T("A. Günther (c) 2012"));
92
93 Set_Description (_TW(
94 "\n"
95 "This tool calculates a topographic wetness index (TWI) following Montgomery & Dietrich (1994) that can be used to estimate the degree of saturation of unconsolidated, permeable "
96 "materials above (more or less) impermeable bedrock. In contrast to the common TOPMODEL (Beven & Kirkby, 1979) - based TWI, this index differs in such that it considers "
97 "hydraulic conductivity to be constant in a soil mantle overlying relatively impermeable bedrock. Also, it uses the sine of the slope rather than its tangens, which is more correct "
98 "and significantly matters for steeper slopes that give raise to landslides. "
99 "For computation, a slope (in radians) and a catchment area (in m2) grid are required. "
100 "Additionally, information on groundwater recharge (m/hr), material hydraulic conductivity (m/hr), "
101 "and depth to potential shear plane (m) are required that can be specified either globally or through grids. "
102 "The tool produces a continuous wetness index (-) where cells with WI values > 1 (overland flow) set to 1, and optionally creates a classified WI grid rendering three saturation classes:.\n"
103 "0): Low moisture (WI smaller 0.1)\n"
104 "1): Partially wet (0.1 smaller WI smaller 1)\n"
105 "2): Saturation zone (WI larger 1)\n"
106 "\n"
107 "References:\n"
108 "<a href=\"http://www.tandfonline.com/doi/abs/10.1080/02626667909491834\">Beven, K.J., Kirkby, M.J. (1979) A physically-based variable contributing area model of basin hydrology. Hydrology Science Bulletin, 24, 43-69.</a>.\n"
109 "\n"
110 "<a href=\"http://www.agu.org/pubs/crossref/1994/93WR02979.shtml\">Montgomery D. R., Dietrich, W. E. (1994) A physically based model for the topographic control on shallow landsliding. Water Resources Research, 30, 1153-1171.</a>.\n"
111
112 ));
113
114 Parameters.Add_Grid(
115 NULL, "DEM", "DEM", "A DEM", PARAMETER_INPUT
116 );
117
118 /*Parameters.Add_Grid(
119 NULL, "B", "Catchment area grid (m2)", "A catchment area grid (in square meters)", PARAMETER_INPUT
120 );
121 */
122 Parameters.Add_Grid(
123 NULL, "Cmin", "Min hydraulic conductivity grid (m/hr) ", "A grid representing minimum material hydraulic conductivity (in m/hr)", PARAMETER_INPUT_OPTIONAL
124 );
125
126 Parameters.Add_Grid(
127 NULL, "Cmax", "Max hydraulic conductivity grid (m/hr) ", "A grid representing maximum material hydraulic conductivity (in m/hr)", PARAMETER_INPUT_OPTIONAL
128 );
129
130 Parameters.Add_Grid(
131 NULL, "Dmin", "Min groundwater recharge grid (m/hr) ", "A grid representing minimum groundwater recharge (in m/hr)", PARAMETER_INPUT_OPTIONAL
132 );
133
134 Parameters.Add_Grid(
135 NULL, "Dmax", "Max groundwater recharge grid (m/hr) ", "A grid representing maximum groundwater recharge (in m/hr)", PARAMETER_INPUT_OPTIONAL
136 );
137
138 Parameters.Add_Grid(
139 NULL, "Emin", "Min material depth grid (m)", "A grid representing minimum depth to potential shear plane (in m)", PARAMETER_INPUT_OPTIONAL
140 );
141
142 Parameters.Add_Grid(
143 NULL, "Emax", "Max material depth grid (m)", "A grid representing maximum depth to potential shear plane (in m)", PARAMETER_INPUT_OPTIONAL
144 );
145
146 Parameters.Add_Value(
147 NULL, "fCmin", "Min global material conductivity (m/hr)", "Constant value if no raster set", PARAMETER_TYPE_Double, 2.7 //Initialisierung eines fixen wertes vs Grid f�r density
148 );
149
150 Parameters.Add_Value(
151 NULL, "fCmax", "Max global material conductivity (m/hr)", "Constant value if no raster set", PARAMETER_TYPE_Double, 2.7 //Initialisierung eines fixen wertes vs Grid f�r density
152 );
153
154 Parameters.Add_Value(
155 NULL, "fDmin", "Min global groundwater recharge (m/hr)", "Constant value if no raster set", PARAMETER_TYPE_Double, 0.001 //Initialisierung eines fixen wertes vs Grid f�r conductivity
156 );
157
158 Parameters.Add_Value(
159 NULL, "fDmax", "Max global groundwater recharge (m/hr)", "Constant value if no raster set", PARAMETER_TYPE_Double, 0.001 //Initialisierung eines fixen wertes vs Grid f�r conductivity
160 );
161
162 Parameters.Add_Value(
163 NULL, "fEmin", "Min global material depth (m)", "Constant value if no raster set", PARAMETER_TYPE_Double, 1.0 //Initialisierung eines fixen wertes vs Grid f�r depth
164 );
165
166 Parameters.Add_Value(
167 NULL, "fEmax", "Max global material depth (m)", "Constant value if no raster set", PARAMETER_TYPE_Double, 1.0 //Initialisierung eines fixen wertes vs Grid f�r depth
168 );
169
170 Parameters.Add_Value(
171 NULL, "fH", "Parameter sampling runs", "Number of sampling cycles",PARAMETER_TYPE_Int, 1 //sampling cycles
172 );
173
174 Parameters.Add_Grid(
175 NULL, "F", "WI values", "Resulting wetness index (-) grid", PARAMETER_OUTPUT
176 );
177
178 Parameters.Add_Grid(
179 NULL, "G", "WI classes", "Classified wetness (-) grid", PARAMETER_OUTPUT_OPTIONAL
180 );
181
182 Parameters.Add_Choice(
183 NULL , "METHOD" , _TL("Catchment Area Calculation"),
184 _TL(""),
185 CSG_String::Format(SG_T("%s|%s|%s|%s|%s|%s|"),
186 _TL("Deterministic 8"),
187 _TL("Rho 8"),
188 _TL("Braunschweiger Reliefmodell"),
189 _TL("Deterministic Infinity"),
190 _TL("Multiple Flow Direction"),
191 _TL("Multiple Triangular Flow Directon")
192 ), 4
193 );
194
195 Parameters.Add_Value(
196 NULL , "PREPROC" , _TL("Preprocessing"),
197 _TL(""),
198 PARAMETER_TYPE_Bool, false
199 );
200 }
201
202
203 ///////////////////////////////////////////////////////////
204 // //
205 // //
206 // //
207 ///////////////////////////////////////////////////////////
208
209 //---------------------------------------------------------
210 enum
211 {
212 WI_NODATA = 0,
213 WI_LOW_MOISTURE,
214 WI_PARTIALLY_WET,
215 WI_SATURATION_ZONE,
216 WI_COUNT
217 };
218
On_Execute(void)219 bool CWETNESS::On_Execute(void)
220 {
221 double fCmin = Parameters("fCmin")->asDouble();
222 double fDmin = Parameters("fDmin")->asDouble();
223 double fEmin = Parameters("fEmin")->asDouble();
224 double fCmax = Parameters("fCmax")->asDouble();
225 double fDmax = Parameters("fDmax")->asDouble();
226 double fEmax = Parameters("fEmax")->asDouble();
227 double fH = Parameters("fH")->asInt();
228
229 CSG_Grid *pDEM, *pB, *pCmin, *pDmin, *pEmin, *pCmax, *pDmax, *pEmax,*pF, *pG;
230
231 pDEM = Parameters("DEM" )->asGrid(); //DEM
232 pCmin = Parameters("Cmin" )->asGrid(); //conductivity
233 pDmin = Parameters("Dmin" )->asGrid(); //recharge
234 pEmin = Parameters("Emin" )->asGrid(); //depth
235 pCmax = Parameters("Cmax" )->asGrid(); //conductivity
236 pDmax = Parameters("Dmax" )->asGrid(); //recharge
237 pEmax = Parameters("Emax" )->asGrid(); //depth
238 pF = Parameters("F" )->asGrid(); //output wetness index
239 pG = Parameters("G" )->asGrid(); //output wetness classes
240
241 //-----------------------------------------------------
242 // get catchment area sizes
243 CSG_Grid B(Get_System(), SG_DATATYPE_Float);
244
245 pB = &B;
246
247 if( Parameters("PREPROC")->asBool() )
248 {
249 CSG_Grid DEM(Get_System(), SG_DATATYPE_Float);
250
251 RUN_TOOL("ta_preprocessor" , 2,
252 SET_PARAMETER("DEM" , pDEM)
253 && SET_PARAMETER("DEM_PREPROC" , &DEM)
254 )
255
256 RUN_TOOL("ta_hydrology" , 0,
257 SET_PARAMETER("ELEVATION" , &DEM)
258 && SET_PARAMETER("FLOW" , pB)
259 && SET_PARAMETER("METHOD" , Parameters("METHOD"))
260 )
261 }
262 else
263 {
264 RUN_TOOL("ta_hydrology" , 0,
265 SET_PARAMETER("ELEVATION" , pDEM)
266 && SET_PARAMETER("FLOW" , pB)
267 && SET_PARAMETER("METHOD" , Parameters("METHOD"))
268 )
269 }
270
271
272 //-----------------------------------------------------
273 for(int y=0; y<Get_NY() && Set_Progress(y); y++)
274 {
275 #pragma omp parallel for
276 for(int x=0; x<Get_NX(); x++)
277 {
278 double a, b, c, d, e, f; //a slope, b catchment area
279 double cmin, dmin, emin;
280 double cmax, dmax, emax;
281 double cc, dd, ee;
282 int cperc, dperc, eperc;
283 int rand_int, h, n;
284
285 b = pB->asDouble(x, y); //Abfrage ob Raster oder Globalwerte:
286 cmin = pCmin ? pCmin->asDouble(x, y) : fCmin;
287 dmin = pDmin ? pDmin->asDouble(x, y) : fDmin;
288 emin = pEmin ? pEmin->asDouble(x, y) : fEmin;
289 cmax = pCmax ? pCmax->asDouble(x, y) : fCmax;
290 dmax = pDmax ? pDmax->asDouble(x, y) : fDmax;
291 emax = pEmax ? pEmax->asDouble(x, y) : fEmax;
292 h = fH;
293
294 if (pDEM->Get_Gradient(x, y, a, ee)==false
295 || (pCmin && pCmin->is_NoData(x, y))
296 || (pCmax && pCmax->is_NoData(x, y))
297 || (pDmin && pDmin->is_NoData(x, y))
298 || (pDmax && pDmax->is_NoData(x, y))
299 || (pEmin && pEmin->is_NoData(x, y))
300 || (pEmax && pEmax->is_NoData(x, y)) )
301 {
302 pF->Set_NoData(x, y);
303
304 if (pG)
305
306 pG->Set_Value(x, y, WI_NODATA);
307 }
308
309 else
310 {
311
312 cperc = ((cmax - cmin) / cmax) * 100; //calculate parameter range %: conductivity
313 if (cperc > 0)
314 {
315 n = 0;
316 cc = 0;
317 while ( n < h) //loop through specified random number iterations:
318 {
319 rand_int = rand() % cperc + 0; //calculate random percentage
320 c = ((cmax/100) * rand_int) + cmin; //calculate value
321 cc = cc + c; //sum
322 n = n + 1;
323 }
324 c = cc / n; // calculate mean from random values
325 }
326 else
327 {
328 c = cmax;
329 }
330
331 dperc = ((dmax - dmin) / dmax) * 100; //calculate parameter range %: recharge
332 if (dperc > 0)
333 {
334 n = 0;
335 dd = 0;
336 while ( n < h) //loop through specified random number iterations:
337 {
338 rand_int = rand() % dperc + 0; //calculate random percentage
339 d = ((dmax/100) * rand_int) + dmin; //calculate value
340 dd = dd + d; //sum
341 n = n + 1;
342 }
343 d = dd / n; // calculate mean from random values
344 }
345 else
346 {
347 d = dmax;
348 }
349
350 eperc = ((emax - emin) / emax) * 100; //calculate parameter range %: depth
351 if (eperc > 0)
352 {
353 n = 0;
354 ee = 0;
355 while ( n < h) //loop through specified random number iterations:
356 {
357 rand_int = rand() % eperc + 0; //calculate random percentage
358 e = ((emax/100) * rand_int) + emin; //calculate value
359 ee = ee + e; //sum
360 n = n + 1;
361 }
362 e = ee / n; // calculate mean from random values
363 }
364 else
365 {
366 e = emax;
367 }
368
369
370 f = (d*(b/pB->Get_Cellsize()))/((c*e)*sin(a)); //Wetness index calculation
371
372 if (f < 1) //calculate wetness index grid
373 pF->Set_Value(x, y, f);
374 else if (f >= 1)
375 pF->Set_Value(x, y, 1);
376 else
377 pF->Set_NoData(x, y);
378
379
380 if (pG) //calculate optional classified grid
381 {
382 if (f <= 0.1)
383 pG->Set_Value(x, y, WI_LOW_MOISTURE);
384 else if( f <= 1 )// if ((f > 0.1) && (f <= 1))
385 pG->Set_Value(x, y, WI_PARTIALLY_WET);
386 else //if (f > 1)
387 pG->Set_Value(x, y, WI_SATURATION_ZONE);
388 }
389 }
390 }
391 }
392
393
394 //-----------------------------------------------------
395
396 CSG_Parameters P;
397
398 if( DataObject_Get_Parameters(pG, P) && P("COLORS_TYPE") && P("LUT") )
399 {
400 int CR_Colors[WI_COUNT] =
401 {
402 SG_GET_RGB(255, 255, 255), // WI_NODATA
403 SG_GET_RGB(255, 255, 0), // WI_LOW_MOISTURE
404 SG_GET_RGB(0, 255, 0), // WI_PARTIALLY_WET
405 SG_GET_RGB(0, 0, 255), // WI_SATURATION_ZONE
406 };
407
408 //-------------------------------------------------
409 CSG_Strings Name, Desc;
410
411 Name += _TL("No data"); Desc += _TL("");
412 Name += _TL("Low moisture (WI <= 0.1"); Desc += _TL("");
413 Name += _TL("Partially wet (WI = 0.1 - 1)"); Desc += _TL("");
414 Name += _TL("Saturation zone (WI > 1)"); Desc += _TL("");
415
416 //-------------------------------------------------
417 CSG_Table *pTable = P("LUT")->asTable();
418
419 pTable->Del_Records();
420
421 for(int i=0; i<WI_COUNT; i++)
422 {
423 CSG_Table_Record *pRecord = pTable->Add_Record();
424
425 pRecord->Set_Value(0, CR_Colors[i]);
426 pRecord->Set_Value(1, Name[i].c_str());
427 pRecord->Set_Value(2, Desc[i].c_str());
428 pRecord->Set_Value(3, i);
429 pRecord->Set_Value(4, i);
430 }
431
432 P("COLORS_TYPE")->Set_Value(1); // Color Classification Type: Lookup Table
433
434 DataObject_Set_Parameters(pG, P);
435 }
436
437 return( true );
438
439 }
440
441
442 ///////////////////////////////////////////////////////////
443 // //
444 // //
445 // //
446 ///////////////////////////////////////////////////////////
447
448 //---------------------------------------------------------
449