model_training.Rmd

---
title: "GSP-CONUS"
author: "Stephen Roecker"
date: "8/4/2020"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```



# Extract geodata

```{r}

library(raster)

# ISRIC CONUS
path <- "D:/geodata/project_data/gsp-sas/1km covariates/ISRIC/CONUS"
lf <- list.files(path)
lf <- lf[grepl(".tif$", lf)]
isric_rs <- stack(file.path(path, lf))
names(isric_rs) <- gsub("_CONUS.tif|.tif", "", lf)
isric_vars <- names(isric_rs)


# Other CONUS
path <- "D:/geodata/project_data/gsp-sas/1km covariates/Other"
lf <- list.files(path)
lf <- lf[grepl(".tif$", lf)]
# find matching extents
# other_extent <- sapply(file.path(path, lf), function(x) paste(extent(raster(x))))
# idx <- names(sort(table(other_extent), decreasing = TRUE)[1])
# idx <- which(other_extent == idx)
# lf  <- lf[idx]

other_rs <- stack(file.path(path, lf))
names(other_rs) <- gsub(".tif", "", lf)
other_vars <- names(other_rs)


# SSURGO CONUS
path <- "D:/geodata/project_data/gsp-sas/1km covariates/SSURGO"
lf <- list.files(path)
lf <- lf[grepl(".tif$", lf)]
# # find matching extents
# other_extent <- sapply(file.path(path, lf), function(x) paste(extent(raster(x))))
# idx <- names(sort(table(other_extent), decreasing = TRUE)[1])
# idx <- which(other_extent == idx)
# lf  <- lf[idx]

ssurgo_rs <- stack(file.path(path, lf))
names(ssurgo_rs) <- gsub(".tif", "", lf)
ssurgo_vars <- names(ssurgo_rs)


rs   <- stack(isric_rs, other_rs, ssurgo_rs)



# Extract geodata
library(sf)

load(file = "C:/Users/stephen.roecker/Nextcloud/projects/2020_gsp-sas/LDM-compact_20200709.RData")

data <- extract(rs, as(s_mps_sf, "Spatial"), df = TRUE, sp = TRUE)@data
data <- within(data, {
  ssurgo_pm_lu3_rs_mode_v2 = as.factor(ssurgo_pm_lu3_rs_mode_v2)
  statsgo_mlra             = as.factor(statsgo_mlra)
  statsgo_pmkind3          = as.factor(statsgo_pmkind3)
  mlra_rs                  = as.factor(mlra_rs)
  nlcd_2016_conus          = as.factor(nlcd_2016_conus)
  })

```



# Data Splitting

```{r}

# split dataset into training and testing
vars <- names(data)[grepl("ph_ptf.|ec_ptf.|esp.", names(data))]
vars <- vars[grepl("030", vars)]
n    <- round(nrow(data) * 0.75)
set.seed(111)
clhs_idx  <- clhs::clhs(data[vars], size = n)
clhs_pkey <- data[clhs_idx, "pedon_key"]

train <- subset(data, pedon_key %in% clhs_pkey)
test  <- subset(data, ! pedon_key %in% clhs_pkey)

# save(data, train, test, clhs_pkey, isric_vars, other_vars, ssurgo_vars, file = "C:/Users/Stephen.Roecker/Nextcloud/projects/2020_gsp-sas/training_data.RData")

```



# Feature Selection

```{r}

load(file = "C:/Users/Stephen.Roecker/Nextcloud/projects/2020_gsp-sas/training_data.RData")

# high correlations
library(caret)

vars <- c(isric_vars, other_vars, ssurgo_vars)
idx <- sapply(data, is.numeric) & names(data) %in% vars
data2 <- na.exclude(data[idx])

vars_0 <- nearZeroVar(data2, names = TRUE)

data3 <- data2[!names(data2) %in% vars_0]
vars_cor <- findCorrelation(x = cor(data3, use = "complete.obs"), names = TRUE, cutoff = 0.75)


# Boruta
sp_vars <- names(data)[grepl("ph_ptf.|ec_ptf.|esp.", names(data)) & grepl("030", names(data))]
vars <- c(isric_vars, other_vars, ssurgo_vars)

set.seed(111)

