lab_data.Rmd

---
title: "Lab Data - GSP"
author: "Stephen Roecker"
date: "July 9, 2020"
output: html_document
editor_options: 
  chunk_output_type: console
---



# LDM Sqlite

```{r ldm}

# https://new.cloudvault.usda.gov/index.php/s/eSoPYbWDBQNX2HP

library(DBI)

con <- dbConnect(RSQLite::SQLite(), "C:/Users/stephen.roecker/Nextcloud/data/ldm/LDM-compact_20200709.sqlite")
# area <- read.csv("C:/Users/stephen.roecker/Nextcloud/data/ldm/lab_area.txt", stringsAsFactors = FALSE)
# dbCreateTable(con, "area", area)
(ldm_names <- dbListTables(con))
ldm <- lapply(ldm_names, function(x) dbReadTable(con , x))
names(ldm) <- ldm_names
dbDisconnect(con)



# horizon table
chem_vars <- c("labsampnum", "ec_predict_one_to_two", "ec_predict_one_to_two_method", "electrical_conductivity_satx", "electrical_cond_satx_method", "ph_h2o", "ph_h2o_method", "ph_saturated_paste", "exchangeable_sodium", "sodium_absorption_ratio", "total_estimated_salts_satx")
phys_vars <- c("labsampnum", "texture_lab", "clay_total", "silt_total", "sand_total", "total_frag_wt_pct_gt_2_mm_ws", "bulk_density_third_bar", "estimated_organic_matter")
l_vars <- c("labsampnum", "layer_key", "pedon_key", "hzn_desgn", "hzn_top", "hzn_bot", "texture_description", "stratified_textures_flag")

h <- merge(ldm$layer[l_vars],
           ldm$physical[! duplicated(ldm$physical$labsampnum), phys_vars],
           by = "labsampnum",
           all.x = TRUE
           )
h <- merge(h, 
           ldm$chemical[chem_vars],
           by = "labsampnum",
           all.x = TRUE
           )



# site table
ncss_vars <- c("site_key", "pedon_key", "pedlabsampnum", "pedoniid", "samp_name", "corr_name", "samp_classdate", "corr_classdate")
site_vars <- c("site_key", "user_site_id", "latitude_std_decimal_degrees", "longitude_std_decimal_degrees")

s <- merge(ldm$nasis_ncss[ncss_vars],
           ldm$nasis_site[site_vars],
           by = "site_key",
           all.x = TRUE
           )
s <- within(s, {
            samp_classdate = strptime(samp_classdate, format = "%Y-%m-%d %H:%M:%S")
            corr_classdate = strptime(corr_classdate, format="%Y-%m-%d %H:%M:%S")
            })
s <- s[with(s, order(corr_classdate, samp_classdate, decreasing = TRUE)), ]



# soil profile collection object
# optionally test for bad horizonation... flag, and remove
  if (TRUE) {
    h.test <- plyr::ddply(h, 'pedon_key', function(d) {
      res <- aqp::hzDepthTests(top=d[['hzn_top']], bottom=d[['hzn_bot']])
      return(all(!res))
    })
    names(h.test)[2] <- "hz_logic_pass"

    # which are the good (valid) ones?
    good.ids <- as.character(h.test$pedon_key[which(h.test$hz_logic_pass)])
    bad.ids  <- as.character(h.test$pedon_key[which(! h.test$hz_logic_pass)])

    # keep the good ones
    h <- h[which(h$pedon_key %in% good.ids), ]

    # keep track of those components with horizonation errors
    if (length(bad.ids) > 0)
      assign('kssl.hz.problems', value=bad.ids)
  }


library(aqp)
h <- h[h$pedon %in% s$pedon_key, ]
s <- s[s$pedon_key %in% h$pedon_key, ]
spc <- h
depths(spc) <- pedon_key ~ hzn_top + hzn_bot
site(spc) <- s


# save(ldm, spc, file = "C:/Users/stephen.roecker/Nextcloud/data/ldm/LDM-compact_20200709.RData")

```



# pedotransfer functions

