01_preprocess_data.R

[smroecker]

I noticed a few minor issues and have provided an alternative example. If you concur, I suggest forwarding an updated example to the other GBS training participants.

Issues

1. assumes the data has been filtered to only include horizons that overlap 0-25cm
2. excludes the Munsell dry value BS criteria (e.g. d_value <= 5)
3. the BS criteria example is complex

Suggestions

1. include a step demonstrating how to subset the data to only those horizons which overlap 0-25cm, if this is not done the BS criteria example fails
2. incorporate the d_value
3. generalize the BS criteria into a function to prevent operator error

Example comparing alternative approach and the original GBS training approach

install aqp GitHub version

remotes::install_github("ncss-tech/aqp", dependencies=FALSE, upgrade=FALSE, build=FALSE)

load lbiraries

library(aqp) library(dplyr)

load soil profile data

dat <- read.csv("https://raw.githubusercontent.com/FAO-GSP/GBSmap/master/input_data/soil_profiles.csv")

allocate Black Soil with aqp

dat2 <- allocate(dat, pedonid = "idp", hztop = "top", hzbot = "bot", OC = "oc", m_chroma = "w_chroma", m_value = "w_value", d_value = "d_value", CEC = "cec", BS = "bsat", to = "FAO Black Soil" )

dat2_prof <- mutate(dat2, z = ifelse(BS2 == TRUE, 2, 0), z = ifelse(BS1 == TRUE, 1, z) )

code from GBS training

x <- dat %>% mutate(bs_oc = if_else(condition = oc >= 1.2, true = 1, false = 0), bs_w_chroma = if_else(condition = w_chroma <= 3, true = 1, false = 0), bs_w_value = if_else(condition = w_value <= 3, true = 1, false = 0), bs_bot = if_else(condition = bot >= 25, true = 1, false = 0), bs_bsat = if_else(condition = bsat >= 50, true = 1, false = 0), bs_cec = if_else(condition = cec >= 25, true = 1, false = 0))

Classify the profiles as BS or no-BS

1st. We create two columns bs1 and bs2 to classify horizons

x <- x %>% mutate(bs1 = ifelse(bs_oc == 1 & bs_w_chroma == 1 & bs_w_value == 1 & bs_cec == 1 & bs_bsat == 1, yes = 1, no = 0), bs2 = ifelse(bs_oc == 1 & bs_w_chroma == 1 & bs_w_value == 1, yes = 1, no = 0))

then, a column with three values are created: 0 for non-BS, 1 for 1st cat. BS, and 2 for 2nd cat. BS

x2 <- x %>% group_by(idp) %>% summarise(x = first(x), y = first(y), n = n(), bs1 = sum(bs1)/n, bs2 = sum(bs2)/n, bottom = sum(bs_bot))

x2 <- x2 %>% transmute(idp = idp, x = x, y = y, z = ifelse(bs1 == 1 & bottom == 1, 1, ifelse(bs2 == 1 & bottom == 1, 2, 0)))

compare aqp & GBS training examples

x3 <- merge(x2, dat2_prof, by = "idp", all.x = TRUE) table(GBS = x3$z.x, aqp = x3$z.y)

idp <- x3[x3$z.x != x3$z.y, "idp"] View(x3[x3$idp %in% idp, ]) View(dat2[dat2$idp %in% idp, ])

[angelini75]

Code fixed https://github.com/FAO-GSP/GBSmap/blob/master/code/01_preprocess_data.R