# pick your variables with idx
idx <- 1:2
fs_b <- lapply(sp_vars[idx], function(x) {
  cat(x, as.character(Sys.time()), "\n")
  
  # remove matching ssurgo variable
  x2     <- strsplit(x, split = "_|\\.")[[1]][1]
  svars2 <- ssurgo_vars[grepl(x2, ssurgo_vars)]
  vars2  <- vars[!vars %in% svars2]
  
  train2 <- data[c(x, vars2)]
  train2 <- na.exclude(train2)
  names(train2) <- gsub("\\.", "_", names(train2))
  
  # tranform y if ec or esp
  y <- if (!grepl("^ph_", x)) log(train2[[1]] +0.1) else train2[[1]]
  
  test <- Boruta::Boruta(x = train2[-1], y = y, maxRuns = 35, doTrace = 1)
})
names(fs_b) <- sp_vars[idx]
# save(fs_b_ph, file = "C:/Users/Stephen.Roecker/Nextcloud/projects/2020_gsp-sas/fs_boruta_ph.RData")
load(file = "C:/Users/Stephen.Roecker/Nextcloud/projects/2020_gsp-sas/fs_boruta_ph.RData")
load(file = "C:/Users/Stephen.Roecker/Nextcloud/projects/2020_gsp-sas/fs_boruta_ec_esp.RData")

plot(fs_b_ph$ph_ptf.000_030_cm)

fs_b_stats <- lapply(fs_b_ph, function(x) {
  s <- Boruta::attStats(x)
  s <- data.frame(variable = row.names(s), varImp = s$medianImp, decision = s$decision, stringsAsFactors = FALSE)
  s$variable <- reorder(as.factor(s$variable), s$varImp, sum)
  s <- s[order(- s$varImp), ]
  row.names(s) <- 1:nrow(s)
  return(s)
  })

library(ggplot2)
ggplot(fs_b_stats$ph_ptf.000_030_cm[1:25, ], aes(x = varImp, y = variable)) +
  geom_point()


# recursive feature elimination
set.seed(111)
library(caret)
rfFuncs$fit <- function(x, y, first, last, ...) {
  train(x, y, 
               method = "ranger", quantreg = TRUE, importance = "permutation",
               # trControl = trainControl(method = "none"), 
               tuneGrid = expand.grid(.mtry = round(sqrt(ncol(y))), .splitrule = "variance", .min.node.size = c(5)
)
               )
}

ctrl.RFE <- rfeControl(functions = rfFuncs,
                       method = "cv",
                       number = 5, 
                       verbose = TRUE)

train3 <- na.exclude(train2)
x = train3[, -1]
y = train3[,  1]
rf.RFE <- rfe(x = x,
              y = y,
              sizes = 5:20,
              rfeControl = ctrl.RFE
              )


# http://www.itgo.me/a/2164573104609813273/r-carets-rfe-error-in-task-1-failed-rfe-is-expecting-184-importance-val
 
# save(data, train, test, test30, clhs_pkey, file = "C:/Users/Stephen.Roecker/Nextcloud/projects/2020_gsp-sas/training_data.RData")

```



# Model Training

```{r, single-rf}

library(ranger)
library(caret)

load(file = "C:/Users/Stephen.Roecker/Nextcloud/projects/2020_gsp-sas/training_data.RData")

# single forest
vars <- c("ph_ptf.000_030_cm",
          isric_vars, 
          other_vars,
          # remove matching ssurgo derivative
          c("statsgo_mlra", "statsgo_pmkind3")
          )
train2 <- na.exclude(train[vars])



# Boruta confirmed variables
idx <- which(names(train2) %in% fs_b_stats$ph_ptf.000_030_cm$variable[fs_b_stats$ph_ptf.000_030_cm$decision == "Confirmed"])
# log transform y if necessary (e.g. log(y + 0.1))
ph_rf_all <- ranger(y = train2$ph_ptf.000_030_cm, x = train2[idx], 
                    quantreg = TRUE,importance = "permutation",
                    respect.unordered.factors = FALSE,
                    seed = 111 
                    )

