1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Tool Library //
9 // grid_analysis //
10 // //
11 //-------------------------------------------------------//
12 // //
13 // soil_water_capacity.cpp //
14 // //
15 // Copyright (C) 2020 by //
16 // Olaf Conrad //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. SAGA is free software; you //
22 // can redistribute it and/or modify it under the terms //
23 // of the GNU General Public License as published by the //
24 // Free Software Foundation, either version 2 of the //
25 // License, or (at your option) any later version. //
26 // //
27 // SAGA is distributed in the hope that it will be //
28 // useful, but WITHOUT ANY WARRANTY; without even the //
29 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
30 // PARTICULAR PURPOSE. See the GNU General Public //
31 // License for more details. //
32 // //
33 // You should have received a copy of the GNU General //
34 // Public License along with this program; if not, see //
35 // <http://www.gnu.org/licenses/>. //
36 // //
37 //-------------------------------------------------------//
38 // //
39 // e-mail: oconrad@saga-gis.de //
40 // //
41 // contact: Olaf Conrad //
42 // Institute of Geography //
43 // University of Hamburg //
44 // Germany //
45 // //
46 ///////////////////////////////////////////////////////////
47
48 //---------------------------------------------------------
49 #include "soil_water_capacity.h"
50
51
52 ///////////////////////////////////////////////////////////
53 // //
54 // //
55 // //
56 ///////////////////////////////////////////////////////////
57
58 //---------------------------------------------------------
59 double CSoil_Water_Capacity::s_Coefficients[4][12] = { // pedotransfer coefficients developed by Hodnett and Tomasella (2002)
60 { -2.294, 0. , -3.526, 0. , 2.44 , 0. , -0.076, -11.331, 0.019, 0. , 0. , 0. }, // ln(alpha)
61 { 62.986, 0. , 0. , -0.833, -0.529, 0. , 0. , 0.593, 0. , 0.007 , -0.014, 0. }, // ln(n)
62 { 81.799, 0. , 0. , 0.099, 0. , -31.42 , 0.018, 0.451, 0. , 0 , 0. , -5.e-04 }, // theta_s
63 { 22.733, -0.164, 0. , 0. , 0. , 0. , 0.235, -0.831, 0. , 0.0018, 0. , 0.0026 } // theta_r
64 };
65
66
67 ///////////////////////////////////////////////////////////
68 // //
69 ///////////////////////////////////////////////////////////
70
71 //---------------------------------------------------------
CSoil_Water_Capacity(bool bGrids)72 CSoil_Water_Capacity::CSoil_Water_Capacity(bool bGrids)
73 : CSG_Tool_Grid(), m_bGrids(bGrids)
74 {
75 Set_Name (CSG_String::Format(m_bGrids ? "%s (%s)" : "%s", _TL("Soil Water Capacity"), _TL("Grid Collections")));
76
77 Set_Author ("O.Conrad (c) 2020");
78
79 Set_Description (_TW(
80 "This tool derives the soil water capacity for the given "
81 "soil moisture potentials (psi) based on pedo-transfer functions.\n"
82 "Suggested psi values for field capacity estimation range between "
83 "60 hPa (pF=1.8) and 316 hPa (pF=2.5). "
84 "For permanent wilting point estimation take a psi "
85 "value of about 15850 hPa (pF=4.2). "
86 "This tool re-implements the R-script AWCPTF by Hengl as well as the "
87 "regression approach by Toth et al. (2015). "
88 "See Hengl et al. (2017), Woesten & Verzandvoort (2013) "
89 "and Toth et al. (2015) for more details. "
90 ));
91
92 Add_Reference("Hengl, T., Mendes de Jesus, J., Heuvelink, G.B.M., Ruiperez Gonzalez, M., Kilibarda, M., Blagotic, A., et al.", "2017",
93 "SoilGrids250m: Global gridded soil information based on machine learning",
94 "PLoS ONE 12(2): e0169748.",
95 SG_T("https://doi.org/10.1371/journal.pone.0169748"), SG_T("doi:10.1371/journal.pone.0169748")
96 );
97
98 Add_Reference("Hodnett, M.G., Tomasella, J.", "2002",
99 "Marked differences between van Genuchten soil water-retention parameters for temperate and tropical soils: A new water-retention pedo-transfer functions developed for tropical soils",
100 "Geoderma 108(3):155-180.",
101 SG_T("https://doi.org/10.1016/S0016-7061(02)00105-2"), SG_T("doi:10.1016/S0016-7061(02)00105-2")
102 );
103
104 Add_Reference("Toth, B., Weynants, M., Nemes, A., Mako, A., Bilas, G., Toth, G.", "2015",
105 "New generation of hydraulic pedotransfer functions for Europe",
106 "European Journal of Soil Science, 66, 226�238.",
107 SG_T("https://doi.org/10.1111/ejss.12192"), SG_T("doi:10.1111/ejss.12192")
108 );
109
110 Add_Reference("Toth, B., Weynants, M., Pasztor, L, Hengl, T.", "2017",
111 "3D soil hydraulic database of Europe at 250 m resolution",
112 "Hydrological Processes, 31:2662�2666.",
113 SG_T("https://doi.org/10.1002/hyp.11203"), SG_T("doi:10.1002/hyp.11203")
114 );
115
116 Add_Reference("Woesten, J.H.M., Verzandvoort, S.J.E., Leenaars, J.G.B., Hoogland, T., Wesseling, J.G.", "2013",
117 "Soil hydraulic information for river basin studies in semi-arid regions",
118 "Geoderma 195:79-86.",
119 SG_T("https://doi.org/10.1016/j.geoderma.2012.11.021"), SG_T("doi:10.1016/j.geoderma.2012.11.021")
120 );
121
122 Add_Reference(
123 "https://github.com/cran/GSIF/blob/master/R/AWCPTF.R",
124 SG_T("R-Script AWCPTF {GSIF} by J. Leenaars, M. Ruiperez Gonzalez and T. Hengl")
125 );
126
127 //-----------------------------------------------------
128 if( m_bGrids )
129 {
130 Parameters.Add_Grids ("", "SAND" , _TL("Sand" ), _TL("[%]" ), PARAMETER_INPUT );
131 Parameters.Add_Grids ("", "SILT" , _TL("Silt" ), _TL("[%]" ), PARAMETER_INPUT );
132 Parameters.Add_Grids ("", "CLAY" , _TL("Clay" ), _TL("[%]" ), PARAMETER_INPUT );
133 Parameters.Add_Grids ("", "CORG" , _TL("Soil Organic Carbon" ), _TL("[permille]"), PARAMETER_INPUT );
134 Parameters.Add_Grids ("", "BULK" , _TL("Bulk Density" ), _TL("[kg/m^3]" ), PARAMETER_INPUT );
135 Parameters.Add_Grids ("", "CEC" , _TL("Cation Exchange Capacity" ), _TL("[cmol/kg]" ), PARAMETER_INPUT );
136 Parameters.Add_Grids ("", "PH" , _TL("pH" ), _TL("" ), PARAMETER_INPUT );
137
138 Parameters.Add_Grids ("", "FC" , _TL("Field Capacity" ), _TL("" ), PARAMETER_OUTPUT );
139 Parameters.Add_Grids ("", "PWP" , _TL("Permanent Wilting Point" ), _TL("" ), PARAMETER_OUTPUT );
140 Parameters.Add_Grids ("", "THETA_S", _TL("Water Capacity at Saturation"), _TL("" ), PARAMETER_OUTPUT_OPTIONAL );
141 }
142 else
143 {
144 Parameters.Add_Grid_or_Const("", "SAND" , _TL("Sand" ), _TL("[%]" ), 15.0, 0., true, 100., true);
145 Parameters.Add_Grid_or_Const("", "SILT" , _TL("Silt" ), _TL("[%]" ), 37.0, 0., true, 100., true);
146 Parameters.Add_Grid_or_Const("", "CLAY" , _TL("Clay" ), _TL("[%]" ), 48.0, 0., true, 100., true);
147 Parameters.Add_Grid_or_Const("", "CORG" , _TL("Soil Organic Carbon" ), _TL("[permille]"), 15.0, 0., true, 1000., true);
148 Parameters.Add_Grid_or_Const("", "BULK" , _TL("Bulk Density" ), _TL("[kg/m^3]" ), 1350.0, 100., true, 2650., true);
149 Parameters.Add_Grid_or_Const("", "CEC" , _TL("Cation Exchange Capacity" ), _TL("[cmol/kg]" ), 45.0, 0., true );
150 Parameters.Add_Grid_or_Const("", "PH" , _TL("pH" ), _TL("" ), 6.4, 0., true, 14., true);
151
152 Parameters.Add_Grid ("", "FC" , _TL("Field Capacity" ), _TL("" ), PARAMETER_OUTPUT );
153 Parameters.Add_Grid ("", "PWP" , _TL("Permanent Wilting Point" ), _TL("" ), PARAMETER_OUTPUT );
154 Parameters.Add_Grid ("", "THETA_S", _TL("Water Capacity at Saturation"), _TL("" ), PARAMETER_OUTPUT_OPTIONAL );
155 }
156
157 //-----------------------------------------------------
158 Parameters.Add_Choice("",
159 "UNIT" , _TL("Output Unit"),
160 _TL(""),
161 CSG_String::Format("%s|%s",
162 _TL("cubic-meter per cubic-meter"),
163 _TL("percentage of volume")
164 ), 1
165 );
166
167 Parameters.Add_Choice("",
168 "FUNCTION" , _TL("Pedo-Transfer Function"),
169 _TL(""),
170 CSG_String::Format("%s|%s",
171 _TL("Hodnett & Tomasella 2002"),
172 _TL("Toth et al. 2015")
173 )
174 );
175
176 //-----------------------------------------------------
177 Parameters.Add_Double("FC" , "PSI_FC" , _TL("Soil Moisture Potential"), _TL("[hPa]"), 316., 0., true);
178 Parameters.Add_Double("PWP", "PSI_PWP", _TL("Soil Moisture Potential"), _TL("[hPa]"), 15850., 0., true);
179
180 Parameters.Add_Bool("",
181 "ADJUST" , _TL("Adjustments"),
182 _TL("Specifies whether to correct values of textures and bulk density to avoid creating nonsensical values."),
183 true
184 );
185
186 Parameters.Add_Bool("",
187 "USERDEF" , _TL("User Defined Coefficients"),
188 _TL(""),
189 false
190 );
191
192 CSG_Table &a = *Parameters.Add_FixedTable("USERDEF",
193 "COEFFICIENTS" , _TL("User Defined Coefficients"),
194 _TL("")
195 )->asTable();
196
197 a.Destroy();
198 a.Set_Name(_TL("User Defined Coefficients"));
199 a.Add_Field("ln(alpha)", SG_DATATYPE_Double);
200 a.Add_Field("ln(n)" , SG_DATATYPE_Double);
201 a.Add_Field("theta s" , SG_DATATYPE_Double);
202 a.Add_Field("theta r" , SG_DATATYPE_Double);
203 a.Set_Count(12);
204
205 for(int j=0; j<12; j++)
206 {
207 for(int i=0; i<4; i++)
208 {
209 a[j].Set_Value(i, s_Coefficients[i][j]);
210 }
211 }
212 }
213
214
215 ///////////////////////////////////////////////////////////
216 // //
217 ///////////////////////////////////////////////////////////
218
219 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)220 int CSoil_Water_Capacity::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
221 {
222 if( pParameter->Cmp_Identifier("FUNCTION") )
223 {
224 pParameters->Set_Enabled("SAND" , pParameter->asInt() == 0);
225 pParameters->Set_Enabled("BULK" , pParameter->asInt() == 0);
226 pParameters->Set_Enabled("CEC" , pParameter->asInt() == 0);
227 pParameters->Set_Enabled("PH" , pParameter->asInt() == 0);
228 pParameters->Set_Enabled("THETA_S" , pParameter->asInt() == 0);
229 pParameters->Set_Enabled("PSI_FC" , pParameter->asInt() == 0);
230 pParameters->Set_Enabled("PSI_PWP" , pParameter->asInt() == 0);
231 pParameters->Set_Enabled("ADJUST" , pParameter->asInt() == 0);
232 pParameters->Set_Enabled("USERDEF" , pParameter->asInt() == 0);
233 }
234
235 if( pParameter->Cmp_Identifier("USERDEF") )
236 {
237 pParameters->Set_Enabled("COEFFICIENTS", pParameter->asBool());
238 }
239
240 return( CSG_Tool_Grid::On_Parameters_Enable(pParameters, pParameter) );
241 }
242
243
244 ///////////////////////////////////////////////////////////
245 // //
246 ///////////////////////////////////////////////////////////
247
248 //---------------------------------------------------------
On_Execute(void)249 bool CSoil_Water_Capacity::On_Execute(void)
250 {
251 switch( Parameters("FUNCTION")->asInt() )
252 {
253 default: return( Get_HodnettTomasella() );
254 case 1: return( Get_Toth () );
255 }
256 }
257
258
259 ///////////////////////////////////////////////////////////
260 // //
261 ///////////////////////////////////////////////////////////
262
263 //---------------------------------------------------------
Get_HodnettTomasella(void)264 bool CSoil_Water_Capacity::Get_HodnettTomasella(void)
265 {
266 double Scale = Parameters("UNIT")->asInt() == 1 ? 100. : 1.;
267
268 double psi_FC = Parameters("PSI_FC" )->asDouble() / 10.; // [hPa] -> [kPa]
269 double psi_PWP = Parameters("PSI_PWP")->asDouble() / 10.; // [hPa] -> [kPa]
270
271 bool bAdjust = Parameters("ADJUST")->asBool();
272
273 //-----------------------------------------------------
274 m_Coefficients.Create(12, 4);
275
276 switch( Parameters("FUNCTION")->asInt() )
277 {
278 default: {
279 for(int j=0; j<12; j++)
280 {
281 for(int i=0; i<4; i++)
282 {
283 m_Coefficients[i][j] = s_Coefficients[i][j];
284 }
285 }
286 break; }
287
288 case 1: {
289 CSG_Table &a = *Parameters("COEFFICIENTS")->asTable();
290
291 if( a.Get_Count() != 12 )
292 {
293 Error_Set(_TL("User defined table needs to provide not less than 12 rows."));
294
295 return( false );
296 }
297
298 for(int j=0; j<12; j++)
299 {
300 for(int i=0; i<4; i++)
301 {
302 m_Coefficients[i][j] = a[j][i];
303 }
304 }
305 break; }
306 }
307
308 //-----------------------------------------------------
309 if( m_bGrids )
310 {
311 CSG_Table Layers;
312
313 #define CHECK_LAYERS(pGrids) if( pGrids && pGrids->Get_NZ() && (!Layers.Get_Count() || Layers.Get_Count() > pGrids->Get_NZ()) ) { Layers.Create(pGrids->Get_Attributes()); }
314 CSG_Grids *pSand = Parameters("SAND" )->asGrids(); CHECK_LAYERS(pSand);
315 CSG_Grids *pSilt = Parameters("SILT" )->asGrids(); CHECK_LAYERS(pSilt);
316 CSG_Grids *pClay = Parameters("CLAY" )->asGrids(); CHECK_LAYERS(pClay);
317 CSG_Grids *pBulk = Parameters("BULK" )->asGrids(); CHECK_LAYERS(pBulk);
318 CSG_Grids *pCorg = Parameters("CORG" )->asGrids(); CHECK_LAYERS(pCorg);
319 CSG_Grids *pCEC = Parameters("CEC" )->asGrids(); CHECK_LAYERS(pCEC );
320 CSG_Grids *ppH = Parameters("PH" )->asGrids(); CHECK_LAYERS(ppH );
321 #undef CHECK_LAYERS
322
323 if( !Layers.Get_Count() )
324 {
325 Error_Set(_TL("Each grid collection must provide at least one level."));
326
327 return( false );
328 }
329
330 #define CREATE_LAYERS(pGrids, Name) if( pGrids ) { pGrids->Create(Get_System(), Layers, 0, SG_DATATYPE_Short, true); pGrids->Set_Scaling(Scale / 1000.); pGrids->Set_Name(Name); }
331 CSG_Grids *pFC = Parameters("FC" )->asGrids(); CREATE_LAYERS(pFC , "FC" );
332 CSG_Grids *pPWP = Parameters("PWP" )->asGrids(); CREATE_LAYERS(pPWP, "PWP" );
333 CSG_Grids *pT_s = Parameters("THETA_S")->asGrids(); CREATE_LAYERS(pT_s, "Theta s");
334 #undef CREATE_LAYERS
335
336 //-------------------------------------------------
337 for(int z=0; z<Layers.Get_Count() && Process_Get_Okay(); z++)
338 {
339 Process_Set_Text("%s [%d/%d]", _TL("processing"), z + 1, Layers.Get_Count());
340
341 for(int y=0; y<Get_NY() && Set_Progress(y); y++)
342 {
343 #pragma omp parallel for
344 for(int x=0; x<Get_NX(); x++)
345 {
346 #define SET_NODATA(x, y, z) {\
347 if( pFC ) pFC ->Set_NoData(x, y, z);\
348 if( pPWP ) pPWP->Set_NoData(x, y, z);\
349 if( pT_s ) pT_s->Set_NoData(x, y, z);\
350 continue;\
351 }
352
353 if( (pSand && pSand->is_NoData(x, y, z))
354 || (pSilt && pSilt->is_NoData(x, y, z))
355 || (pClay && pClay->is_NoData(x, y, z))
356 || (pBulk && pBulk->is_NoData(x, y, z))
357 || (pCorg && pCorg->is_NoData(x, y, z))
358 || (pCEC && pCEC ->is_NoData(x, y, z))
359 || (ppH && ppH ->is_NoData(x, y, z)) )
360 {
361 SET_NODATA(x, y, z);
362 }
363
364 //-------------------------------------
365 double Sand = pSand->asDouble(x, y, z);
366 double Silt = pSilt->asDouble(x, y, z);
367 double Clay = pClay->asDouble(x, y, z);
368 double Bulk = pBulk->asDouble(x, y, z);
369 double Corg = pCorg->asDouble(x, y, z);
370 double CEC = pCEC ->asDouble(x, y, z);
371 double pH = ppH ->asDouble(x, y, z);
372
373 //-------------------------------------
374 double alpha, n, theta_s, theta_r;
375
376 if( !Get_HodnettTomasella(alpha, n, theta_s, theta_r, Sand, Silt, Clay, Bulk, Corg, CEC, pH, bAdjust) )
377 {
378 SET_NODATA(x, y, z);
379 }
380
381 if( pFC ) pFC ->Set_Value(x, y, z, Scale * van_Genuchten(psi_FC , alpha, n, theta_s, theta_r));
382 if( pPWP ) pPWP->Set_Value(x, y, z, Scale * van_Genuchten(psi_PWP, alpha, n, theta_s, theta_r));
383 if( pT_s ) pT_s->Set_Value(x, y, z, Scale * theta_s);
384
385 #undef SET_NODATA
386 }
387 }
388 }
389 }
390
391 //-----------------------------------------------------
392 else
393 {
394 CSG_Grid *pSand = Parameters("SAND" )->asGrid(); double cSand = Parameters("SAND")->asDouble();
395 CSG_Grid *pSilt = Parameters("SILT" )->asGrid(); double cSilt = Parameters("SILT")->asDouble();
396 CSG_Grid *pClay = Parameters("CLAY" )->asGrid(); double cClay = Parameters("CLAY")->asDouble();
397 CSG_Grid *pBulk = Parameters("BULK" )->asGrid(); double cBulk = Parameters("BULK")->asDouble();
398 CSG_Grid *pCorg = Parameters("CORG" )->asGrid(); double cCorg = Parameters("CORG")->asDouble();
399 CSG_Grid *pCEC = Parameters("CEC" )->asGrid(); double cCEC = Parameters("CEC" )->asDouble();
400 CSG_Grid *ppH = Parameters("PH" )->asGrid(); double cpH = Parameters("PH" )->asDouble();
401
402 CSG_Grid *pFC = Parameters("FC" )->asGrid();
403 CSG_Grid *pPWP = Parameters("PWP" )->asGrid();
404 CSG_Grid *pT_s = Parameters("THETA_S")->asGrid();
405
406 //-------------------------------------------------
407 for(int y=0; y<Get_NY() && Set_Progress(y); y++)
408 {
409 #pragma omp parallel for
410 for(int x=0; x<Get_NX(); x++)
411 {
412 #define SET_NODATA(x, y) {\
413 if( pFC ) pFC ->Set_NoData(x, y);\
414 if( pPWP ) pPWP->Set_NoData(x, y);\
415 if( pT_s ) pT_s->Set_NoData(x, y);\
416 continue;\
417 }
418
419 if( (pSand && pSand->is_NoData(x, y))
420 || (pSilt && pSilt->is_NoData(x, y))
421 || (pClay && pClay->is_NoData(x, y))
422 || (pBulk && pBulk->is_NoData(x, y))
423 || (pCorg && pCorg->is_NoData(x, y))
424 || (pCEC && pCEC ->is_NoData(x, y))
425 || (ppH && ppH ->is_NoData(x, y)) )
426 {
427 SET_NODATA(x, y);
428 }
429
430 //-----------------------------------------
431 double Sand = pSand ? pSand->asDouble(x, y) : cSand;
432 double Silt = pSilt ? pSilt->asDouble(x, y) : cSilt;
433 double Clay = pClay ? pClay->asDouble(x, y) : cClay;
434 double Bulk = pBulk ? pBulk->asDouble(x, y) : cBulk;
435 double Corg = pCorg ? pCorg->asDouble(x, y) : cCorg;
436 double CEC = pCEC ? pCEC ->asDouble(x, y) : cCEC ;
437 double pH = ppH ? ppH ->asDouble(x, y) : cpH ;
438
439 //-----------------------------------------
440 double alpha, n, theta_s, theta_r;
441
442 if( !Get_HodnettTomasella(alpha, n, theta_s, theta_r, Sand, Silt, Clay, Bulk, Corg, CEC, pH, bAdjust) )
443 {
444 SET_NODATA(x, y);
445 }
446
447 if( pFC ) pFC ->Set_Value(x, y, Scale * van_Genuchten(psi_FC , alpha, n, theta_s, theta_r));
448 if( pPWP ) pPWP->Set_Value(x, y, Scale * van_Genuchten(psi_PWP, alpha, n, theta_s, theta_r));
449 if( pT_s ) pT_s->Set_Value(x, y, Scale * theta_s);
450
451 #undef SET_NODATA
452 }
453 }
454 }
455
456 //-----------------------------------------------------
457 return( true );
458 }
459
460
461 ///////////////////////////////////////////////////////////
462 // //
463 ///////////////////////////////////////////////////////////
464
465 //---------------------------------------------------------
van_Genuchten(double psi,double alpha,double n,double theta_s,double theta_r)466 inline double CSoil_Water_Capacity::van_Genuchten(double psi, double alpha, double n, double theta_s, double theta_r)
467 {
468 return( theta_r + (theta_s - theta_r) / pow(1. + pow(alpha * psi, n), 1. - 1. / n) );
469 }
470
471 //---------------------------------------------------------
Get_HodnettTomasella(double & alpha,double & n,double & theta_s,double & theta_r,double Sand,double Silt,double Clay,double Bulk,double Corg,double CEC,double pH,bool bAdjust)472 inline bool CSoil_Water_Capacity::Get_HodnettTomasella(double &alpha, double &n, double &theta_s, double &theta_r, double Sand, double Silt, double Clay, double Bulk, double Corg, double CEC, double pH, bool bAdjust)
473 {
474 if( bAdjust ) // standardize sand, silt, clay and check bulk density...
475 {
476 double Sum = (Clay + Silt + Sand) / 100.;
477
478 if( Sum <= 0. )
479 {
480 return( false );
481 }
482
483 Clay /= Sum; Silt /= Sum; Sand /= Sum;
484
485 if( Bulk < 100. ) { Bulk = 100.; } else if( Bulk > 2650. ) { Bulk = 2650.; } // density of quartz [kg/m^3]
486 }
487
488 //-----------------------------------------------------
489 double X[4];
490
491 for(int i=0; i<4; i++)
492 {
493 X[i] = // Note: Formula available from [http://www.sciencedirect.com/science/article/pii/S001670611200417X]
494 ( m_Coefficients[i][ 0]
495 + m_Coefficients[i][ 1] * Sand
496 + m_Coefficients[i][ 2] * Silt
497 + m_Coefficients[i][ 3] * Clay
498 + m_Coefficients[i][ 4] * Corg / 10.
499 + m_Coefficients[i][ 5] * Bulk / 1000.
500 + m_Coefficients[i][ 6] * CEC
501 + m_Coefficients[i][ 7] * pH
502 + m_Coefficients[i][ 8] * Silt*Silt
503 + m_Coefficients[i][ 9] * Clay*Clay
504 + m_Coefficients[i][10] * Sand*Silt
505 + m_Coefficients[i][11] * Sand*Clay
506 ) / 100.;
507 }
508
509 alpha = exp(X[0]);
510 n = exp(X[1]);
511
512 theta_s = X[2] > 100. ? 100. : X[2];
513 theta_r = X[3] < 0. ? 0. : X[3];
514
515 return( true );
516 }
517
518
519 ///////////////////////////////////////////////////////////
520 // //
521 ///////////////////////////////////////////////////////////
522
523 //---------------------------------------------------------
Get_Toth(void)524 bool CSoil_Water_Capacity::Get_Toth(void)
525 {
526 double Scale = Parameters("UNIT")->asInt() == 1 ? 100. : 1.;
527
528 //-----------------------------------------------------
529 if( m_bGrids )
530 {
531 CSG_Table Layers;
532
533 #define CHECK_LAYERS(pGrids) if( pGrids && pGrids->Get_NZ() && (!Layers.Get_Count() || Layers.Get_Count() > pGrids->Get_NZ()) ) { Layers.Create(pGrids->Get_Attributes()); }
534 CSG_Grids *pSilt = Parameters("SILT" )->asGrids(); CHECK_LAYERS(pSilt);
535 CSG_Grids *pClay = Parameters("CLAY" )->asGrids(); CHECK_LAYERS(pClay);
536 CSG_Grids *pCorg = Parameters("CORG" )->asGrids(); CHECK_LAYERS(pCorg);
537 #undef CHECK_LAYERS
538
539 if( !Layers.Get_Count() )
540 {
541 Error_Set(_TL("Each grid collection must provide at least one z level."));
542
543 return( false );
544 }
545
546 #define CREATE_LAYERS(pGrids, Name) if( pGrids ) { pGrids->Create(Get_System(), Layers, 0, SG_DATATYPE_Short, true); pGrids->Set_Scaling(Scale / 1000.); pGrids->Set_Name(Name); }
547 CSG_Grids *pFC = Parameters("FC" )->asGrids(); CREATE_LAYERS(pFC , "FC" );
548 CSG_Grids *pPWP = Parameters("PWP")->asGrids(); CREATE_LAYERS(pPWP, "PWP");
549 #undef CREATE_LAYERS
550
551 //-------------------------------------------------
552 for(int z=0; z<Layers.Get_Count() && Process_Get_Okay(); z++)
553 {
554 Process_Set_Text("%s [%d/%d]", _TL("processing"), z + 1, Layers.Get_Count());
555
556 for(int y=0; y<Get_NY() && Set_Progress(y); y++)
557 {
558 #pragma omp parallel for
559 for(int x=0; x<Get_NX(); x++)
560 {
561 #define SET_NODATA(x, y, z) {\
562 if( pFC ) pFC ->Set_NoData(x, y, z);\
563 if( pPWP ) pPWP->Set_NoData(x, y, z);\
564 continue;\
565 }
566
567 if( (pSilt && pSilt->is_NoData(x, y, z))
568 || (pClay && pClay->is_NoData(x, y, z))
569 || (pCorg && pCorg->is_NoData(x, y, z)) )
570 {
571 SET_NODATA(x, y, z);
572 }
573
574 double Silt = pSilt->asDouble(x, y, z);
575 double Clay = pClay->asDouble(x, y, z);
576 double Corg = pCorg->asDouble(x, y, z);
577
578 double FC, PWP;
579
580 if( !Get_Toth(FC, PWP, Silt, Clay, Corg) )
581 {
582 SET_NODATA(x, y, z);
583 }
584
585 if( pFC ) pFC ->Set_Value(x, y, z, Scale * FC );
586 if( pPWP ) pPWP->Set_Value(x, y, z, Scale * PWP);
587
588 #undef SET_NODATA
589 }
590 }
591 }
592 }
593
594 //-----------------------------------------------------
595 else
596 {
597 CSG_Grid *pSilt = Parameters("SILT")->asGrid(); double cSilt = Parameters("SILT")->asDouble();
598 CSG_Grid *pClay = Parameters("CLAY")->asGrid(); double cClay = Parameters("CLAY")->asDouble();
599 CSG_Grid *pCorg = Parameters("CORG")->asGrid(); double cCorg = Parameters("CORG")->asDouble();
600
601 CSG_Grid *pFC = Parameters("FC" )->asGrid();
602 CSG_Grid *pPWP = Parameters("PWP")->asGrid();
603
604 //-------------------------------------------------
605 for(int y=0; y<Get_NY() && Set_Progress(y); y++)
606 {
607 #pragma omp parallel for
608 for(int x=0; x<Get_NX(); x++)
609 {
610 #define SET_NODATA(x, y) {\
611 if( pFC ) pFC ->Set_NoData(x, y);\
612 if( pPWP ) pPWP->Set_NoData(x, y);\
613 continue;\
614 }
615
616 if( (pSilt && pSilt->is_NoData(x, y))
617 || (pClay && pClay->is_NoData(x, y))
618 || (pCorg && pCorg->is_NoData(x, y)) )
619 {
620 SET_NODATA(x, y);
621 }
622
623 //-----------------------------------------
624 double Silt = pSilt ? pSilt->asDouble(x, y) : cSilt;
625 double Clay = pClay ? pClay->asDouble(x, y) : cClay;
626 double Corg = pCorg ? pCorg->asDouble(x, y) : cCorg;
627
628 double FC, PWP;
629
630 if( !Get_Toth(FC, PWP, Silt, Clay, Corg) )
631 {
632 SET_NODATA(x, y);
633 }
634
635 if( pFC ) pFC ->Set_Value(x, y, Scale * FC );
636 if( pPWP ) pPWP->Set_Value(x, y, Scale * PWP);
637
638 #undef SET_NODATA
639 }
640 }
641 }
642
643 //-----------------------------------------------------
644 return( true );
645 }
646
647 //---------------------------------------------------------
Get_Toth(double & FC,double & PWP,double Silt,double Clay,double Corg)648 inline bool CSoil_Water_Capacity::Get_Toth(double &FC, double &PWP, double Silt, double Clay, double Corg)
649 {
650 Corg = 1. / (1. + Corg);
651
652 FC = 0.2449 - 0.1887 * Corg + 0.004527 * Clay + 0.001535 * Silt + 0.001442 * Silt * Corg - 0.00005110 * Silt * Clay + 0.0008676 * Clay * Corg;
653 PWP = 0.09878 + 0.002127 * Clay - 0.0008366 * Silt - 0.07670 * Corg + 0.00003853 * Silt * Clay + 0.002330 * Clay * Corg + 0.0009498 * Silt * Corg;
654
655 return( true );
656 }
657
658
659 ///////////////////////////////////////////////////////////
660 // //
661 // //
662 // //
663 ///////////////////////////////////////////////////////////
664
665 //---------------------------------------------------------
666 //# R-Script by Tomislav Hengl version 0.5-5: https://github.com/cran/GSIF/blob/master/R/AWCPTF.R
667 //
668 //# Note: Formula available from [http://www.sciencedirect.com/science/article/pii/S001670611200417X]
669 //
670 //AWCPTF <- function(SNDPPT, SLTPPT, CLYPPT, ORCDRC, BLD=1400, CEC, PHIHOX, h1=-10, h2=-20, h3=-31.6, pwp=-1585, PTF.coef, fix.values=TRUE, print.coef=TRUE){
671 // ## pedotransfer coefficients developed by Hodnett and Tomasella (2002)
672 // if(missing(PTF.coef)){
673 // PTF.coef <- data.frame(
674 // lnAlfa = c(-2.294, 0, -3.526, 0, 2.44, 0, -0.076, -11.331, 0.019, 0, 0, 0),
675 // lnN = c(62.986, 0, 0, -0.833, -0.529, 0, 0, 0.593, 0, 0.007, -0.014, 0),
676 // tetaS = c(81.799, 0, 0, 0.099, 0, -31.42, 0.018, 0.451, 0, 0, 0, -5e-04),
677 // tetaR = c(22.733, -0.164, 0, 0, 0, 0, 0.235, -0.831, 0, 0.0018, 0, 0.0026)
678 // )
679 // }
680 // ## standardize sand silt clay:
681 // if(fix.values){
682 // sum.tex <- CLYPPT+SLTPPT+SNDPPT
683 // CLYPPT <- CLYPPT/(sum.tex)*100
684 // SLTPPT <- SLTPPT/(sum.tex)*100
685 // SNDPPT <- SNDPPT/(sum.tex)*100
686 // BLD[BLD<100] <- 100
687 // BLD[BLD>2650] <- 2650 ## weight of quartz
688 // }
689 // ## rows:
690 // clm <- data.frame(SNDPPT, SLTPPT, CLYPPT, ORCDRC/10, BLD*0.001, CEC, PHIHOX, SLTPPT^2, CLYPPT^2, SNDPPT*SLTPPT, SNDPPT*CLYPPT)
691 // alfa <- apply(clm, 1, function(x){ exp((PTF.coef$lnAlfa[1] + sum(PTF.coef$lnAlfa[-1] * x))/100) })
692 // N <- apply(clm, 1, function(x){ exp((PTF.coef$lnN[1] + sum(PTF.coef$lnN[-1] * x))/100) })
693 // tetaS <- apply(clm, 1, function(x){ (PTF.coef$tetaS[1] + sum(PTF.coef$tetaS[-1] * x))/100 })
694 // tetaR <- apply(clm, 1, function(x){ (PTF.coef$tetaR[1] + sum(PTF.coef$tetaR[-1] * x))/100 })
695 // ## change negative of tetaR to 0
696 // tetaR[tetaR < 0] <- 0
697 // tetaS[tetaS > 100] <- 100
698 // m <- 1-1/N
699 // tetah1 <- tetaR + (tetaS-tetaR)/((1+(alfa*-1*h1)^N))^m
700 // tetah2 <- tetaR + (tetaS-tetaR)/((1+(alfa*-1*h2)^N))^m
701 // tetah3 <- tetaR + (tetaS-tetaR)/((1+(alfa*-1*h3)^N))^m
702 // WWP <- tetaR + (tetaS-tetaR)/((1+(alfa*-1*pwp)^N))^m
703 // if(fix.values){
704 // ## if any of the tetah values is smaller than WWP, then replace:
705 // sel <- which(WWP > tetah1 | WWP > tetah2 | WWP > tetah3)
706 // if(length(sel)>0){
707 // WWP[sel] <- apply(data.frame(tetah1[sel], tetah2[sel], tetah3[sel]), 1, function(x){min(x, na.rm=TRUE)})
708 // warning(paste("Wilting point capacity for", length(sel), "points higher than h1, h2 and/or h3"))
709 // }
710 // }
711 // AWCh1 <- tetah1 - WWP
712 // AWCh2 <- tetah2 - WWP
713 // AWCh3 <- tetah3 - WWP
714 // out <- data.frame(AWCh1=signif(AWCh1,3), AWCh2=signif(AWCh2,3), AWCh3=signif(AWCh3,3), WWP=signif(WWP,3), tetaS=signif(tetaS,3))
715 // if(print.coef==TRUE){
716 // attr(out, "coef") <- as.list(PTF.coef)
717 // attr(out, "PTF.names") <- list(variable=c("ai1", "sand", "silt", "clay", "oc", "bd", "cec", "ph", "silt^2", "clay^2", "sand*silt", "sand*clay"))
718 // }
719 // return(out)
720 //}
721
722
723 ///////////////////////////////////////////////////////////
724 // //
725 // //
726 // //
727 ///////////////////////////////////////////////////////////
728
729 //---------------------------------------------------------
730