From 2c75726e1a6d74e4888c159bedd1f9d4992c4744 Mon Sep 17 00:00:00 2001 From: "Andrew G. Brown" Date: Thu, 7 Dec 2023 07:44:52 -0800 Subject: [PATCH] simple sda example/comparison v.s. stored Ksat/texture class --- misc/sda-exmaple.R | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 misc/sda-exmaple.R diff --git a/misc/sda-exmaple.R b/misc/sda-exmaple.R new file mode 100644 index 0000000..bb46800 --- /dev/null +++ b/misc/sda-exmaple.R @@ -0,0 +1,35 @@ +library(soilDB) +library(rosettaPTF) + +x <- SDA_query("SELECT TOP 100000 component.cokey, chorizon.chkey, + texture, + wtenthbar_l, wtenthbar_r, wtenthbar_h, + wthirdbar_l, wthirdbar_r, wthirdbar_h, + wfifteenbar_l, wfifteenbar_r, wfifteenbar_h, + ---wsatiated_l, wsatiated_r, wsatiated_h, + ksat_l, ksat_r, ksat_h, + awc_l, awc_r, awc_h, + sandtotal_l, sandtotal_r, sandtotal_h, + silttotal_l, silttotal_r, silttotal_h, + claytotal_l, claytotal_r, claytotal_h, + dbthirdbar_l, dbthirdbar_r, dbthirdbar_h + FROM component + INNER JOIN chorizon ON component.cokey = chorizon.cokey + INNER JOIN chtexturegrp ON chorizon.chkey = chtexturegrp.chkey AND rvindicator = 'Yes'") +x$id <- 1:nrow(x) +res <- run_rosetta(x[,c("sandtotal_r", "silttotal_r", "claytotal_r", + "dbthirdbar_r", "wthirdbar_r", "wfifteenbar_r")]) +x <- merge(x, res, by = "id") +View(x) +x$ksat_pred_r <- 10^(x$log10_Ksat_mean) / 24 # cm/day -> cm/hr +x$ksat_pred_ratio <- x$ksat_pred_r / x$ksat_r +plot(x$ksat_r, x$ksat_pred_r, xlab="SSURGO (cm/h)",ylab="ROSETTA (cm/h)", + xlim=c(0,100), ylim=c(0,100)) +m<-lm(x$ksat_pred_r~x$ksat_r) +abline(m, lty=2) +abline(0,1) +summary(m) +x2 <- subset(x, !grepl("-|SR|BR|WB|DUR|VAR|PT|PM|GR|CB|FRAG|PEAT|SG|IND|CEM|GYP|CPF|MARL|G|CE|MUCK", x$texture)) +plot(x2$ksat_pred_r~factor(x2$texture), ylim=c(0,200)) +plot(x2$ksat_r~factor(x2$texture),ylim=c(0,200)) +plot(x2$ksat_pred_ratio~factor(x2$texture))