# top 50 variables
idx <- which(names(train2) %in% fs_b_stats$ph_ptf.000_030_cm$variable[1:50])
ph_rf_t50 <- ranger(y = train2$ph_ptf.000_030_cm, x = train2[idx], 
                    quantreg = TRUE,importance = "permutation",
                    respect.unordered.factors = NULL,
                    seed = 111 
                    )

# top 25 variables
idx <- which(
  names(train2) %in%
    c(as.character(fs_b_stats$ph_ptf.000_030_cm$variable[1:23]),  
      "statsgo_mlra", "statsgo_pmkind3")
  )
ph_rf_t25 <- ranger(ph_ptf.000_030_cm ~ ., data = train2[c(1, idx)], 
                    importance = "permutation", quantreg = TRUE, # importance = "impurity",
                    mtry = 8, min.node.size = 6,
                    respect.unordered.factors = NULL,
                    seed = 111 
                    )


# uncorrelated variables
idx <- names(train2) %in% fs_b_stats$ph_ptf.000_030_cm$variable[fs_b_stats$ph_ptf.000_030_cm$decision == "Confirmed"] & sapply(train2, is.numeric)
vars_cor <- findCorrelation(x = cor(train2[idx], use = "complete.obs"), names = TRUE, cutoff = 0.75)
vars <- c(vars_cor, ssurgo_vars[1:2], "nlcd_2016_conus")
ph_rf_cor <- ranger(y = train2$ph_ptf.000_030_cm, x = train2[vars_cor], 
                    quantreg = TRUE,importance = "permutation",
                    seed = 111 
                    )




# variable importance plots
vip::vip(ph_rf_t25, 50)


# residual plots
library(DALEX)

ph_rf <- ph_rf_t25
ph_rf_e <- explain(ph_rf, data = train2[-1], y = train2$ph_ptf.000_030_cm)
ph_rf_p <- model_performance(ph_rf_e)
ph_rf_d <- model_diagnostics(ph_rf_e)
plot(ph_rf_d, variable = "y", yvariable = "y_hat") + geom_abline(intercept = 0, slope = 1)
plot(ph_rf_d, variable = "y", yvariable = "residuals")



# model tuning
# pH 0-30 cm
vars <- c("ph_ptf.000_030_cm",
          isric_vars, 
          other_vars,
          # remove matching ssurgo derivative
          c("statsgo_mlra", "statsgo_pmkind3")
          )
train2 <- na.exclude(train[vars])

idx <- which(names(train2) %in%
               c(as.character(fs_b_stats$ph_ptf.000_030_cm$variable[1:23]),
                 c("statsgo_mlra", "statsgo_pmkind3")
                 ))
ph_rf_030_v2 <- train(y = train2$ph_ptf.000_030_cm, x = train2[idx], 
               method = "ranger", quantreg = TRUE, seed = 111,
               importance = "permutation",
               trControl = trainControl(method = "cv", number = 10, returnResamp = "all", savePredictions = TRUE, search = "random", verboseIter = FALSE)
               )

# pH 30-100 cm
vars <- c("ph_ptf.030_100_cm",
          isric_vars, 
          other_vars,
          # remove matching ssurgo derivative
          c("statsgo_mlra", "statsgo_pmkind3")
          )
train2 <- na.exclude(train[vars])

idx <- which(names(train2) %in%
               c(as.character(fs_b_stats$ph_ptf.030_100_cm$variable[1:23]),
                 c("statsgo_mlra", "statsgo_pmkind3")
                 ))
 
ph_rf_100_v2 <- train(y = train2$ph_ptf.030_100_cm, x = train2[idx], 
               method = "ranger", quantreg = TRUE, seed = 111,
               importance = "permutation",
               # na.action = na.exclude,
               trControl = trainControl(method = "cv", number = 10, returnResamp = "all", savePredictions = TRUE, search = "random", verboseIter = FALSE)
               )

print(ph_rf_030_v2)
print(ph_rf_100_v2)

# save(ph_rf_030_v2, ph_rf_100_v2, file = "C:/Users/Stephen.Roecker/Nextcloud/projects/2020_gsp-sas/ph_rf.RData")

predfun <- function(model, ...) predict(model, ...)$predictions

idx <- which(
  names(rs) %in% ph_rf_030_v2$finalModel$forest$independent.variable.names
  )
