ssurgo_derivatives.Rmd

---
title: "SSURGO derivatives"
author: "Stephen Roecker"
date: "7/14/2020"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}

knitr::opts_chunk$set(echo = TRUE)

# load pedotransfer functions
load(file = "C:/Users/stephen.roecker/Box/data/LDM-compact_20200709.RData")

```



# MLRA overlaps

```{r mlra}
library(dplyr)

nmo <- read.csv("C:/Users/stephen.roecker/OneDrive - USDA/data/mukey_mlra_overlap.csv", stringsAsFactors = FALSE)
nmo$mukey <- as.character(nmo$mukey)
nmo <- nmo[with(nmo, order(mukey, - membership)), ]
# nmo <- nmo[! duplicated(nmo$mukey), ]
nmo <- nmo %>%
  group_by(mukey) %>%
  summarize(muacres = sum(area_ac),
            membership = membership[which.max(membership)],
            mlra    = mlra[which.max(membership)]
            ) %>%
  ungroup() %>%
  as.data.frame()

```



# fetchGDB

```{r fetchGDB}
library(aqp)
library(soilDB)

f_gnatsgo <- fetchGDB(dsn ="G:/archive/gNATSGO_CONUS_FY20.gdb", WHERE = "areasymbol != 'US'")
mu_us <- get_mapunit_from_GDB(dsn = "G:/archive/gNATSGO_CONUS_FY20.gdb", stats = TRUE)

f_statsgo <- fetchGDB(dsn = "D:/geodata/soils/wss_gsmsoil_US_20161013.gdb", WHERE = "areasymbol = 'US'", childs = TRUE)

f_us2 <- f_natsgo[! f_natsgo$mukey %in% f_statsgo$mukey]
f_us  <- aqp::combine(list(f_us2, f_statsgo))

missing_mukeys <- mu_us[is.na(mu_us$n_component), "mukey"]

# save(f_us, mu_us, f_statsgo, file = "C:/Users/stephen.roecker/OneDrive - USDA/projects/gsp-sas/gnatsgo_fy20.RData")
load(file = "C:/Users/stephen.roecker/OneDrive - USDA/projects/gsp-sas/gnatsgo_fy20.RData")
# load(file = "C:/Users/stephen.roecker/OneDrive - USDA/projects/gsp-sas/gsp_salinity_fetchGDB.RData")


# missing mukeys
# f_sda <- fetchSDA(WHERE = paste0("mukey IN ('", paste0(missing_mukeys, collapse = "', '"), "')"), duplicates = TRUE)
# s <- site(f_sda)[siteNames(f_us)]
# h <- horizons(f_sda)[horizonNames(f_us)]
# f_mis <- h
# depths(f_mis) <- cokey ~ hzdept_r + hzdepb_r
# site(f_mis) <- s
# 
# f_us <- aqp::combine(list(f_us, f_mis))


# generate unique (?) mukeys (aka nationalmusym)
# s     <- site(f_us)
# s     <- merge(s, mu_us[c("mukey", "muname")], by = "mukey", all.x = TRUE)
# s$idx <- with(s, paste(muname, compname, comppct_r, compkind, majcompflag, localphase, drainagecl, erocl, slope_r))
# s_sub <- s[! duplicated(s$idx), ]
# horizons(f_us) <- s[c("cokey", "idx")]
# idx <- which(! duplicated(horizons(f_us)$idx))
# f_us_sub <- f_us[, idx]

```



# Aggregate chorizon data to map units

```{r}

library(aqp)
library(dplyr)

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

co_agg <- aqp::segment(aqp::horizons(f_us), intervals = c(0, 30, 100), hzdepcols = c("hzdept_r", "hzdepb_r")) %>%
  mutate(ec_predict_one_to_two = ec_r,
         ph_h2o     = ph1to1h2o_r,
         clay_total = claytotal_r,
         )
# pedotransfer functions
co_agg$ph_ptf <- predict(ph_lm,  co_agg)
# co_agg$ec_ptf <- predict(ec_glm, co_agg, type = "response")
# aggregate components
co_agg <- co_agg %>%
  mutate(hzdepb_r = ifelse(is.na(hzdepb_r), hzdept_r + 1, hzdepb_r),
         hzdepm_r = hzdepb_r - hzdept_r
         ) %>%
  # values come back NA if weights (w) are missing
  group_by(cokey, segment_id) %>%
  summarize(ec_r  = weighted.mean(ec_r, w = hzdepm_r, na.rm = TRUE),
            ph_r  = log10(weighted.mean(1/10^-ph_ptf, w = hzdepm_r, na.rm = TRUE)),
            sar_r = weighted.mean(sar_r, w = hzdepm_r, na.rm = TRUE)
            ) %>%
  ungroup() %>%
  right_join(select(aqp::site(f_us), mukey, cokey, comppct_r), by = "cokey")