```{r pt}

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

library(aqp)
library(dplyr)

s <- as.data.frame(s_mps_sf)
h <- horizons(spc)
h <- within(h, {
  hzn_mid = (hzn_bot - hzn_top) /  2
  texture = factor(texture_lab, levels = aqp::SoilTextureLevels())
  hdm = NA
  hds = NA
  hdm[grepl("C", hzn_desgn)] = "C"
  hdm[grepl("B", hzn_desgn)] = "B"
  hdm[grepl("A", hzn_desgn)] = "A"
  hdm[grepl("O", hzn_desgn)] = "O"
  hds[grepl("k", hzn_desgn)] = "k"
  hds[grepl("kk", hzn_desgn)] = "kk"
  hds[grepl("k", hzn_desgn)] = "k"
  hds[grepl("y", hzn_desgn)] = "y"
  hds[grepl("z", hzn_desgn)] = "z"
  hds[grepl("n", hzn_desgn)] = "n"
})



# calculate weighted averages
h_seg <- segment(h, intervals = c(0, 30, 100), hzdepcols = c("hzn_top", "hzn_bot"))

test <- group_by(h_seg, pedon_key, segment_id) %>%
    summarize(clay = weighted.mean(clay_total, w = hzn_mid, na.rm = TRUE),
              ec  = weighted.mean(ec_predict_one_to_two, w = hzn_mid, na.rm = TRUE),
              ph  = log10(weighted.mean(1/10^-ph_h2o, w = hzn_mid, na.rm = TRUE)),
              sar = weighted.mean(sodium_absorption_ratio, w = hzn_mid, na.rm = TRUE),
              esp = weighted.mean(exchangeable_sodium, w = hzn_mid, na.rm = TRUE)
            ) %>%
  ungroup() %>%
  filter(!is.na(segment_id)) %>%
  as.data.frame()
test$segment_id <- paste0("i", test$segment_id)

test_w <- reshape(test,
                  direction = "wide",
                  idvar     = "pedon_key",
                  timevar   = "segment_id",
                  v.names   = c("clay", "ec", "ph", "sar", "esp")
                  )
attributes(test_w)$reshapeWide <- NULL



# sample
s <- test_w
vars <- c("pedon_key", "labsampnum", "hzn_desgn", "hzn_top", "hzn_bot", "texture_lab", "clay_total", "ec_predict_one_to_two", "electrical_conductivity_satx", "ph_h2o", "ph_saturated_paste")

set.seed(111)
idx <- s$pedon_key %in% unique(h[complete.cases(h[vars]), ]$pedon_key)
s_sub <- s[idx, ]
n    <- round(nrow(s_sub) * 0.75)
clhs_idx  <- clhs::clhs(s_sub[-1], size = n)
clhs_pkey <- s_sub[clhs_idx, "pedon_key"]


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



# model
# pH
library(rms)

ph_ols <- ols(ph_saturated_paste ~ rcs(ph_h2o, parms = c(0, 5.4, 6.9, 10.7)), data = train)
latex(ph_ols)

ph_lm <- lm(ph_saturated_paste ~ splines::ns(ph_h2o, 3), data = train, y = TRUE)
plot(ph_lm)
termplot(ph_lm, partial.resid = TRUE)
visreg::visreg(ph_lm)
summary(ph_lm)
plot(ph_lm$y, ph_lm$fitted.values); abline(0, 1)

ph_test <- predict(ph_lm, newdata = test) 
caret::R2(ph_test, test$ph_saturated_paste, na.rm = TRUE)
caret::RMSE(ph_test, test$ph_saturated_paste, na.rm = TRUE)


# EC
train <- train[row.names(train) != "288352", ]
train$ec_t <- log(train$ec_predict_one_to_two + 0.1)


ec_lrm <- lrm(electrical_conductivity_satx ~ rcs(log(ec_predict_one_to_two + 0.1), 3) + clay_total, data = train)

ec_glm <- glm(electrical_conductivity_satx ~ splines::ns(log(ec_predict_one_to_two + 0.1), 3) + clay_total, data = train, y = TRUE, family = Gamma(link = log))

ec_glm2 <- glm(electrical_conductivity_satx ~ splines::ns(log(ec_predict_one_to_two + 0.1), 3), data = train, y = TRUE, family = Gamma(link = log))



# glm 1
plot(ec_glm)
termplot(ec_glm, partial.resid = TRUE)
visreg::visreg(ec_glm, scale = "response", partial = TRUE)
with(ec_glm, plot(log(y + 0.1), log(fitted.values + 0.1))); abline(0, 1)
modEvA::Dsquared(ec_glm, adjust = TRUE)
with(ec_glm, caret::R2(fitted.values, y, na.rm = TRUE))
summary(ec_glm)

test$ec_t <- log(test$ec_predict_one_to_two + 0.1)
ec_test <- predict(ec_glm, newdata = test, type = "response") 
caret::R2(ec_test, test$electrical_conductivity_satx, na.rm = TRUE)
caret::RMSE(ec_test, test$electrical_conductivity_satx, na.rm = TRUE)
plot(log(test$electrical_conductivity_satx + 0.1), log(ec_test + 0.1)); abline(0, 1)


# glm 2
plot(ec_glm2)
termplot(ec_glm2, partial.resid = TRUE)
visreg::visreg(ec_glm2, scale = "response", partial = TRUE)
with(ec_glm, plot(log(y + 0.1), log(fitted.values + 0.1))); abline(0, 1)
modEvA::Dsquared(ec_glm, adjust = TRUE)
with(ec_glm2, caret::R2(fitted.values, y, na.rm = TRUE))
summary(ec_glm2)

test$ec_t <- log(test$ec_predict_one_to_two + 0.1)
ec_test <- predict(ec_glm2, newdata = test, type = "response") 
caret::R2(ec_test, test$electrical_conductivity_satx, na.rm = TRUE)
caret::RMSE(ec_test, test$electrical_conductivity_satx, na.rm = TRUE)
plot(log(test$electrical_conductivity_satx + 0.1), log(ec_test + 0.1)); abline(0, 1)


# rf
ec_rf <- randomForest(electrical_conductivity_satx ~ ec_predict_one_to_two + clay_total, data = train, na.action = na.exclude)
caret::R2(ec_rf$predicted, ec_rf$y, na.rm = TRUE)
caret::RMSE(ec_rf$predicted, ec_rf$y, na.rm = TRUE)
plot(log(ec_rf$y + 0.1), log(ec_rf$predicted + 0.1)); abline(0, 1)
visreg::visreg(ec_rf)


ec_test_rf <- predict(ec_rf, newdata = test)
caret::R2(ec_test_rf, test$electrical_conductivity_satx, na.rm = TRUE)
caret::RMSE(ec_test_rf, test$electrical_conductivity_satx, na.rm = TRUE)
plot(log(ec_test_rf + 0.1), log(test$electrical_conductivity_satx + 0.1)); abline(0, 1)


# soilassessment pedotransfer functions
library(soilassessment)

h <- within(h, {
  hzn_top     = as.integer(hzn_top)
  hzn_bot     = as.integer(hzn_bot)
  texture_lab = tolower(texture_lab)
  tex = NA
  tex[texture_lab %in% "c"]    = 1
  tex[texture_lab %in% "sic"]  = 2
  tex[texture_lab %in% "sc"]   = 8
  tex[texture_lab %in% "cl"]   = 7
  tex[texture_lab %in% "sicl"] = 3
  tex[texture_lab %in% "scl"]  = 9
  tex[texture_lab %in% "l"]    = 11
  tex[texture_lab %in% "sil"]  = 4
  tex[texture_lab %in% c("sl", "vfsl", "fsl")] = 5
  tex[texture_lab %in% c("ls", "vfls", "fls")] = 10
  tex[texture_lab %in% "s"]    = 12
  tex[texture_lab %in% "cos"]  = 13
  tex[texture_lab %in% c("fs", "vfs")] = 15
  ec_ptf = ECconversion1(ec = ec_predict_one_to_two, oc = estimated_organic_matter, clay = clay_total, texture = tex, soilsolution = "1:2", method = "sonmez")
  # ec_ptf = ifelse(is.na(electrical_conductivity_satx), ec_ptf, electrical_conductivity_satx)
})
plot(log(h$electrical_conductivity_satx + 0.1), log(h$ec_ptf + 0.1)); abline(0, 1)
caret::R2(h$electrical_conductivity_satx, h$ec_ptf, na.rm = TRUE)


# apply pedotransfer functions
h$ph_ptf <- predict(ph_lm, newdata = h)
h$ec_ptf <- predict(ec_glm, newdata = h, type = "response")
h <- within(h, {
  ph_ptf = ifelse(is.na(ph_saturated_paste), ph_ptf, ph_saturated_paste)
  ec_ptf = ifelse(is.na(electrical_conductivity_satx), ec_ptf, electrical_conductivity_satx)
  sar    = sodium_absorption_ratio
  esp    = exchangeable_sodium
  tss    = total_estimated_salts_satx
})

replaceHorizons(spc) <- h

# save(ldm, spc, clhs_pkey, ec_glm, ec_glm2, ph_lm, file = "C:/Users/stephen.roecker/Nextcloud/data/ldm/LDM-compact_202000923.RData")

```



