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