# aggregate mapunits
mu_agg_h <- co_agg %>%
  mutate(comppct_r = ifelse(is.na(comppct_r), 0, comppct_r)) %>%
  group_by(mukey, segment_id) %>%
  summarize(ec_r  = weighted.mean(ec_r, w = comppct_r, na.rm = TRUE),
            ph_r = log10(weighted.mean(1/10^-ph_r, w = comppct_r, na.rm = TRUE)),
            sar_r = weighted.mean(sar_r, w = comppct_r, na.rm = TRUE)
            ) %>%
  ungroup() %>%
  filter(!is.na(segment_id)) %>%
  as.data.frame()
mu_agg_h$segment_id <- paste0("i", mu_agg_h$segment_id)
mu_agg_h2 <- reshape(mu_agg_h,
                     direction = "wide",
                     idvar     = "mukey",
                     timevar   = "segment_id",
                     v.names   = c("ec_r", "ph_r", "sar_r")
                     )
attributes(mu_agg_h2)$reshapeWide <- NULL

# landform
mu_agg_landform <- aqp::site(f_us) %>%
  filter(!is.na(landform)) %>%
  group_by(mukey, landform) %>%
  summarize(comppct_r = sum(comppct_r, na.rm = TRUE)) %>%
  arrange(mukey, - comppct_r) %>%
  filter(!duplicated(mukey)) %>%
  select(- comppct_r)

# pmkind
# find components with missing pmkind
mis_mukeys <- aqp::site(f_us) %>% 
  group_by(mukey) %>% 
  summarize(na = all(is.na(pmkind))) %>% 
  ungroup() %>% 
  filter(na == TRUE)
mis_compname <- unique(f_us[f_us$mukey %in% mis_mukeys$mukey, ]$compname)
ss <- soilDB::get_soilseries_from_NASIS()
mis_compname <- mis_compname[tolower(mis_compname) %in% tolower(gsub("^ | $", "", ss$soilseriesname))]

# replace components with missing pmkind from same named components from gNATSGO
mis_compname_pmkind <- aqp::site(f_us) %>%
  mutate(compname = tolower(compname)) %>%
  # filter(compname %in% tolower(mis_compname)) %>%
  group_by(compname, pmkind) %>%
  summarize(comppct_r = sum(comppct_r, na.rm = TRUE)) %>%
  group_by(compname) %>%
  arrange(compname, - comppct_r) %>%
  filter(!duplicated(compname)) %>%
  select(- comppct_r)

# replace components with missing pmkind from same named components from gNATSGO
s_tmp <- aqp::site(f_us)
s_tmp$compname <- gsub("^ | $", "", tolower(s_tmp$compname))
mis_compname <- subset(s_tmp, 
                       compname %in% tolower(ss$soilseriesname)
                       & !compname %in% unique(s_tmp$compname[!is.na(s_tmp$pmkind)]) 
                       # & !compname %in% tolower(mis_compname_pmkind$compname)
                       )
osd <- lapply(unique(mis_compname$compname), function(x) soilDB::fetchOSD(x, extended = TRUE))
names(osd) <- unique(mis_compname$compname)
osd_pmkind <- lapply(osd, function(x) x$pmkind)
idx <- sapply(osd_pmkind, is.data.frame)
osd_pmkind <- do.call("rbind", lapply(osd_pmkind[idx], function(x) x$pmkind[1]))
osd_pmkind <- data.frame(compname = rownames(osd_pmkind), pmkind = osd_pmkind[, 1, drop = TRUE])


# compute dominant condition of parent material by mukey
s <- aqp::site(f_us) %>% 
  mutate(compname = tolower(compname)) %>%
  left_join(mis_compname_pmkind, by = "compname") %>%
  left_join(osd_pmkind, by = "compname") %>%
  mutate(pmkind.x = ifelse(is.na(pmkind.x), pmkind.y, pmkind.x),
         pmkind.x = ifelse(is.na(pmkind.x), pmkind,   pmkind.x),
         pmkind.y = NULL,
         pmkind   = NULL
         )