rs2 <- rs[[idx]]
ph_030_r <- predict(rs2, ph_rf_030_v2$finalModel, fun = predfun, index = 1, progress = "text", overwrite = TRUE, filename = "D:geodata/project_data/gsp-sas/ph_030_qrf_mlra_pm_NULL_final4.tif")


idx <- which(
  names(rs) %in% ph_rf_100_v2$finalModel$forest$independent.variable.names
  )
rs2 <- rs[[idx]]
ph_100_r <- predict(rs2, ph_rf_100_v2$finalModel, fun = predfun, index = 1, progress = "text", overwrite = TRUE, filename = "D:geodata/project_data/gsp-sas/ph_100_qrf_mlra_pm_NULL_final4.tif")


```


# Fill gaps and mask

```{r}

library(raster)

ph_030_r <- raster("D:geodata/project_data/gsp-sas/ph_030_qrf_mlra_pm_NULL_final4.tif")
ph_100_r <- raster("D:geodata/project_data/gsp-sas/ph_100_qrf_mlra_pm_NULL_final4.tif")

ph_030_r0 <- raster("D:geodata/project_data/gsp-sas/1km covariates/SSURGO/ssurgo_ptph30_rs_med_v2.tif")
ph_100_r0 <- raster("D:geodata/project_data/gsp-sas/1km covariates/SSURGO/ssurgo_ptph100_rs_med_v2.tif")


ph_030_r2 <- merge(ph_030_r, ph_030_r0, overlap = TRUE, filename = "D:geodata/project_data/gsp-sas/ph_030_qrf_mlra_pm_NULL_final4_fill.tif", progress = "text", overwrite = TRUE)
ph_100_r2 <- merge(ph_100_r, ph_100_r0, overlap = TRUE, filename = "D:geodata/project_data/gsp-sas/ph_100_qrf_mlra_pm_NULL_final4_fill.tif", progress = "text", overwrite = TRUE)

# mask
lc <- raster("D:/geodata/project_data/gsp-sas/1km covariates/ISRIC/CONUS_redundant/LCEE10__CONUS.tif")
m <- c(
  0,   209, 1, 
  209, 210, NA,
  210, 300, 1
  )
m <- matrix(m, ncol = 3, byrow = TRUE)
lc <- reclassify(lc, m, filename = "D:/geodata/project_data/gsp-sas/1km covariates/ISRIC/CONUS_redundant/test.tif", overwrite = TRUE, progress = TRUE)

mask(ph_030_r2, lc, filename = "D:geodata/project_data/gsp-sas/ph_030_qrf_mlra_pm_NULL_final4_fill_mask.tif", overwrite = TRUE)
mask(ph_100_r2, lc, filename = "D:geodata/project_data/gsp-sas/ph_100_qrf_mlra_pm_NULL_final4_fill_mask.tif", overwrite = TRUE)

```


# Model Uncertainty

```{r}

load(file = "C:/Users/Stephen.Roecker/Nextcloud/projects/2020_gsp-sas/training_data.RData")
load(file = "D:/geodata/project_data/gsp-sas/deliverables/models/ph_rf.RData")
load(file = "D:/geodata/project_data/gsp-sas/deliverables/models/ec_rf.RData")
load(file = "D:/geodata/project_data/gsp-sas/deliverables/models/esp_rf.RData")
load(file = "C:/Users/stephen.roecker/Nextcloud/projects/2020_gsp-sas/LDM-compact_20200709.RData")


# plot
library(ggplot2)

models <- list(ph_rf_030_v2, ph_rf_100_v2, ec_rf_030, ec_rf_100, esp_rf_030_tune, esp_rf_100_tune)
vars <- list(
  "ph_ptf.000_030_cm", "ph_ptf.030_100_cm", 
  "ec_ptf.000_030_cm", "ec_ptf.030_100_cm",
  "esp.000_030_cm",    "esp.030_100_cm"
  )
