rosetta-soil 0.3.2

Installation

pip install rosetta-soil

Quickstart

>>> from rosetta import rosetta

# required ordering for data records:
# [sa (%), si (%), cl (%), bd (g/cm3), th33, th1500]
# sa, si, and cl are required; others optional

>>> data = [
        [30,30,40,1.5,0.3,0.1],
        [20,60,20],
        [55,25,20,1.1]
    ]

>>> mean, stdev, codes = rosetta(3, data)

>>> with np.printoptions(precision=4, suppress=True):
...     print(mean)

[[ 0.0687  0.3839  0.0037  1.5123 10.4522  1.0203  0.2242]
 [ 0.0899  0.4301  0.0038  1.4993 15.8995  0.8909  0.1726]
 [ 0.0913  0.485   0.0097  1.4172 84.7834  2.9015 -0.3463]]

>>> with np.printoptions(precision=4, suppress=True):
...  print(stdev)

[[ 0.0136  0.015   0.0012  0.1176  4.816   0.4968  1.4017]
 [ 0.0067  0.0088  0.0007  0.0454  3.1789  0.449   1.6119]
 [ 0.0128  0.0131  0.0022  0.0577 27.0357  1.8203  1.2248]]

>>> print(codes)

[5 2 3]

Background

The Rosetta pedotransfer functions predict seven parameters for the van Genuchten model of unsaturated soil hydraulic properties

• theta_r : residual volumetric water content

• theta_s : saturated volumetric water content

• alpha : retention shape parameter (1/cm)

• n : retention shape parameter

• ksat : saturated hydraulic conductivity (cm/d)

• k0 : unsaturated conductivity matching point (cm/d)

• L : unsaturated conductivity exponent

Rosetta provides four models for predicting the seven parameters from soil characterization data. The models differ in the required input data

where

• sa, si, cl are percentages of sand, silt and clay

• BD is soil bulk density (g/cm3)

• th33 is the soil volumetric water content at 33 kPa

• th1500 is the soil volumetric water content at 1500 kPa

Three versions of Rosetta are available for predicting theta_r, theta_s alpha, n, and ksat. The versions effectively represent three alternative calibrations or trainings of the four Rosetta models. The references that should be cited when using Rosetta versions 1, 2, and 3 are, respectively:

[1] Schaap, M.G., Leij, F.J., and Van Genuchten, M.T. 2001. ROSETTA: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. Journal of Hydrology 251(3-4): 163-176. doi: 10.1016/S0022-1694(01)00466-8

[2] Schaap, M.G., A. Nemes, and M.T. van Genuchten. 2004. Comparison of Models for Indirect Estimation of Water Retention and Available Water in Surface Soils. Vadose Zone Journal 3(4): 1455-1463. doi: 10.2136/vzj2004.1455

[3] Zhang, Y. and Schaap, M.G. 2017. Weighted recalibration of the Rosetta pedotransfer model with improved estimates of hydraulic parameter distributions and summary statistics (Rosetta3). Journal of Hydrology 547: 39-53. doi: 10.1016/j.jhydrol.2017.01.004

The parameters k0 and L are predicted with a neural network model from [1] that uses as inputs values for the retention parameters theta_r, theta_s, alpha, and n. In the current implementation, the same model is to predict k0 and L regardless of the version of Rosetta used to estimate the retention paramters.

Usage

Soil data should be orgainzed as a collection of data records,

>>> data = [
...     [30, 30, 40, 1.5, 0.3, 0.1],
...     [20, 60, 20],
...     [55, 25, 20, 1.1]
... ]

where each element contains soil data in this order:

[%sand, %silt, %clay, bulk density, th33, th1500]

Sand, silt, and clay are required in each element; the others are optional.

Two functions, rosetta and rosesoil, are provided for predicting soil hydraulic parameters. The two functions have the same arguments and return the same results, but differ in the format of the returned results. The function rosetta returns the predicted parameters in numpy arrays, whereas rosesoil returns RosettaResults object.