names(s)[names(s) == "pmkind.x"] <- "pmkind"

mu_agg_pmkind <- s %>%
  filter(!is.na(pmkind)) %>%
  group_by(mukey, pmkind) %>%
  summarize(comppct_r = sum(comppct_r, na.rm = TRUE)) %>%
  arrange(mukey, - comppct_r) %>%
  filter(!duplicated(mukey)) %>%
  select(- comppct_r)


# pmorigin
mu_agg_pmorigin <- site(f_us) %>%
  filter(!is.na(pmorigin)) %>%
  group_by(mukey, pmorigin) %>%
  summarize(comppct_r = sum(comppct_r, na.rm = TRUE)) %>%
  arrange(mukey, - comppct_r) %>%
  filter(!duplicated(mukey)) %>%
  select(- comppct_r)


# join
mu_agg <- mu_us["mukey"] %>%
  left_join(mu_agg_h2,       by = "mukey") %>%
  left_join(mu_agg_landform, by = "mukey") %>%
  left_join(mu_agg_pmorigin, by = "mukey") %>%
  left_join(mu_agg_pmkind, by = "mukey") %>%
  left_join(nmo, by = "mukey")


f_us$drainagecl2 <- factor(f_us$drainagecl, rev(levels(f_us$drainagecl)), ordered = TRUE)

mu_agg_pd_pct <- site(f_us) %>%
  group_by(mukey, drainagecl2) %>%
  summarize(comppct_r = sum(comppct_r, na.rm = TRUE)) %>%
  arrange(mukey, - comppct_r, drainagecl2) %>%
  ungroup() %>%
  group_by(mukey) %>%
  summarize(pd_pct = sum(comppct_r[drainagecl2 <= "Somewhat poorly drained"], na.rm = TRUE),
            drainagecl = drainagecl2[1]
            ) %>%
  ungroup()

write.csv(mu_agg_pd_pct, "mu_agg_drainagecl_fy20.csv", row.names = FALSE)

# save(f_us, f_statsgo, mu_us, mu_agg, osd, file = "gsp_salinity_fetchGDB_20221017.RData")

```



# Aggregate map units to landform, pmorigin, and pmkind

```{r}

setwd("C:/Users/stephen.roecker/OneDrive - USDA/projects/gsp-sas/")

load(file = "gsp_salinity_fetchGDB_20221017.RData")

library(dplyr)

# aggregate by pmkind, pmorigin and landform
mu_pmorigin <- mu_agg %>%
  # pmkind, pmorigin, and landform
  filter(muacres > 0) %>%
  group_by(mlra, pmorigin) %>%
  summarize(pct = ifelse(any(ec_r > 2), sum(muacres[ec_r > 2], na.rm = TRUE) / sum(muacres, na.rm = TRUE), 0)) %>%
  ungroup() %>%
  mutate(pct = round(pct * 100))
  # filter(ec_r > 2)


mu_pmkind <- mu_agg %>%
  # pmkind, pmorigin, and landform
  filter(muacres > 0) %>%
  group_by(mlra, pmkind) %>%
  summarize(p = ifelse(any(ec_r > 2), sum(muacres[ec_r > 2], na.rm = TRUE) / sum(muacres, na.rm = TRUE), 0)) %>%
  ungroup() %>%
  mutate(p = round(p * 100))
# filter(ec_r > 2)


mu_landform <- mu_agg %>%
  # pmkind, pmorigin, and landform
  filter(muacres > 0) %>%
  group_by(mlra, landform) %>%
  summarize(p = ifelse(any(ec_r > 2), sum(muacres[ec_r > 2], na.rm = TRUE) / sum(muacres, na.rm = TRUE), 0)) %>%
  ungroup() %>%
  mutate(p = round(p * 100))
# filter(ec_r > 2)


# mu_agg$ec_c <- cut(mu_agg$ec_r, breaks = c(0, 2, 4, 8, 16))
# ec_t <- as.data.frame(xtabs(~ pmkind + ec_c, data = mu))
# ec_t <- tidyr::spread(ec_t, ec_c, Freq)

```



# Generate SSURGO derivatives

```{r}

