Forum: open-discussion

RE: Soil properties DB [ Reply ]
By: Julien Moeys on 2010-10-06 21:07

Well, it seems the attachement has not been attached. Here is the code:

# This code is provided as a quick example, and without garantee
# concerning bugs and units convertions!
require("soilwaterfun")
require("soilwaterptf")
require("soiltexture") # on r-forge or CRAN

# 0. Defines "a-priori" organic carbon, bulk density and topsoil /
# sub-soil parameter
om <- 2 # %
bulkD <- 1.2 # kg.dm-3
topSoil <- 1 # 0 for subsoil

# 1. Calculates classes centroid for the USDA texture
# triangle: (soiltexture package)
soilProp <- TT.get("USDA.TT")
soilProp <- do.call(
"rbind",
lapply(
X = soilProp[["tt.polygons"]],
FUN = function(X){
pol <- soilProp[["tt.points"]][ X[["points"]], ]
#
centr <- as.numeric( TT.polygon.centroids( # Not documented
pol.x = pol[,"CLAY"],
pol.y = pol[,"SAND"]
) ) #
#
c("CLAY"=centr[1],"SILT"=1-centr[1]-centr[2],"SAND"=centr[2]) * 100
} #
) #
) #

# 2. Calculates the properties of each classes (soilwaterptf)
soilProp <- data.frame(
"clay" = soilProp[,"CLAY"],
"silt" = soilProp[,"SILT"],
"bulkD" = bulkD,
"om" = om,
"topSoil" = topSoil
) #

soilProp <- data.frame(
soilProp,
as.data.frame( ptf.wosten( soilprop = soilProp ) ),
"thetaR" = 0.010
) #
# NB: thetaR estimation has not been implemented yet (from classes
# PTFs). But let fix it to 0.010 (not so good for coarse soils)

# 3. Plot the result (soilwaterfun package)
par("ask"=TRUE,mfrow=c(2,1))
for( textClass in rownames(soilProp) )
{ #
curve(
fun.vangenuchten.theta.h(
h = -x, # [m]
alpha = soilProp[textClass,"alpha"], # [m-1]
n = soilProp[textClass,"n"],
thetaS = soilProp[textClass,"thetaS"],
thetaR = soilProp[textClass,"thetaR"]
), #
xlim = c(0.01,158), # [m]
ylim = c(0,0.7),
col = "red",
log = "x",
xlab = "-h",
ylab = expression( theta ),
main = paste("Water retention for USDA texture class",textClass)
) #
#
curve(
fun.mualem.vangenuchten.K.h(
h = -x,
Ks = soilProp[textClass,"kSat"], # [mm.h-1]
alpha = soilProp[textClass,"alpha"], # [m-1]
n = soilProp[textClass,"n"]
), #
xlim = c(0.01,158),
ylim = c(1e-7,max(soilProp[,"kSat"])),
col = "red",
log = "xy",
xlab = "-h",
ylab = expression( K ~ "[mm.h"^-1 * "]" ),
main = paste("Water retention for USDA texture class",textClass)
) #
} #

# Now the same exercise, but this time using a table containing
# the class values of US soils (After the Rosetta software)
# Schaap MG, Leij FJ, van Genuchten MT, 2001. ROSETTA:
# a computer program for estimating soil hydraulic parameters
# with hierarchical pedotransfer functions. Journal of Hydrology,
# 251:163-176 (issue 3-4).
# NB: the table is undocumented at the moment in this package.

soilProp2 <- read.csv(
"http://r-forge.r-project.org/scm/viewvc.php/*checkout*/pkg/soilwaterptf/inst/rosetta_class_ptf_values.csv?revision=11&root=soilwater",
stringsAsFactors=FALSE
) #
soilProp2b <- t( soilProp2[-1,-(1:2)] )
mode(soilProp2b) <- "numeric"
colnames(soilProp2b) <- soilProp2[-1,2]

par("ask"=TRUE,mfrow=c(2,1)) # Notice the differences in thetaR
# # with Wosten 1999
for( textClass in rownames(soilProp2b) )
{ #
curve(
fun.vangenuchten.theta.h(
h = -x, # [m]
alpha = (10^soilProp2b[textClass,"log.alpha"]) * 100, # [m-1]
n = 10^soilProp2b[textClass,"log.n"],
thetaS = soilProp2b[textClass,"theta.s"],
thetaR = soilProp2b[textClass,"theta.r"]
), #
xlim = c(0.01,158), # [m]
ylim = c(0,0.7),
col = "red",
log = "x",
xlab = "-h",
ylab = expression( theta ),
main = paste("Water retention for USDA texture class",textClass)
) #
#
curve(
fun.mualem.vangenuchten.K.h(
h = -x,
Ks = (10^soilProp2b[textClass,"log.Ks"])*(10/24), # [mm.h-1]
alpha = (10^soilProp2b[textClass,"log.alpha"]) * 100, # [m-1]
n = 10^soilProp2b[textClass,"log.n"],
b = soilProp2b[textClass,"L"]
), #
xlim = c(0.01,158),
ylim = c(1e-7,max(soilProp[,"kSat"])),
col = "red",
log = "xy",
xlab = "-h",
ylab = expression( K ~ "[mm.h"^-1 * "]" ),
main = paste("Water retention for USDA texture class",textClass)
) #
} #

RE: Soil properties DB [ Reply ]
By: Julien Moeys on 2010-10-06 21:05

[forum:3459]

Hi Mikhail. Interesting question.

It is risky to define any "typical" soils in such a package. What is a typical soil? For which region of the world? (etc.)

But both Wosten et al. 1999 and Schaap, Leij & van Genuchten, 2001 (Rosetta software) have provided "classes" pedotransfer functions that could be considered as "typical" properties for given texture classes (FAO and USDA). These are not implemented yet in soilwaterptf but...

it is actually possible to "play" with the functions and data provided with the package soilwaterfun, soilwaterptf and soiltexture to "plot" the typical theta(h) and K(h).

I attach a detailed R script that show how this can be done.

all the best

Julien

Soil properties DB [ Reply ] By: Mikhail Titov on 2010-10-06 18:24	[forum:3458]
It would be nice to store some properties for at least some typical/generic soils for test purposes. Saying that I mean the possibility to get typical/generic curves for some clay, sand, etc.