# mpspline & weighted averages

```{r mps}

# mpspline
library(mpspline2)


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


gsp_vars <- c("ph_ptf", "ec_ptf", "sar", "esp", "tss")
s_mps_l <- lapply(gsp_vars, function(x) {
  
  vars <- c("pedon_key", "hzn_top", "hzn_bot", x)
  h2 <- horizons(spc)[vars]
  h2 <- h2[! is.na(h2[, x]), ]
  
  mps <- mpspline2::mpspline(
    h2,
    var_name = x, 
    lam = 0.8, 
    d   = c(0, 30, 100, 150)
    )
  mps_df <- lapply(mps, function(x2) {
    # d <- c("000_030_cm", "030_100_cm", "100_150_cm")
    data.frame(pedon_key = x2$pedon_key, var = x2$est_dcm[1:3]) ->.;
    names(.)[2] <- x
    .$depths <- names(x2$est_dcm)[1:3]
    return(.)
    })
  mps_df <- do.call("rbind", mps_df)
  
  # tidyr::spread(mps_df, depths, var)
  var <- names(mps_df)[2]
  mps_df <- reshape(mps_df,
                    direction = "wide",
                    idvar   = "pedon_key",
                    v.names = var,
                    timevar = "depths"
                    )
  attributes(mps_df)$reshapeWide <- NULL
  
  return(mps_df)
  })
names(s_mps_l) <- gsp_vars


# weighted average
s_wa <- aqp::segment(horizons(spc), intervals = c(0, 30, 100, 150), hzdepcols = c("hzn_top", "hzn_bot")) %>%
  mutate(hzn_bot = ifelse(is.na(hzn_bot), hzn_top + 1, hzn_bot),
         hzn_mid = hzn_bot - hzn_top,
         pedon_key = as.character(pedon_key)
         ) %>%
  # values come back NA if weights (w) are missing
  group_by(pedon_key, segment_id) %>%
  summarize(ec_wa = weighted.mean(ec_ptf, w = hzn_mid, na.rm = TRUE),
            ph_wa = log10(weighted.mean(1/10^-ph_ptf, w = hzn_mid, na.rm = TRUE)),
            sar_wa = weighted.mean(sar, w = hzn_mid, na.rm = TRUE),
            tss_wa = weighted.mean(tss, w = hzn_mid, na.rm = TRUE)
            ) %>%
  ungroup() %>%
  right_join(select(aqp::site(spc), pedon_key), by = "pedon_key") %>% 
  # reshape() doesn't work with tibbles
  as.data.frame()
s_wa_w <- reshape(data = s_wa, direction = "wide", 
                   idvar = "pedon_key", 
                   v.names = c("ph_wa", "ec_wa", "sar_wa", "tss_wa"), 
                   timevar = "segment_id"
                   )
attributes(s_wa_w)$reshapeWide <- NULL
s_wa_w <- s_wa_w[!grepl(".NA$", names(s_wa_w))]


# merge mpspline results
s <- site(spc)
s_mps <- merge(s,     s_mps_l$ph_ptf, by = "pedon_key", all.x = TRUE, sort = FALSE)
s_mps <- merge(s_mps, s_mps_l$ec_ptf, by = "pedon_key", all.x = TRUE, sort = FALSE)
s_mps <- merge(s_mps, s_mps_l$esp, by = "pedon_key", all.x = TRUE, sort = FALSE)
s_mps <- merge(s_mps, s_mps_l$sar, by = "pedon_key", all.x = TRUE, sort = FALSE)
s_mps <- merge(s_mps, s_mps_l$tss, by = "pedon_key", all.x = TRUE, sort = FALSE)
s_mps <- merge(s_mps, s_wa_w, by = "pedon_key", all.x = TRUE, sort = FALSE)


# sf object
library(sf)
s_mps_sf <- subset(s_mps, complete.cases(latitude_std_decimal_degrees, longitude_std_decimal_degrees))
s_mps_sf <- st_as_sf(s_mps_sf, 
                 coords = c("longitude_std_decimal_degrees", "latitude_std_decimal_degrees"),
                 crs = 4326
                 )

idx  <- grepl("^ph_ptf|^ec_ptf|^esp", names(s_mps_sf))
s_mps_sf2 <- s_mps_sf
s_mps_sf2[idx] <- lapply(st_drop_geometry(s_mps_sf2[idx]), function(x) {ifelse(is.na(x), - 9999, x)})
write_sf(s_mps_sf2, dsn = "D:/geodata/project_data/gsp-sas/s_mps_sf.shp", layer = "ldm")

# save(ldm, spc, clhs_pkey, ph_lm, ec_glm, s_mps_l, s_wa_w, s_mps_sf, file = "C:/Users/stephen.roecker/Nextcloud/projects/2020_gsp-sas/LDM-compact_20200709.RData")

```