mu_agg <- within(mu_agg, {
  pmkind2 = as.character(pmkind)
  # Mine
  pmkind2[grepl("Human|mine spoil|Mine|Dredge spoils", pmkind2)] = "Mine spoil or earthy fill"
  pmkind2[grepl("Marine|Glaciomarine|Lagoonal deposits|Fluviomarine deposits|Estuarine deposits|Greensands", pmkind2)] = "Marine deposits"
  pmkind2[grepl("mixed", pmkind2)] = NA
  # Others
  pmkind2[grepl("Residuum|grus|igneous|Saprolite|sandstone|limestone|shale|conglomerate", pmkind2, ignore.case = TRUE)] = "Residuum"
  # Glacial
  pmkind2[grepl("till|Till|Drift", pmkind2)] = "Till"
  pmkind2[grepl("Glaciofluvial deposits|Outwash", pmkind2)] = "Outwash"
  pmkind2[grepl("Glaciolacustrine deposits", pmkind2)] = "Lacustrine deposits"
  # Colluvium
  pmkind2[grepl("Valley side alluvium|Slope alluvium|Creep deposits|Pedisediment|Cryoturbate|Colluvium", pmkind2)] = "Colluvium"
  pmkind2[grepl("spread|slide|Earthflow|flow deposits|Debris|Talus|rock|Scree|mass movement", pmkind2, ignore.case = TRUE)] = "Mass movement deposits"
  # Ash
  pmkind2[grepl("Pyroclastic|Tephra|Ash flow|volcanic ash|Lahar deposits|Pyroclastic flow|Tuff|Scoria|Cinders|Pumice", pmkind2, ignore.case = TRUE)] = "Volcanic ash"
  # Alluvium
  pmkind2[grepl("Overbank deposits|Alluvium", pmkind2)] = "Alluvium"
  # Lacustrine
  pmkind2[grepl("Backswamp deposit|Diatomaceous earth|Diamicton|Marl|Lacustrine", pmkind2)] = "Lacustrine deposits"
  # Other
  pmkind2[grepl("Loess", pmkind2, ignore.case = TRUE)] = "Loess"
  pmkind2[grepl("Eolian|Beach sand", pmkind2)] = "Eolian deposits"
  # Organic
  pmkind2[grepl("organic|Organic|Coprogenic material", pmkind2)] = "Organic material"
  })


# save(f_us, mu_us, mu_agg, osd, file = "C:/Users/stephen.roecker/OneDrive - USDA/projects/gsp-sas/gsp_salinity_fetchGDB_20220909.csv.RData")
write.csv(mu_agg, file = "C:/Users/stephen.roecker/OneDrive - USDA/projects/gsp-sas/gnatsgo_gsp_muaggatt_20221003.csv", row.names = FALSE)

```



# rastify

```{r}

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

mu_agg <- within(mu_agg, {
  mukey   = as.numeric(mukey)
  pmkind3 = as.integer(as.factor(pmkind2))
})
names(mu_agg)[1] <- "ID"

mlra <- "D:/geodata/project_data/gsp-sas/1km covariates/Other/mlra_rs.tif"
r <- raster("D:/geodata/project_data/gsp-sas/1km covariates/SSURGO/statsgo_pmkind3.tif")
gdalUtilities::gdalwarp(
    srcfile = mlra,
    dstfile = gsub(".tif", "_v2.tif", mlra),
    t_srs = proj4string(r), 
    te = bbox(r),
    r = "near",
    tr = res(r),
    ot = "Int32"
)


library(raster)

gnatsgo <- "D:/geodata/soils/gnatsgo_fy20_30m.tif"
r <- raster(gnatsgo)
crs(r) <- sf::st_crs(5070)$proj4string
gdalUtilities::gdalwarp(
    srcfile = gnatsgo,
    dstfile = gsub("30m.tif", "120m.tif", gnatsgo),
    t_srs = proj4string(r), 
    te = bbox(r),
    r = "mode",
    tr = c(120, 120),
    ot = "Int32",
    overwrite = TRUE
)

r <- raster("D:/geodata/soils/gnatsgo_fy20_30m.tif")
# r2 <- r[1:100, 1:100, drop = FALSE]
# r2 <- ratify(r2)
# rat <- levels(r2)[[1]]
levels(r) <- mu_agg

vars <- c("pmkind3", "ph_r.i0-30", "ph_r.i30-100", "ec_r.i0-30", "ec_r.i30-100")

# ratify
lapply(vars[4:5], function(x) {
  cat(x, as.character(Sys.time()), "\n")
  # beginCluster(type = "SOCK")
  deratify(r, att = x, 
           filename = paste0("D:/geodata/soils/gnatsgo_fy20_300m_", x, ".tif"),
           options = c("COMPRESS=DEFLATE"), 
           overwrite = TRUE, progress = "text" 
           )
  # endCluster()
  cat(x, as.character(Sys.time()), "\n")
  })


