Tidy Soil Data
library(aqp)
library(dplyr)
fp <- "C:/Users/Stephen.Roecker/NextCloud/projects/thesis-wvu"
s <- read.csv(file.path(fp, "siteAttributes.csv"), stringsAsFactors = FALSE)
s[1:4] <- lapply(s[1:4], as.numeric)
h <- read.csv(file.path(fp, "horizonAttributes.csv"), stringsAsFactors = FALSE)
h$CaMg <- with(h, Ca + Mg)
spc <- h
depths(spc) <- upedonid ~ hzdept + hzdepb
h_s <- aqp::slice(spc, 0:150 ~ fragvol + clay + sand + C + pH + CaMg, just.the.data = TRUE)
h_s$dep_int <- cut(h_s$hzdepb,
breaks = c(0, 15, 60, 100, 150),
labels = c("0-15", "15-60", "60-100", "100-150")
)
h_di <- group_by(h_s, upedonid, dep_int) %>%
summarize(fragvol = mean(fragvol, na.rm = TRUE),
clay = mean(clay, na.rm = TRUE),
C = sum(C),
pH = mean(pH, na.rm = TRUE),
CaMg = sum(CaMg),
hzthk = sum(!is.na(.pctMissing))
)%>%
mutate(fragvol = ifelse(fragvol <= 0.1, 0.1, fragvol),
C = ifelse(C == 0, NA, C),
CaMg = ifelse(CaMg == 0, NA, CaMg)
)
sh_di <- merge(s, h_di, by = "upedonid", all.x = TRUE)
Depth Plots
library(ggplot2)
library(GGally)
# Depth Plot
h_slab <- slab(spc, ~ fragvol + clay + C + pH + CaMg)
ggplot(h_slab, aes(x = bottom, y = p.q50)) +
geom_line() +
geom_ribbon(aes(ymin = p.q5, ymax = p.q95, x = bottom), alpha = 0.2) +
geom_ribbon(aes(ymin = p.q25, ymax = p.q75, x = bottom), alpha = 0.2) +
xlim(150, 0) +
facet_wrap(~ variable, scales = "free_x") +
coord_flip() +
xlab("depth (cm)") +
ylab("5th, 25th, Median, 75th, and 95th Quantiles") +
ggtitle("Depth Plot of Soil Properties")
# Scatter Plot Matrix
h$CaMg_log <- log(h$CaMg + 0.1)
sh_di$CaMg_log <- log(sh_di$CaMg + 0.1)
vars <- c("fragvol", "clay", "C", "pH", "CaMg_log")
ggpairs(h[vars])
ggpairs(sh_di[vars])
Sample Geodata and Correlate with Soil Depth Intervals
library(sp)
library(sf)
library(raster)
pts <- sh_di
coordinates(pts) <- ~ utm_easting + utm_northing
proj4string(pts) <- CRS("+init=epsg:26917")
pts2 <- st_sf(sh_di,
geometry = st_sfc(st_multipoint(as.matrix(sh_di[3:4]))),
crs = "+init=epsg:26917"
)
# stack SAGA rasters
load(file = "C:/workspace2/github/thesis-wvu/trunk/scale/geodata_df.RData")
sg <- {
subset(geodata, ! var %in% c("elev", "cucon") & ! grepl("slopeR|slopeD", var) & ! ws %in% c(27, 45, 63, 81)) ->.;
split(., .$res) ->.;
lapply(., function(x) {
cat("raster stacking and extracting", x$res[1], "\n")
rs = stack(x$sdat)
sg = as.data.frame(raster::extract(rs, pts, sp = TRUE))
# compute northness
idx = grepl("aspect", names(sg))
sg[idx] = lapply(sg[idx], function(x) abs(180 - x))
return(sg)
}) ->.;
}
sg2 <- {
lapply(sg, function(x) {
cat("correlating, stacking, and tidying", x$res[1]
# split by depth interval
split(x, x$dep_int) ->.;
# compute correlation and convert to long format
lapply(., function(x2) {
fragvol = data.frame(hzthk = x2$hzthk, fragvol = log(x2$fragvol + 0.1), x2[16:ncol(x2)])
clay = data.frame(hzthk = x2$hzthk, clay = x2$clay, x2[16:ncol(x2)])
C = data.frame(hzthk = x2$hzthk, C = x2$C, x2[16:ncol(x2)])
pH = data.frame(hzthk = x2$hzthk, pH = x2$pH, x2[16:ncol(x2)])
CaMg = data.frame(hzthk = x2$hzthk, CaMg = log(x2$CaMg + 0.1), x2[16:ncol(x2)])
# fragvol
idx = complete.cases(fragvol)
frag_cor = cov.wt(x = fragvol[idx, ], wt = fragvol$hzthk[idx], cor = TRUE)$cor[- c(1:2), 2]
# clay
idx = complete.cases(clay)
clay_cor = cov.wt(x = clay[idx, ], wt = clay$hzthk[idx], cor = TRUE)$cor[- c(1:2), 2]
# C
idx = complete.cases(C)
C_cor = cov.wt(x = C[idx, ], wt = C$hzthk[idx], cor = TRUE)$cor[- c(1:2), 2]
# pH
idx = complete.cases(pH)
pH_cor = cov.wt(x = pH[idx, ], wt = pH$hzthk[idx], cor = TRUE)$cor[- c(1:2), 2]
# CaMg
idx = complete.cases(CaMg)
CaMg_cor = cov.wt(x = CaMg[idx, ], wt = CaMg$hzthk[idx], cor = TRUE)$cor[- c(1:2), 2]
test = data.frame(values = c(fragvol = frag_cor, clay = clay_cor, C = C_cor, pH = pH_cor, CaMg = CaMg_cor))
test$ind = row.names(test)
# append depth interval
test$dep_int = x2$dep_int[1]
return(test)
}) -> .;
do.call("rbind", .) ->test;
# split property.filename grouping variable into separate columns
test = cbind(
test,
data.frame(
do.call("rbind",
strsplit(test$ind, "\\.|_")
),
stringsAsFactors = FALSE
)
)
names(test)[4:7] = c("prop", "source", "area", "var")
# split source and variable columns into additional columns
test = within(test, {
res = substr(source, 5, 8)
res = as.numeric(substr(res, 1, nchar(res) - 1))
source = substr(source, 1, 4)
ws = sub("slope|aspect|cupro|cutan", "", var)
ws = as.numeric(substr(ws, 1, nchar(ws) - 1))
ns = ws * res
res = factor(res, levels = sort(unique(res)))
dep_int = factor(paste(dep_int, "cm"), levels = c("0-15 cm", "15-60 cm", "60-100 cm", "100-150 cm"))
var[grepl("slope", var)] = "slope"
var[grepl("aspect", var)] = "northness"
var[grepl("cupro", var)] = "cupro"
var[grepl("cutan", var)] = "cutan"
})
}) ->.;
do.call("rbind", .) ->.;
}
save(sg2, file = "C:/Users/Stephen.Roecker/NextCloud/projects/thesis-wvu/soil_correlations.RData")
Create Line Plots
All DEM Grid and Neighborhood Size Combinations
load(file = "C:/Users/Stephen.Roecker/NextCloud/projects/thesis-wvu/soil_correlations.RData")
filter(sg2, ! res %in% c("6", "45", "81")) %>%
ggplot(aes(x = ns, y = values, col = var, shape = res)) +
geom_point(size = 2.5, alpha = 0.5) +
geom_line(lwd = 1, alpha = 0.5) +
geom_hline(yintercept = 0) +
facet_grid(prop ~ dep_int) +
theme(aspect.ratio = 1) +
scale_shape_discrete(name = "grid size\n(meters)") +
scale_color_discrete(name="terrain\nattribute")+
xlab("neighborhood size (meters)") +
ylab("correlation coefficient (r)") +
labs(caption = "log transformation applied to fragvol and CaMg") +
ggtitle("Correlation Between Soil Properties and Terrain Attributes\nOver 4 Depth Intervals and Several Grid and Neighborhood Sizes")
9-meter DEM and Neighborhood Size Combinations
filter(sg2, res == "9") %>%
ggplot(aes(x = ns, y = values, col = var, shape = res)) +
geom_point(size = 2.5) +
geom_line(lwd = 1) +
geom_hline(yintercept = 0) +
facet_grid(prop ~ dep_int) +
theme(aspect.ratio = 1) +
scale_shape_discrete(name = "grid size\n(meters)") +
scale_color_discrete(name="terrain\nattribute")+
xlab("neighborhood size (meters)") +
ylab("correlation coefficient (r)") +
labs(caption = "log transformation applied to fragvol and CaMg") +
ggtitle("Correlation Between Soil Properties and Terrain Attributes\nOver 4 Depth Intervals and Several Neighborhood Sizes")
All DEM Grid Sizes and 3x3 Window Size Combinations
filter(sg2, ws == 3) %>%
ggplot(aes(x = ns, y = values, col = var)) +
geom_point(aes(shape = res), size = 2.5) +
geom_line(lwd = 1) +
geom_hline(yintercept = 0) +
facet_grid(prop ~ dep_int) +
theme(aspect.ratio = 1) +
scale_shape_discrete(name = "grid size\n(meters)") +
scale_color_discrete(name="terrain\nattribute")+
xlab("neighborhood size (meters)") +
ylab("correlation coefficient (r)") +
labs(caption = "log transformation applied to fragvol and CaMg") +
ggtitle("Correlation Between Soil Properties and Terrain Attributes\nOver 4 Depth Intervals and Several Grid Sizes")
