Skip to content

Commit

Permalink
Update mps2-eval.R
Browse files Browse the repository at this point in the history
  • Loading branch information
@dylanbeaudette
dylanbeaudette committed Nov 1, 2023
1 parent 9eb0ca0 commit b79019b
Showing 1 changed file with 106 additions and 17 deletions.
123 changes: 106 additions & 17 deletions misc/sandbox/mps2-eval.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ library(sharpshootR)
library(mpspline2)
library(reshape2)
library(lattice)
library(tactile)

## thoughts:
# * variable depths option from EAS
Expand All @@ -13,20 +14,94 @@ library(lattice)
# * slab2SPC() or option to slab() to re-shape into an SPC


# make some example data
x <- list(
id = 'A',
depths = c(25, 33, 100, 150, 175),
name = c('A', 'Bw', 'Bt1', 'Bt2', 'BC'),
p1 = c(12, 15, 35, 20, 15),
soil_color = c('10YR 3/3', '10YR 4/4', '10YR 4/6', '5G 6/2', '5BG 4/4')
)

x <- quickSPC(x)

# fake site data
site(x)$fake_site_attr <- 1:length(x)

# check source data
par(mar=c(0,0,3,1))
plotSPC(x, color='p1', col.palette = hcl.colors(10, 'viridis'))

s <- seq(0, 1, by = 0.25)
m <- lapply(s, function(i) {
.spc <- spc2mpspline(x, var_name = 'p1', lam = i, method = 'est_1cm')
profile_id(.spc) <- sprintf("%s-%s", profile_id(.spc), format(i, digits = 2))
return(.spc)
})

z <- combine(m)
z$p1 <- z$p1_spline

## don't need these anymore
# helper functions
# source('https://raw.githubusercontent.com/ncss-tech/aqp/master/misc/sandbox/mps-functions.R')
xx <- combine(z, x)

# fewer moving parts
options(stringsAsFactors = FALSE)
par(mar = c(0, 0, 3, 2))
plotSPC(xx, color = 'p1', col.palette = hcl.colors(10, 'viridis'), divide.hz = FALSE, width = 0.35, cex.names = 0.8, name = NA)
mtext(side = 3, at = 2, line = -1.5, text = 'value of lambda')

# o <- order(xx$p1_rmse)
# plotSPC(xx, color = 'p1', col.palette = hcl.colors(10, 'viridis'), divide.hz = FALSE, width = 0.35, plot.order = o)

d <- dice(xx)
d$p1diff <- rep(NA, times = nrow(d))

d$p1diff <- profileApply(d, function(i) {
d[1, ]$p1 - i$p1
})

plotSPC(d, color = 'p1diff', col.palette = hcl.colors(10, 'viridis'), divide.hz = FALSE, width = 0.35, cex.names = 0.8, name = NA)
mtext(side = 3, at = 2, line = -2, text = 'value of lambda')



tps <- tactile.theme(superpose.line = list(lwd = 2))

# compare depth-functions by method, no aggregation
xyplot(cbind(top, bottom) ~ p1, id=xx$id, groups = factor(id),
data=as(xx, 'data.frame'),
par.settings=tps,
auto.key=list(columns=3, lines=TRUE, points=FALSE),
strip=strip.custom(bg=grey(0.85)),
ylim=c(185, -5), as.table=TRUE, panel=panel.depth_function,
scales = list(alternating = 1),
asp = 1
)

xyplot(cbind(top, bottom) ~ p1 | id, id=xx$id,
data=as(xx, 'data.frame'),
par.settings=tps,
auto.key=list(columns=4, lines=TRUE, points=FALSE),
strip=strip.custom(bg=grey(0.85)),
ylim=c(185, -5), as.table=TRUE, panel=panel.depth_function,
layout = c(6, 1),
scales = list(alternating = 1)
)



## test effect of lambda
slab(xx, id ~ p1, slab.structure = c(0, 50), slab.fun = mean, na.rm = TRUE)
slab(xx, id ~ p1, slab.structure = c(0, 50), slab.fun = sum, na.rm = TRUE)




####################### more complex examples / eval ########################

# make some example data
ids <- LETTERS[1:6]

set.seed(10101)
x <- lapply(ids, random_profile, n=c(6, 7, 8), n_prop=2, method='LPP', SPC=TRUE)
x <- lapply(ids, random_profile, n = c(6, 7, 8), n_prop = 2, method = 'LPP', SPC = TRUE, lpp.u = 40)
x <- combine(x)