# upscale
lapply(vars[4:5], function(x) {
  cat(x, as.character(Sys.time()), "\n")
  # beginCluster(type = "SOCK")
  tmp <- paste0("D:/geodata/soils/gnatsgo_fy20_120m_", x, ".tif")
  
  gdalUtilities::gdalwarp(
    srcfile = tmp,
    dstfile = gsub("120m", "1km", tmp),
    t_srs = proj4string(r), 
    te = bbox(r),
    r = "average",
    tr = c(1000, 1000),
    overwrite = TRUE
    )
  # endCluster()
  cat(x, as.character(Sys.time()), "\n")
  })


# reproject
r <- raster("D:/geodata/project_data/gsp-sas/1km covariates/ISRIC/CONUS/SLPMRG5_CONUS.tif")

lapply(vars[4:5], function(x) {
  cat(x, as.character(Sys.time()), "\n")
  # beginCluster(type = "SOCK")
  tmp <- paste0("D:/geodata/soils/gnatsgo_fy20_1km_", x, ".tif")
  
  gdalUtilities::gdalwarp(
    srcfile = tmp,
    dstfile = gsub(".tif", "_wgs84.tif", tmp),
    t_srs = proj4string(r), 
    te = bbox(r),
    r = "bilinear",
    tr = res(r),
    overwrite = TRUE
    )
  # endCluster()
  cat(x, as.character(Sys.time()), "\n")
  })


# library(terra)
# mat <- matrix(mu_agg[c("ID", "ec_r.i0-30")])
# mat <- cbind(mat[[1]] - 0.5, mat[[1]] + 0.5, mat[[2]])
# classify(x = r, rcl = mat, filename = "D:/geodata/soils/gnatsgo_fy20_30m_ec000030_v2.tif", overwrite = TRUE, wopt = list(gdal = c("COMPRESS=DEFLATE")))


pmkind <- raster("D:/geodata/project_data/gsp-sas/1km covariates/SSURGO/ssurgo_pm_lu3_rs_mode_v2.tif")
pmkind2 <- projectRaster(pmkind, crs = "+init=epsg:5070", method = "ngb")



lv <- unique(mu_agg$pmkind2)
lv[3] <- "NA"


library(tmap)
tmap_pmkind <- tm_shape(pmkind2, projection = 5070, raster.downsample = T) + 
  tm_raster(palette = RColorBrewer::brewer.pal(13, name = "Paired"), labels = levels(as.factor(mu_agg$pmkind2)), breaks = 1:14, title = "Generalized Parent \n Material Kind") + 
  tm_legend(legend.outside = TRUE) +
  c
tmap_pmkind + png(filename = "gnatsgo_pmkind3.png", width = 8, height = 5, units = "in", res = 300)
dev.off()


```


# STASTGO

```{r}

library(raster)
library(terra)

# rasterize STATSGO
g_r <- rast("D:/geodata/project_data/gsp-sas/1km covariates/SSURGO/mlra_rs.tif")
names(g_r) <- "statsgo_mukey"

statsgo <- vect("D:/geodata/soils/wss_gsmsoil_US_[2016-10-13]/gsmsoilmu_a_us.shp")
idx <- as.integer(statsgo$MUKEY)
g_r <- rasterize(statsgo, g_r, field = idx)
writeRaster(g_r, filename = "D:/geodata/project_data/gsp-sas/1km covariates/SSURGO_move/statsgo_mukey.tif", wopt = list(datatype = "INT4U"), overwrite = TRUE)


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


# MLRA & STATSGO
gdalUtils::gdal_translate(
  src_dataset = "D:/geodata/project_data/gsp-sas/1km covariates/SSURGO/mlra_rs.tif", 
  dst_dataset = "mlra.sdat",
  ot = "UInt32",
  of = "SAGA"
  )
gdalUtils::gdal_translate(
  src_dataset = "D:/geodata/project_data/gsp-sas/1km covariates/SSURGO_move/statsgo_mukey.tif", 
  dst_dataset = "mukey.sdat",
  ot = "UInt32",
  of = "SAGA"
  )
# compute zonal stats in SAGA (it's fastest)
zs <- read.table("C:/Users/stephen.roecker/Nextcloud/projects/2020_gsp-sas/mukey_mlra_zonal.txt", header = TRUE, sep = "\t")
zs <- zs[order(zs$mukey, - zs$Count.UCU), ]
zs <- zs[! duplicated(zs$mukey), ]
names(zs)[1] <- "ID"

