Data

The dataset of SOC stock by horizons for 96 pedon locations in Washington state contains 481 separate observations.

library(terra)
library(aqp)
library(ithir)
library(mpspline2)
library(dplyr)

setwd("/Users/Anthony/OneDrive - UW/University of Washington/Data and Modeling/")

wa_dat <- read.csv("SOIL CARBON/ANALYSIS/R/WA_SOC_controls/All_WA_horizons_spec_geoagelith.csv")
str(wa_dat)

Next we calculate top, center, and bottoms of the soil horizons

wa_dat_horC <- wa_dat |>
    dplyr::group_by(sample_ID) |>
    mutate(top = case_when(is.na(depth_cm - lag(depth_cm)) ~
        0, .default = lag(depth_cm)), bottom = depth_cm, center = abs(top -
        (top - bottom)/2)) |>
    dplyr::select(sample_ID, top, center, bottom, carbon_perc,
        carbon_stock_g_cm2)
str(wa_dat_horC)

Next, create a subset of the data with only the predictor variables


wa_hor_dat_sub <- wa_dat |>
    dplyr::select(-c(X, depth_cm, BD_g_cm3, carbon_perc, carbon_stock_g_cm2)) |>
    group_by(sample_ID) |>
    filter(row_number() == 1) |>
    dplyr::select(sample_ID, everything()) |>
    dplyr::arrange(site)
str(wa_hor_dat_sub)

In order to generate mapped predictions of SOC stock, we use adjust the dataset from observations of horizons of varying depths to observations of standard depth increments. Specifically, we use the mass-preserving spline method of Bishop et al (1999) (doi: 10.1016/S0016-7061(99)00003-8) to generate SOC stock estimates for intervals of 0-30cm, 30-60cm, and 60-100cm

wa_spl <- mpspline_tidy(obj = wa_dat_horC, var_name = "carbon_stock_g_cm2",
    d = c(0, 30, 60, 100))
wa_spl_dat_0_30_60_100 <- wa_spl$est_dcm |>
    left_join(y = wa_hor_dat_sub, by = join_by(sample_ID), relationship = "many-to-one") |>
    rename(SOC_stock_spline = SPLINED_VALUE, upper_depth = UD,
        lower_depth = LD)
str(wa_spl_dat_0_30_60_100)
write.csv(wa_spl_dat_0_30_60_100, file = "SOIL CARBON/ANALYSIS/R/WA_SOC_controls/wa_spl_dat_0_30_60_100.csv")