# fake site data
Expand Down Expand Up @@ -127,23 +202,23 @@ plotSPC(z, color='p1', divide.hz = FALSE, width = 0.25, cex.names = 0.66, max.de

# plot by group
par(mar=c(0, 0, 3, 1))
plotSPC(z, color='p1', max.depth=175, name=NA, divide.hz=FALSE, width = 0.3, col.palette = hcl.colors(10), n.legend = 8)
plotSPC(z, color='p1', max.depth = max(x), name=NA, divide.hz=FALSE, width = 0.3, col.palette = hcl.colors(10), n.legend = 8)

groupedProfilePlot(z, groups = 'original_id', color='p1', max.depth = max(x), group.name.offset = -10, name=NA, divide.hz=FALSE, width=0.3, col.palette = hcl.colors(10), n.legend = 8)

groupedProfilePlot(z, groups = 'provenance', color='p1', max.depth = max(x), group.name.offset = -10, name=NA, divide.hz=FALSE, width = 0.3, col.palette = hcl.colors(10), n.legend = 8)

tps <- list(superpose.line=list(lwd=2, col=c('royalblue', 'darkgreen', 'firebrick', 'orange')))

# compare depth-functions by method, no aggregation
xyplot(cbind(top, bottom) ~ p1 | original_id, id=z$id, groups=provenance,
data=as(z, 'data.frame'),
par.settings=tps,
auto.key=list(columns=4, lines=TRUE, points=FALSE),
strip=strip.custom(bg=grey(0.85)),
ylim=c(175, -5), as.table=TRUE, panel=panel.depth_function,
ylim=c(125, -5), as.table=TRUE, panel=panel.depth_function,
layout = c(6, 1),
scales = list(x = list(relation = 'free'))
scales = list(x = list(relation = 'free')),
ylab = 'Depth (cm)',
xlab = 'simulated property'
)


Expand All @@ -152,7 +227,7 @@ xyplot(cbind(top, bottom) ~ p1 | original_id, id=z$id, groups=provenance,
z.agg <- slab(z, fm = provenance ~ p1)

xyplot(top ~ p.q50, groups=provenance, data=z.agg, ylab='Depth', asp=1.5,
lower=z.agg$p.q25, upper=z.agg$p.q75, ylim=c(200,-5),
lower=z.agg$p.q25, upper=z.agg$p.q75, ylim=c(125,-5),
xlab='median bounded by 25th and 75th percentiles',
sync.colors=TRUE, alpha=0.25,
par.settings=tps,
Expand All @@ -164,7 +239,7 @@ xyplot(top ~ p.q50, groups=provenance, data=z.agg, ylab='Depth', asp=1.5,
)

xyplot(top ~ p.q50, groups=provenance, data=z.agg, ylab='Depth', asp=1.5,
ylim=c(200,-5),
ylim=c(125,-5),
xlab='median bounded by 25th and 75th percentiles',
par.settings=tps,
auto.key=list(columns=4, lines=TRUE, points=FALSE),
Expand All @@ -174,11 +249,11 @@ xyplot(top ~ p.q50, groups=provenance, data=z.agg, ylab='Depth', asp=1.5,
)

xyplot(top ~ p.q50 | provenance, data=z.agg, ylab='Depth', asp=1.5,
lower=z.agg$p.q25, upper=z.agg$p.q75, ylim=c(175,-5),
lower=z.agg$p.q25, upper=z.agg$p.q75, ylim=c(125, -5),
xlab='median bounded by 25th and 75th percentiles',
sync.colors=TRUE, alpha=0.5,
strip=strip.custom(bg=grey(0.85)),
par.settings=list(superpose.line=list(lwd=2, col=c('RoyalBlue', 'firebrick'))),
par.settings = tps,
panel=panel.depth_function,
prepanel=prepanel.depth_function,
cf=z.agg$contributing_fraction,
Expand All @@ -188,11 +263,25 @@ xyplot(top ~ p.q50 | provenance, data=z.agg, ylab='Depth', asp=1.5,

# * only using a single property.. need to update EAS code to do more than a single property
# * truncate comparisons to 110cm
d <- NCSP(z, vars = 'p1', maxDepth = 110)
d <- NCSP(z, vars = 'p1', maxDepth = 100)

# interesting
par(mar = c(0, 0, 0, 3))
plotProfileDendrogram(z, cluster::diana(d), dend.y.scale = 175, scaling.factor = 0.85, y.offset = 2, color = 'p1', divide.hz = FALSE, width = 0.3, max.depth = 150, cex.names = 0.8)
plotProfileDendrogram(z, cluster::diana(d), dend.y.scale = 175, scaling.factor = 0.95, y.offset = 2, color = 'p1', divide.hz = FALSE, width = 0.3, max.depth = 125, cex.names = 0.8)


## TODO: viz via nMDS and betadispersion


## what about simple summaries within a fixed depth interval?
# nearly identical

slab(z, provenance ~ p1, slab.structure = c(0, 50), slab.fun = mean, na.rm = TRUE)

slab(z, provenance ~ p1, slab.structure = c(0, 50), slab.fun = sum, na.rm = TRUE)






0 comments on commit b79019b

Please sign in to comment.