g_r <- raster("D:/geodata/project_data/gsp-sas/1km covariates/SSURGO_move/statsgo_mukey.tif")
levels(g_r) <- zs
raster::deratify(g_r, att = "mlra", 
                 filename = "D:/geodata/project_data/gsp-sas/1km covariates/SSURGO/statsgo_mlra.tif",
                 options = c("COMPRESS=DEFLATE"), 
                 overwrite = TRUE, progress = "text"
                 )



# SDA
library(soilDB)

mu <- get_mapunit_from_SDA(WHERE = "areasymbol = 'US'")
idx <- c(0, rep(1500, 10) * 1:10)
mu$idx <- as.character(cut(1:nrow(mu), breaks = idx))
f <- by(mu, mu$idx, function(x){
  f <- fetchSDA(WHERE = paste0("mukey IN ('", paste0(x$mukey, collapse = "', '"), "')"), childs = FALSE, duplicates = TRUE)
})
f <- lapply(f, function(x) x)
f <- aqp::combine(f)



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

co_agg <- aqp::segment(aqp::horizons(f), intervals = c(0, 30, 100), hzdepcols = c("hzdept_r", "hzdepb_r")) %>%
  mutate(ph_h2o     = ph1to1h2o_r,
         )
# pedotransfer functions
co_agg$ph_ptf <- predict(ph_lm,  co_agg)
# aggregate components
co_agg <- co_agg %>%
  mutate(hzdepb_r = ifelse(is.na(hzdepb_r), hzdept_r + 1, hzdepb_r),
         hzdepm_r = hzdepb_r - hzdept_r
         ) %>%
  # values come back NA if weights (w) are missing
  group_by(cokey, segment_id) %>%
  summarize(ph_r  = log10(weighted.mean(1/10^-ph_ptf, w = hzdepm_r, na.rm = TRUE)),
            ec_r  = weighted.mean(ec_r, w = hzdepm_r, na.rm = TRUE)
            ) %>%
  ungroup() %>%
  right_join(select(aqp::site(f), mukey, cokey, comppct_r), by = "cokey")


# aggregate mapunits
mu_agg_h <- co_agg %>%
  mutate(comppct_r = ifelse(is.na(comppct_r), 0, comppct_r)) %>%
  group_by(mukey, segment_id) %>%
  summarize(ph_r = log10(weighted.mean(1/10^-ph_r, w = comppct_r, na.rm = TRUE)),
            ec_r = weighted.mean(ec_r, w = comppct_r, na.rm = TRUE)
            ) %>%
  ungroup() %>%
  filter(!is.na(segment_id)) %>%
  as.data.frame()
mu_agg_h$segment_id <- paste0("i", mu_agg_h$segment_id)
mu_agg_h2 <- reshape(mu_agg_h,
                     direction = "wide",
                     idvar     = "mukey",
                     timevar   = "segment_id",
                     v.names   = c("ph_r", "ec_r")
                     )
attributes(mu_agg_h2)$reshapeWide <- NULL



# extract dominant parent material
mis_compname_pmkind$compname <- tolower(mis_compname_pmkind$compname)
mu_agg_pmkind <- co %>%
  mutate(compname = tolower(compname),
         mukey    = NULL
         ) %>%
  left_join(select(mu, mukey, nationalmusym), by = "nationalmusym") %>%
  left_join(select(mis_compname_pmkind, compname, pmkind), by = "compname") %>%
  filter(!is.na(pmkind)) %>%
  group_by(mukey, pmkind) %>%
  summarize(comppct_r = sum(comppct_r, na.rm = TRUE)) %>%
  arrange(mukey, - comppct_r) %>%
  filter(!duplicated(mukey)) %>%
  select(- comppct_r) %>%
  right_join(select(mu, mukey), by = "mukey") %>%
  as.data.frame()