rosetta and rosesoil take as arguments the Rosetta version to be used (1, 2, or 3) and the array-like soil data. A third optional argument, estimate_type, specifies the parameter estimate to be returned for alpha, npar, Ksat, and K0:

estimate_type: Literal['linear', 'log', 'geo']
   'linear' (default):
       Estimate is arithmetic mean.
   'log':
       Estimate is mean of log10 transformed parameter value.
   'geo':
       Estimate is geometric mean.

The function rosetta returns a 3-tuple of numpy arrays

>>> from rosetta import rosetta

>>> data = [
...     [30,30,40,1.5,0.3,0.1],
...     [20,60,20],
...     [55,25,20,1.1]
... ]

>>> mean, stdev, codes = rosetta(3, data)

>>> with np.printoptions(precision=4, suppress=True):
...     print(mean)

[[ 0.0687  0.3839  0.0037  1.5123 10.4522  1.0203  0.2242]
 [ 0.0899  0.4301  0.0038  1.4993 15.8995  0.8909  0.1726]
 [ 0.0913  0.485   0.0097  1.4172 84.7834  2.9015 -0.3463]]

>>> with np.printoptions(precision=4, suppress=True):
...  print(stdev)

[[ 0.0136  0.015   0.0012  0.1176  4.816   0.4968  1.4017]
 [ 0.0067  0.0088  0.0007  0.0454  3.1789  0.449   1.6119]
 [ 0.0128  0.0131  0.0022  0.0577 27.0357  1.8203  1.2248]]

>>> print(codes)

[5 2 3]

In the above, mean is a 2D array with shape (nsamples, 7). The ith row holds predicted soil hydraulic parameters for ith entry in data. The columns are:

stdev is 2D array holding the corresponding parameter standard deviations. Note these are standard deviations for the bootstrap resamples, not the standard error of the mean estimate.

codes is a 1D array with the ith entry indicating the Rosetta model and input data used to predict the ith row of mean and stdev.

The function rosesoil works the same way:

>>> from rosetta import rosesoil


>>> data = [
...     [30,30,40,1.5,0.3,0.1],
...     [20,60,20],
...     [55,25,20,1.1]
... ]

>>> result = rosesoil(3, data)

>>> print(result)

RosettaResults([RosettaResult(sand=30.0000, silt=30.0000, clay=40.0000,
rhob=1.5000, th33=0.3000, th1500=0.1000, version=3, estimate_type='linear',
code=5, thr=0.0687, ths=0.3839, alpha=0.0037, npar=1.5123, ksat=10.4522,
k0=1.0203, lpar=0.2242, thr_std=0.0136, ths_std=0.0150, alpha_std=0.0012,
npar_std=0.1176, ksat_std=4.8160, k0_std=0.4968, lpar_std=1.4017),
RosettaResult(sand=20.0000, silt=60.0000, clay=20.0000, rhob=None,
th33=None, th1500=None, version=3, estimate_type='linear', code=2,
thr=0.0899, ths=0.4301, alpha=0.0038, npar=1.4993, ksat=15.8995, k0=0.8909,
lpar=0.1726, thr_std=0.0067, ths_std=0.0088, alpha_std=0.0007,
npar_std=0.0454, ksat_std=3.1789, k0_std=0.4490, lpar_std=1.6119),
RosettaResult(sand=55.0000, silt=25.0000, clay=20.0000, rhob=1.1000,
th33=None, th1500=None, version=3, estimate_type='linear', code=3,
thr=0.0913, ths=0.4850, alpha=0.0097, npar=1.4172, ksat=84.7834, k0=2.9015,
lpar=-0.3463, thr_std=0.0128, ths_std=0.0131, alpha_std=0.0022,
npar_std=0.0577, ksat_std=27.0357, k0_std=1.8203, lpar_std=1.2248)])


To get the results as a list of dictionaries, use `.asdicts()`