fun <- function(variable, model) {
  idx     <- c("pedon_key", variable, model$finalModel$forest$independent.variable.names)
  test2   <- na.exclude(train[idx])
  
  test2$obs  <- test2[[variable]]
  test2$pred <- predict(model, newdata = test2)
  if (grepl(c("^ec|^esp"), variable)) {
    test2$pred <- exp(test2$pred) - 0.1
  }
  
  return(test2)
}
test2 <- mapply(FUN = fun, vars, models, SIMPLIFY = FALSE)
names(test2) <- unlist(vars)

# pH
ggplot(test2[[1]], aes(x = obs, y = pred)) +
  geom_point(alpha = 0.1) +
  geom_abline(col = "blue", size = 1)

# EC
ggplot(test2[[5]], aes(x = obs, y = pred)) +
  geom_point(alpha = 0.1) +
  geom_abline(col = "blue", size = 1) +
  scale_x_log10() + scale_x_continuous(breaks = c(0, 2, 8, 16, 200)) + xlim(0, 50) + ylim(0, 16) # EC



# accuracy metrics stratia
stratia_f <- function(predicted, measured) {

  data.frame(
    Bias     =      mean( predicted - measured,    na.rm = TRUE),
    RMSE     = sqrt(mean((predicted - measured)^2, na.rm = TRUE)),
    Rsquared =      cor(  predicted,  measured, use = "pairwise.complete.obs")^2,
    NSE      = 1 - (
      sum((measured - predicted)^2) /
      sum((measured - mean(measured))^2)
      )
    # NSE = hydroGOF::NSE(predicted, measured)
  )}


models <- list(ph_rf_030_v2, ph_rf_100_v2, ec_rf_030, ec_rf_100, esp_rf_030_tune, esp_rf_100_tune)
vars   <- list("pH", "EC", "ESP")[c(1, 1, 2, 2, 3, 3)]
deps   <- list(30, 100)[rep(1:2, 3)]
variables <- list(
  "ph_ptf.000_030_cm", "ph_ptf.030_100_cm", 
  "ec_ptf.000_030_cm", "ec_ptf.030_100_cm",
  "esp.000_030_cm",    "esp.030_100_cm"
  )
fun <- function(variable, var, depth, model, dataset, data) {
  
  idx     <- c("pedon_key", variable, model$finalModel$forest$independent.variable.names)
  data2   <- na.exclude(data[idx])
  data2$p <- predict(model, newdata = data2)
  # m <- strsplit(var, "_|\\.")[[1]][1]

  obs  <- data2[[variable]]
  pred <- data2$p
  if (grepl(c("^ec|^esp"), variable)) {
    pred <- exp(pred) - 0.1
  }

  temp = cbind(dataset = dataset,
               var,
               depth,
               stratia_f(predicted = pred, measured = obs)
               )

  return(temp)
}
test_stratia <- mapply(FUN = fun, variables, vars, deps, models, list("test"), list(test), SIMPLIFY = FALSE)
test_stratia <- do.call("rbind", test_stratia)
idx <- 4:7
test_stratia[idx] <- apply(test_stratia[idx], 2, round, 2)

train_stratia <- mapply(FUN = fun, variables, vars, deps, models, list("train"), list(train), SIMPLIFY = FALSE)
train_stratia <- do.call("rbind", train_stratia)
idx <- 4:7
train_stratia[idx] <- apply(train_stratia[idx], 2, round, 2)

# combine
val2 <- rbind(train_stratia, test_stratia)
val2 <- val2[order(val2$var, val2$depth), ]


write.csv(ph_030_stratia, file = "C:/Users/stephen.roecker/Nextcloud/projects/2020_gsp-sas/ph_030_validmodel_stats.csv")
write.csv(ph_100_stratia, file = "C:/Users/stephen.roecker/Nextcloud/projects/2020_gsp-sas/ph_100_validmodel_stats.csv")



# compare distributions
pkg <- c("sf", "sp", "raster", "soilassessment")
sapply(pkg, require, character.only = TRUE)

test_sf <- s_mps_sf[s_mps_sf$pedon_key %in% clhs_pkey, ]
test_sp <- as(test_sf, "Spatial")
vars <- c(ph_030 = "D:geodata/project_data/gsp-sas/ph_030_qrf_mlra_pm_NULL_final.tif",
          ph_100 = "D:geodata/project_data/gsp-sas/ph_100_qrf_mlra_pm_NULL_final.tif")