mu_agg_pmkind <- within(mu_agg_pmkind, {
  pmkind2 = as.character(pmkind)
  pmkind2[is.na(pmkind2)] =  "z_missing"
  # Mine
  pmkind2[grepl("Human|mine spoil|Mine|Dredge spoils", pmkind2)] = "Mine spoil or earthy fill"
  pmkind2[grepl("Marine|Glaciomarine|Lagoonal deposits|Fluviomarine deposits|Estuarine deposits|Greensands", pmkind2)] = "Marine deposits"
  pmkind2[grepl("mixed", pmkind2)] = NA
  # Others
  pmkind2[grepl("Residuum|grus|igneous|Saprolite|sandstone|limestone|shale|conglomerate", pmkind2, ignore.case = TRUE)] = "Residuum"
  # Glacial
  pmkind2[grepl("till|Till|Drift", pmkind2)] = "Till"
  pmkind2[grepl("Glaciofluvial deposits|Outwash", pmkind2)] = "Outwash"
  pmkind2[grepl("Glaciolacustrine deposits", pmkind2)] = "Lacustrine deposits"
  # Colluvium
  pmkind2[grepl("Valley side alluvium|Slope alluvium|Creep deposits|Pedisediment|Cryoturbate|Colluvium", pmkind2)] = "Colluvium"
  pmkind2[grepl("spread|slide|Earthflow|flow deposits|Debris|Talus|rock|Scree|mass movement", pmkind2, ignore.case = TRUE)] = "Mass movement deposits"
  # Ash
  pmkind2[grepl("Pyroclastic|Tephra|Ash flow|volcanic ash|Lahar deposits|Pyroclastic flow|Tuff|Scoria|Cinders|Pumice", pmkind2, ignore.case = TRUE)] = "Volcanic ash"
  # Alluvium
  pmkind2[grepl("Overbank deposits|Alluvium", pmkind2)] = "Alluvium"
  # Lacustrine
  pmkind2[grepl("Backswamp deposit|Diatomaceous earth|Diamicton|Marl|Lacustrine", pmkind2)] = "Lacustrine deposits"
  # Other
  pmkind2[grepl("Loess", pmkind2, ignore.case = TRUE)] = "Loess"
  pmkind2[grepl("Eolian|Beach sand", pmkind2)] = "Eolian deposits"
  # Organic
  pmkind2[grepl("organic|Organic|Coprogenic material", pmkind2)] = "Organic material"
  })


mu_agg_pmkind <- within(mu_agg_pmkind, {
  mukey   = as.numeric(mukey)
  pmkind3 = as.integer(factor(pmkind2, levels = sort(unique(c(mu_agg$pmkind2, "z_missing")))))
})
names(mu_agg_pmkind)[1] <- "ID"



# deratify STATSGO
library(raster)

names(mu_agg_h2)[1] <- "ID"
mu_agg_h2[is.na(mu_agg_h2$`ph_r.i0-30`),   "ph_r.i0-30"] <- 7
mu_agg_h2[is.na(mu_agg_h2$`ph_r.i30-100`), "ph_r.i30-100"] <- 7
summary(mu_agg_h2)

g_r <- raster("D:/geodata/project_data/gsp-sas/1km covariates/SSURGO_move/statsgo_mukey.tif")
levels(g_r) <- mu_agg_h2

vars <- c("pmkind3", "ph_r.i0-30", "ph_r.i30-100", "ec_r.i0-30", "ec_r.i30-100")[4:5]
lapply(vars, function(x) {
  cat(x, as.character(Sys.time()), "\n")
  # beginCluster(type = "SOCK")
  tmp <- paste0("D:/geodata/project_data/gsp-sas/1km covariates/SSURGO/statsgo_", x, ".tif")
  
  raster::deratify(g_r, att = x, 
           filename = tmp,
           options = c("COMPRESS=DEFLATE"), 
           overwrite = TRUE, progress = "text" 
           )
  # endCluster()
  cat(x, as.character(Sys.time()), "\n")
  })


natsgo_r  <- raster("D:/geodata/project_data/gsp-sas/1km covariates/SSURGO/ssurgo_pm_lu3_rs_mode.tif")
statsgo_r <- raster("D:/geodata/project_data/gsp-sas/1km covariates/SSURGO/statsgo_pmkind3.tif")
pmkind <- merge(natsgo_r, statsgo_r, overlap = TRUE, filename = "D:/geodata/project_data/gsp-sas/1km covariates/SSURGO/pmkind3.tif", progress = "text", overwrite = TRUE)


natsgo_ph030_r  <- raster("D:/geodata/project_data/gsp-sas/1km covariates/SSURGO/ssurgo_ptph30_rs_med.tif")
statsgo_ph030_r <- raster("D:/geodata/project_data/gsp-sas/1km covariates/SSURGO/statsgo_ph_r.i0-30.tif")
merge(natsgo_ph030_r, statsgo_ph030_r, overlap = TRUE, filename = "D:/geodata/project_data/gsp-sas/1km covariates/SSURGO/ssurgo_ptph30_rs_med_v2.tif", progress = "text", overwrite = TRUE)