>>> print(result.asdicts())

[{'sand': 30.0, 'silt': 30.0, 'clay': 40.0, 'rhob': 1.5, 'th33': 0.3,
'th1500': 0.1, 'version': 3, 'estimate_type': 'linear', 'code': 5, 'thr':
0.06872133198419333, 'ths': 0.3839050853475144, 'alpha':
0.0036914540698334396, 'npar': 1.5122960137476906, 'ksat':
10.45220672557941, 'k0': 1.020283781688175, 'lpar': 0.22415802649553881,
'thr_std': 0.013635805076704855, 'ths_std': 0.014976665455951114,
'alpha_std': 0.0012009965476429596, 'npar_std': 0.11756697611413282,
'ksat_std': 4.816027639984864, 'k0_std': 0.4967836368490979, 'lpar_std':
1.4016987918707104}, {'sand': 20.0, 'silt': 60.0, 'clay': 20.0, 'rhob':
None, 'th33': None, 'th1500': None, 'version': 3, 'estimate_type':
'linear', 'code': 2, 'thr': 0.08994502219206939, 'ths': 0.4301366480210401,
'alpha': 0.0038032951212636345, 'npar': 1.4992978643076722, 'ksat':
15.899549552218094, 'k0': 0.8908630387641998, 'lpar': 0.17264425930013833,
'thr_std': 0.006710949117601514, 'ths_std': 0.008790220586404331,
'alpha_std': 0.0006666891009750212, 'npar_std': 0.04539463666512193,
'ksat_std': 3.178862086739814, 'k0_std': 0.4489771681065423, 'lpar_std':
1.6118830509323745}, {'sand': 55.0, 'silt': 25.0, 'clay': 20.0, 'rhob':
1.1, 'th33': None, 'th1500': None, 'version': 3, 'estimate_type': 'linear',
'code': 3, 'thr': 0.09130753033144609, 'ths': 0.485031958049669, 'alpha':
0.009748718589139088, 'npar': 1.4171964842691995, 'ksat': 84.7833726853498,
'k0': 2.9015059553301388, 'lpar': -0.3463378623447919, 'thr_std':
0.012777797583517441, 'ths_std': 0.013068706563916467, 'alpha_std':
0.0022148332082555515, 'npar_std': 0.05767077407189314, 'ksat_std':
27.035739850886362, 'k0_std': 1.8202787951521058, 'lpar_std':
1.2248193041204392}]

>>> from rosetta import Rosetta

>>> rose33 = Rosetta(rosetta_version=3, model_code=3)

>>> from rosetta import Rosetta
>>> rose33 = Rosetta(rosetta_version=3, model_code=3)
>>> data =[
...     # sa, si, cl, bd
...     [30, 30, 40, 1.5],
...     [55, 25, 20, 1.1]
... ]
>>> retc_boot, ksat_boot = rose33.predict(data)
>>> print(retc_boot.shape)
(1000, 2, 4)
>>> print(ksat_boot.shape)
(1000, 2, 1)

>>> import numpy as np
>>> print(np.mean(retc_boot, axis=0))

[[ 0.11535773  0.41791199 -2.06713851  0.11201021]
[ 0.09130753  0.48503196 -2.02238809  0.15107161]]

>>> from rosetta import UnsaturatedK
>>> unsatk = UnsaturatedK()
>>> retc_params = [
...     # theta_r, theta_s, alpha, npar
...     [0.12, 0.42, 0.008, 1.29],
...     [0.09, 0.49, 0.009, 1.41]
... ]
>>> unsat_boot = unsatk.predict(retc_params)
>>> print(unsat_boot.shape)
(100, 2, 2)
>>> print(np.mean(unsat_boot, axis=0))
[[-0.04057941, -1.03302704],
[ 0.35788353, -0.25048901]]

Acknowledgments

The values for the neural network weights and biases were obtained from Marcel Schaap and Yonggen Zhang at the University of Arizona.