rs <- stack(vars)
rs <- as(rs, "SpatialGridDataFrame")

featureRep(rs["ph_030"], test_sp)
featureRep(rs["ph_100"], test_sp)

print(rbind(
  raster = summary(rs$ph_030)[1:6], 
  test   = summary(test2$ph_030_p)
  ), digits = 2)

print(rbind(
  raster = summary(rs$ph_100)[1:6], 
  test   = summary(test2$ph_100_p)
  ), digits = 2)


# bootstrap prediction interval
library(caret)
source("https://raw.githubusercontent.com/ncss-tech/gsp-sas/master/predUncertain.R")

vars <- ph_rf_030$finalModel$forest$independent.variable.names
idx <- which(names(rs) %in% vars)
rs2 <- rs[[idx]]

data2 <- as(s_mps_sf["ph_ptf.000_030_cm"], "Spatial")
data2 <- data2[!is.na(data2$ph_ptf.000_030_cm), ]
ex    <- extract(rs2, data2)
idx   <- complete.cases(as.data.frame(ex))
data2 <- data2[idx, ]

ph_pred_uncert <- predUncertain(data2, rs2, 3, 95, "ranger")


```



# SSURGO vs DSM accuracy

```{r}

load(file = "C:/Users/stephen.roecker/OneDrive - USDA/projects/gsp-sas/training_data.RData")
load(file = "C:/Users/stephen.roecker/OneDrive - USDA/projects/gsp-sas/ph_rf.RData")
load(file = "C:/Users/stephen.roecker/OneDrive - USDA/projects/gsp-sas/LDM-compact_20200709.RData")

stratia_f <- function(measured, predicted) {
  data.frame(
    Bias     =      mean( predicted - measured,    na.rm = TRUE),
    RMSE     = sqrt(mean((predicted - measured)^2, na.rm = TRUE)),
    Rsquared =      cor(  predicted,  measured, use = "complete.obs")^2,
    NSE      = 1 -  sum(  predicted - measured, na.rm = TRUE)^2 /
      sum((predicted - mean(measured, na.rm = TRUE))^2, na.rm = TRUE)
  )}

# SSURGO vs KSSL at 1000-meters
library(raster)
library(sf)

ssurgo_fn <- c(
  ec_030  = "gnatsgo_fy20_1km_ec_r.i0-30_wgs84.tif",
  ec_100  = "gnatsgo_fy20_1km_ec_r.i30-100_wgs84.tif",
  ph_030  = "ssurgo_ptph30_rs_med_v2.tif",
  ph_100  = "ssurgo_ptph100_rs_med_v2.tif",
  sar_030 = "ssurgo_ptsar30_rs_med.tif",
  sar_100 = "ssurgo_ptsar100_rs_med.tif"
  )
ssurgo_rs <- stack(file.path("D:/geodata/project_data/gsp-sas/1km covariates/SSURGO", ssurgo_fn))
names(ssurgo_rs) <- names(ssurgo_fn)

ssurgo_ex <- extract(ssurgo_rs, as(s_mps_sf, "Spatial"), df = TRUE, sp = TRUE)@data

ssurgo_stratia <- rbind(
  cbind(property = "EC", interval = "0-30",    stratia_f(ssurgo_ex$ec_ptf.000_030_cm,  ssurgo_ex$ec_030)),
  cbind(property = "EC", interval = "30-100",  stratia_f(ssurgo_ex$ec_ptf.030_100_cm,  ssurgo_ex$ec_100)),
  cbind(property = "pH", interval = "0-30",    stratia_f(ssurgo_ex$ph_ptf.000_030_cm,  ssurgo_ex$ph_030)),
  cbind(property = "pH", interval = "30-100",  stratia_f(ssurgo_ex$ph_ptf.030_100_cm,  ssurgo_ex$ph_100)),
  cbind(property = "SAR", interval = "0-30",   stratia_f(ssurgo_ex$sar.000_030_cm,     ssurgo_ex$sar_030)),
  cbind(property = "SAR", interval = "30-100", stratia_f(ssurgo_ex$sar.030_100_cm,     ssurgo_ex$sar_100))
)