# Map

```{r}
library(ggmap)
library(sf)
library(ggplot2)
library(USAboundaries)

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

st <- us_states() 
st <- subset(st, !state_name %in% c("Alaska", "Hawaii", "Puerto Rico"))
st <- st_transform(st, crs = 5070)

# conus <- read_sf(dsn = "D:/geodata/soils/SSURGO_CONUS_FY19.gdb", layer = "SAPOLYGON")
# conus <- st_transform(conus, crs = 4326)
# conus$state <- substr(conus$AREASYMBOL, 1, 2)
# # conus <- conus %>% group_by(state) %>% summarize(state = state[1])
# conus <- aggregate(conus, list(conus$state))
# bb    <- st_bbox(conus)
# bb    <- make_bbox(lon = bb[c(3, 1)], lat = bb[c(2, 4)])
# gmap  <- get_map(bb, maptype = "terrain", source = "osm")

s   <- s_mps_sf
s   <- st_transform(s_mps_sf, crs = 5070)
idx <- st_intersects(s, st) 
idx <- sapply(idx, function(x) length(x) > 0)
s   <- s[idx, ]
s   <- within(s, {
  property = NA
  property = ifelse(!is.na(ph_ptf.000_030_cm), "pH", property)
  property = ifelse(complete.cases(ph_ptf.000_030_cm, s$esp.000_030_cm), "pH & ESP", property)
  property = ifelse(complete.cases(ph_ptf.000_030_cm, esp.000_030_cm, ec_ptf.000_030_cm), "pH & ESP & EC", property)
})
s$property <- factor(s$property, ordered = TRUE)
s <- subset(s, !is.na(property))
s <- s[order(s$property), ]
# s_t <- as.data.frame(cbind(st_coordinates(s), pH = s$ph_ptf.000_030_cm))
# 
# brks <- c(0, 3.5, 4.5, 5.1, 5.6, 6.1, 6.6, 7.4, 7.9, 8.5, 9.5, 14)
# s_t$ph_lev <- cut(s_t$pH, 
#                        breaks = brks,
#                        labels = paste(brks[-length(brks)], brks[-1], sep = "-")
#                        )

gg_gsp <- ggplot() +
  geom_sf(data = st, fill = NA) +
  geom_sf(data = s, aes(col = property), size = 0.2) +
  # scale_color_manual(values = c("blue", "orange")) +
  guides(color = guide_legend(override.aes = list(size = 2))) +
  scale_fill_viridis_d() +
  # ggtitle("Pedon locations") +
  theme_bw()
ggsave(gg_gsp, file = "test.png", dpi = 300)
  
# geom_hex(data = s_t, aes(x = X, y = Y, fill = pH)) + 
  # scale_fill_viridis_c()

library(tmap)

tm_shape(st) + tm_borders() +
tm_shape(s_pH) + tm_dots(col = "blue")


```