natsgo_ph100_r  <- raster("D:/geodata/project_data/gsp-sas/1km covariates/SSURGO/ssurgo_ptph100_rs_med.tif")
statsgo_ph100_r <- raster("D:/geodata/project_data/gsp-sas/1km covariates/SSURGO/statsgo_ph_r.i30-100.tif")
merge(natsgo_ph100_r, statsgo_ph100_r, overlap = TRUE, filename = "D:/geodata/project_data/gsp-sas/1km covariates/SSURGO/ssurgo_ptph100_rs_med_v2.tif", progress = "text", overwrite = TRUE)
  
```



# Reclassify MLRAs

```{r mlra}

# SSURGO CONUS
library(raster)
path <- "D:/geodata/project_data/gsp-sas/1km covariates/SSURGO"
lf <- list.files(path)
lf <- lf[grepl(".tif$", lf)]
ssurgo_rs <- stack(file.path(path, lf))
names(ssurgo_rs) <- gsub(".tif", "", lf)
ssurgo_vars <- names(ssurgo_rs)

# classify SSURGO derivatives
samp <- na.exclude(raster::sampleRegular(ssurgo_rs, 30000, sp = TRUE)@data)
samp.sc <- scale(samp[-1])
samp.dist <- dist(samp.sc)

samp.sil <- sapply(2:35, function(i) {
  cat("cluster", i, "\n")
  k <- kmeans(samp.sc, centers = i, iter.max = 100)
  k.sil <- summary(silhouette(k$cluster, samp.dist))$avg.width
})
plot(2:35, samp.sil, type = "l")
points(2:35, samp.sil)
samp.k <- kmeans(samp.sc, center = 19, iter.max = 100)
samp.k <- data.frame(mlra = as.character(samp$mlra_rs), cluster = as.character(samp.k$cluster))


# determine dominant cluster per MLRA
samp.k.p <- prop.table(table(samp.k$mlra, samp.k$cluster), margin = 1)
samp.k.p <- as.data.frame(samp.k.p, stringsAsFactors = FALSE)
names(samp.k.p) <- c("mlra", "cluster", "pct")
samp.k.p <- samp.k.p[order(samp.k.p$mlra, - samp.k.p$pct), ]
samp.k.p <- samp.k.p[!duplicated(samp.k.p$mlra), ]
names(samp.k.p)[1] <- "ID"
samp.k.p[1:2] <- lapply(samp.k.p[1:2], as.integer)


# reclassify MLRA layer
mlra_r <- ssurgo_rs$mlra_rs
mlra_r <- ratify(mlra_r)
rat <- levels(mlra_r)[[1]]
rat <- merge(rat, samp.k.p, by = "ID", all.x = TRUE, sort = FALSE)
levels(mlra_r) <- rat
mlra_r2 <- deratify(mlra_r, att = "cluster", filename = "D:/geodata/project_data/gsp-sas/1km covariates/SSURGO/mlra_rs_k19.tif")

```



# OCONUS

```{r}

library(raster)

fp <- "D:/geodata/project_data/gsp-sas/OCONUS/SSURGO_properties_no_ptf/Raster/"
lf <- list.files(path = fp, pattern = ".tif$", full.names = TRUE)

sapply(lf, function(x) {
  if (grepl("_pH", x)) {
    cat("applying pH ptf to ", x, "\n")
    r <- raster(x)
    x2 <- names(r)
    names(r) <- "ph_h2o"
    predict(r, ph_lm, filename = paste0("D:/geodata/project_data/gsp-sas/OCONUS/", x2, "_ptf.tif"), format = "GTiff", progress = "text", overwrite = TRUE)
    }
  })


predfun <- function(model, data) {
  v <- predict(model, data, type = "response")
}

sapply(lf, function(x) {
  if (grepl("_ec", x)) {
    cat("applying ec ptf2 to ", x, "\n")
    r <- raster(x)
    x2 <- names(r)
    names(r) <- "ec_predict_one_to_two"
    predict(r, ec_glm2, type = "response", filename = paste0("D:/geodata/project_data/gsp-sas/OCONUS/", x2, "_ptf.tif"), format = "GTiff", progress = "text", overwrite = TRUE, index = 1)
    }
  })

```