diff --git a/misc/estimateAWC-vs-ROSETTA-centroids.R b/misc/estimateAWC-vs-ROSETTA-centroids.R new file mode 100644 index 000000000..902d6a9c0 --- /dev/null +++ b/misc/estimateAWC-vs-ROSETTA-centroids.R @@ -0,0 +1,59 @@ +library(aqp) +library(soilDB) +library(lattice) +library(tactile) + +data("ROSETTA.centroids") + +x <- ROSETTA.centroids + +# hmm.. why does this give NA? +estimateAWC('si', omcl = 1) + + +m <- data.frame( + names = SoilTextureLevels(which = 'names'), + codes = SoilTextureLevels(which = 'codes') +) + +idx <- match(x$texture, m$names) +x$txtcl <- m$codes[idx] + + +x$est.awc <- estimateAWC(x$txtcl, omcl = rep(1, times = nrow(x))) + + +xyplot( + est.awc ~ awc, + data = x, + xlim = c(-0.025, 0.23), + ylim = c(-0.025, 0.23), + xlab = 'AWC via ROSETTA Centroids', + ylab = 'Estimated AWC', + par.settings = tactile.theme(), + asp = 1, + scales = list(alternating = 3, tick.number = 10), + pch = 16, + panel = function(...) { + panel.grid(-1, -1) + panel.xyplot(...) + panel.text(x = x$awc, y = x$est.awc, labels = x$texture, cex = 0.75 , pos = 1) + panel.abline(lm(est.awc ~ awc, data = x)) + }) + + + +x$diff <- with(x, awc - est.awc) + +dotplot( + texture ~ diff, + data = x, + xlab = 'ROSETTA AWC - Estimated AWC', + par.settings = tactile.theme(), + cex = 1.5, + panel = function(...) { + panel.abline(v = 0, lty = 2) + panel.dotplot(...) + }) + +