dsm_fn <- c(
  ec_030  = "ec030_final_fill_mask.tif",
  ec_100  = "ec100_final_fill_mask.tif",
  ph_030  = "ph_030_qrf_mlra_pm_NULL_final4_fill_mask.tif",
  ph_100  = "ph_100_qrf_mlra_pm_NULL_final4_fill_mask.tif",
  esp_030 = "notr_esp_030_t25_nomlra_nafw.tif",
  esp_100 = "notr_esp_100_t25_nomlra_nafw.tif"
  )
ssurgo_fn <- c(
  ssurgo_ph_030 = "ssurgo_ptph30_rs_med.tif"
)

dsm_rs <- stack(file.path("D:/geodata/project_data/gsp-sas/predictions_v2", dsm_fn))
names(dsm_rs) <- names(dsm_fn)
ssurgo_r <- raster(file.path("D:/geodata/project_data/gsp-sas/1km covariates/SSURGO", ssurgo_fn))

dsm_ex <- extract(dsm_rs, as(s_mps_sf, "Spatial"), df = TRUE, sp = TRUE)@data
ssurgo_ex <- extract(ssurgo_r, as(s_mps_sf, "Spatial"))
dsm_ex <- cbind(dsm_ex, ssurgo_ph_030 = ssurgo_ex)

dsm_stratia <- rbind(
  cbind(property = "EC",  interval = "0-30",   stratia_f(dsm_ex$ec_ptf.000_030_cm,  dsm_ex$ec_030)),
  cbind(property = "EC",  interval = "30-100", stratia_f(dsm_ex$ec_ptf.030_100_cm,  dsm_ex$ec_100)),
  cbind(property = "pH",  interval = "0-30",   stratia_f(dsm_ex$ph_ptf.000_030_cm,  dsm_ex$ph_030)),
  cbind(property = "pH",  interval = "30-100", stratia_f(dsm_ex$ph_ptf.030_100_cm,  dsm_ex$ph_100)),
  cbind(property = "ESP", interval = "0-30",   stratia_f(dsm_ex$esp.000_030_cm, dsm_ex$esp_030)),
  cbind(property = "ESP", interval = "30-100", stratia_f(dsm_ex$esp.030_100_cm, dsm_ex$esp_100))
)


stratia <- rbind(
  cbind(source = "gNATSGO", ssurgo_stratia),
  cbind(source = "DSM",    dsm_stratia)
)
stratia[4:7] <- round(stratia[4:7], 2)
stratia$interval <- factor(stratia$interval, levels = c("0-30", "30-100")[2:1])
stratia <- dplyr::arrange(stratia, rev(interval), property)

vars <- c("Bias", "RMSE", "Rsquared", "NSE")
test <- reshape(stratia,
                direction = "long",
                timevar = "metric", times = vars,
                v.names = "value",  varying = vars
                )
ggplot(test, aes(x = value, y = interval, shape = source)) + 
  geom_point() + 
  facet_grid(property ~ metric, scales = "free_x") #, ncol = 4)

dim <- c(min(dsm_ex$ph_ptf.000_030_cm, na.rm = TRUE), 
         max(dsm_ex$ph_ptf.000_030_cm, na.rm = TRUE)
)

dsm_gg <- ggplot(dsm_ex, aes(x = ph_ptf.000_030_cm, y = ph_030)) + 
  geom_hex(show.legend = FALSE) + 
  geom_abline() +
  scale_fill_viridis_c() +
  xlim(dim) + ylim(dim) +
  coord_fixed() +
  xlab("pH measured") + ylab("pH predicted") +
  ggtitle("DSM")
  

ssurgo_gg <- ggplot(dsm_ex, aes(x = ph_ptf.000_030_cm, y = ssurgo_ph_030)) + 
  geom_hex(show.legend = FALSE) + 
  geom_abline() +
  scale_fill_viridis_c() +
  xlim(dim) + ylim(dim) +
  coord_fixed() +
  xlab("pH measured") + ylab("pH predicted") +
  ggtitle("SSURGO")

gridExtra::grid.arrange(dsm_gg, ssurgo_gg, ncol = 2) # , top = grid::textGrob("Predictions vs KSSL", gp = grid::gpar(fontsize = 15)))

```