Functional functions :)

transforme_om_oc, Model linear relation between organic matter and organic carbon and estimate organic carbon values from organic matter data
estimate_h, checks for space between samples and, if any, divide this space between the previous and next sample to return sample thickness withouth gaps in the core

estimate_stocks, Estimates the carbon stock from soil core data until a specific depth, 100 cm by default. If the core do not reach the standardization depth it extrapolate the stock from a linear model between the organic carbon accumulated mass and depth.

test_extrapolations, subset those cores that reach the standardization depth and estimates the stock (observed stock), estimate the stock from the linear relation of organic carbon accumulated mass and depth using the 90, 75, 50 and 25% top length of the indicated standardization depth. Compares the observed stock with the estimated stocks by extrapolation.

estimate_fluxes, estimate the average organic carbon flux to the soil in a indicated time frame (by default last 100 years) from the organic carbon concentration and ages obtained from a age-depth or age-accumulated mass model

R/estimate_flux.R

@@ -0,0 +1,83 @@
+#' Estimate the average organic carbon flux
+#'
+#' @description estimate the average organic carbon flux to the soil in a indicated time frame (by default last 100 years) from the organic carbon concentration and ages obtained from a age-depth or age-accumulated mass model
+#'
+#' @param df A [data.frame] with, at least, columns: Core.ID, DMin (minimum depth of the sample), DMax (maximum depth of the sample), DBD (dry bulk density), POC (organic carbon %), Age (age of the sample obtained from a age-depth or age-accumulated mass model)
+#' @param TimeFrame standardization time frame, by default 100 years
+#'
+#' @return [data.frame] with columns Core.id, F.WC (organic carbon fluxes at the whole core), A.Max (maximum age of the core), and Flux (average organic carbon flux at the indicated time frame)
+#' @export
+#'
+#' @examples
+
+
+estimate_flux<- function(df=NULL,TimeFrame=100) {
+
+  # class of the dataframe
+  if (is.data.frame(df)==FALSE) {stop("The data provided is not class data.frame, please check data and transforme")}
+  if (is.numeric(TimeFrame)==FALSE) {stop("The TimeFrame provided is not class numeric, please check data and transforme")}
+
+  # name of the columns
+  if ("CoreID" %in% colnames(df)==FALSE) {stop("There is not column named CoreID. Please, check necessary columns in functions documentation")}
+  if ("DMin" %in% colnames(df)==FALSE) {stop("There is not column named Min.D. Please, check necessary columns in functions documentation")}
+  if ("DMax" %in% colnames(df)==FALSE) {stop("There is not column named Max.D. Please, check necessary columns in functions documentation")}
+  if ("DBD" %in% colnames(df)==FALSE) {stop("There is not column named DBD. Please, check necessary columns in functions documentation")}
+  if ("fOC" %in% colnames(df)==FALSE) {stop("There is not column named fOC. Please, check necessary columns in functions documentation")}
+  if ("Age" %in% colnames(df)==FALSE) {stop("There is not column named Age. Please, check necessary columns in functions documentation")}
+
+  # class of the columns
+  if (is.numeric(df$DMin)==FALSE) {stop("Minimum depth data is not class numeric, please check")}
+  if (is.numeric(df$DMax)==FALSE) {stop("Maximum depth data is not class numeric, please check")}
+  if (is.numeric(df$DBD)==FALSE) {stop("Dry Bulk Density data is not class numeric, please check")}
+  if (is.numeric(df$fOC)==FALSE) {stop("Organic carbon data is not class numeric, please check")}
+  if (is.numeric(df$Age)==FALSE) {stop("Age data is not class numeric, please check")}
+
+
+
+  #select those cores with chronological models
+  df<-df[!is.na(df$Age),]
+  df<-df[!is.na(df$fOC),]
+
+  df<-estimate_h (df)
+
+  X<-split(df, df$CoreID)
+
+  BCF <- data.frame(CoreID=character(),
+                    F.WC=numeric(),
+                    A.Max=numeric(),
+                    Flux=numeric())
+
+
+  for(i in 1:length(X)) {
+    BCF[i,1]<-names(X[i])
+    Data<-as.data.frame(X[i])
+    colnames(Data)<-colnames(df)
+
+    if(nrow(Data)<3) next
+
+    else{
+
+      Data <-Data |> dplyr::mutate (OCgcm2 = DBD*(fOC/100)*h)
+
+      #estimation of the average carbon flux for the whole core (OC stock/Max Age)
+      BCF[i,2]<-(sum(Data[,which(colnames(Data)=="OCgcm2")]))/max(Data$Age)
+      BCF[i,3]<-max(Data$Age)
+
+      #estimation of the average carbon flux for the selected TimeFrame (OC stock last 100 yrs/TimeFrame)
+
+      if (max(Data$Age)==TimeFrame) {
+
+        BCF[i,4]<-((sum(Data[c(1:(nrow(Data))),which(colnames(Data)=="OCgcm2")]))/TimeFrame)
+
+      } else {
+        Data<-Data[c(1:(length(which(Data$Age <=TimeFrame))+1)),]
+
+        BCF[i,4]<-(
+          (sum(Data[c(1:(nrow(Data)-1)),which(colnames(Data)=="OCgcm2")])+
+                       (Data[nrow(Data),which(colnames(Data)=="OCgcm2")]/((max(Data$Age)-Data[(nrow(Data)-1),which(colnames(Data)=="Age")])))
+                        *(TimeFrame-Data[(nrow(Data)-1),which(colnames(Data)=="Age")]))/TimeFrame)
+      }
+
+    }}
+  return(BCF)
+}

R/estimate_h.R

@@ -0,0 +1,90 @@
+#' Estimation of the thickness of the sample
+#'
+#' @description checks for space between samples and, if any, divide this space between the previous and next sample to return sample thickness withouth gaps in the core
+#'
+#' @param df A [data.frame] with with, at least, columns CoreID, DMin (minimum depth of the sample)and DMax (maximum depth of the sample)
+#'
+#' @return the initial [data.frame] with three additional columns: EMin (estimated minimum depth of the sample), EMax (estimated maximum depth of the sample) and h (estimated thikness of the sample)
+#' @export
+#'
+#' @examples
+estimate_h <- function(df = NULL) {
+
+  # class of the dataframe
+  if (is.data.frame(df)==FALSE) {stop("The data provided is not class data.frame, please chaeck data and transforme")}
+
+  # name of the columns
+  if ("CoreID" %in% colnames(df)==FALSE) {stop("There is not column named CoreID. Please, check necessary columns in functions documentation")}
+  if ("DMin" %in% colnames(df)==FALSE) {stop("There is not column named DMin. Please, check necessary columns in functions documentation")}
+  if ("DMax" %in% colnames(df)==FALSE) {stop("There is not column named DMax. Please, check necessary columns in functions documentation")}
+
+  # class of the columns
+  if (is.numeric(df$DMin)==FALSE) {stop("Minimum depth data is not class numeric, please check")}
+  if (is.numeric(df$DMax)==FALSE) {stop("Maximum depth data is not class numeric, please check")}
+
+  #check for NAs in depth columns
+  if (sum(is.na(df$DMin))>0){stop("Samples minimun depth column has NAs, please check")}
+  if (sum(is.na(df$DMax))>0){stop("Samples maximun depth column has NAs, please check")}
+
+
+  # create individual data frames per each core
+
+  df$CoreID <- factor(df$CoreID, levels=unique(df$CoreID))
+  X<-split(df, df$CoreID)
+
+
+  columns<-c("EMin","EMax","h")
+  Fdf2 = data.frame(matrix(nrow = 0, ncol = length(columns)))
+  colnames(Fdf2) = columns
+
+
+  for(i in 1:length(X)) {
+
+    Data<-as.data.frame(X[i])
+    colnames(Data)<-colnames(df)
+
+    #check if there is spaces between samples (e.g, first sample ends at 5 cm and next starts at 7)
+    space<- c()
+
+    for (j in 1:(nrow(Data)-1)) {
+
+      # if there are no spaces between samples min and maximun depth of samples remain the same
+      if (Data[j,which(colnames(Data)=="DMax")] == Data[j+1,which(colnames(Data)=="DMin")]) {
+        space[j]<-FALSE} else {space[j]<-TRUE}}
+
+    if (any(space==TRUE)) {
+      # if there are spaces between samples it estimate the medium point between the maximum depth of the sample and the minimun depth of the next sample
+      # and divide that distance between both samples
+      Data <- cbind(Data, EMin=NA, EMax=NA)
+      Data[1,"EMin"]<-0
+      Data[nrow(Data),"EMax"]<-Data[nrow(Data),"DMax"]
+      for (j in 1:(nrow(Data)-1)) {
+        if(space[j]==TRUE) {
+          Data[j,"EMax"]<-Data[j,"DMax"]+((Data[j+1,"DMin"]-Data[j,"DMax"])/2)
+          Data[j+1,"EMin"]<-Data[j,"DMax"]+((Data[j+1,"DMin"]-Data[j,"DMax"])/2)} else {
+            Data[j,"EMax"]<-Data[j,"DMax"]
+            Data[j+1,"EMin"]<-Data[j+1,"DMin"]}}
+
+    }  else{
+      Data <- cbind(Data, EMin=NA, EMax=NA)
+      Data$EMin<-Data$DMin
+      Data$EMax<-Data$DMax
+
+    }
+
+    Data <- cbind(Data, h=NA)
+
+    #estimation of the thickness of the sample (h) from the new minimun and max depth of the sample
+
+    Data<- Data |> dplyr::mutate (h = EMax-EMin)
+
+    temp<-cbind(Data$EMin, Data$EMax, Data$h)
+    colnames(temp)<-colnames(Fdf2)
+    Fdf2<-rbind(Fdf2, temp)
+
+  }
+  Fdf<-cbind(df, Fdf2)
+
+  return(Fdf)
+}
+

R/estimate_stock.R

@@ -0,0 +1,94 @@
+#' Organic Carbon Stock estimation
+#'
+#' @description Estimates the carbon stock from soil core data until a specific depth, 100 cm by default. If the core do not reach the
+#' standardization depth it extrapolate the stock from a linear model between the organic carbon accumulated mass and depth.
+#'
+#' @param df A [data.frame] with, at least, columns CoreID, DMin (minimum depth of the sample), DMax (maximum depth of the sample), DBD (dry bulk density), fOC (organic carbon %)
+#' @param Depth standardization soil depth, by default 100 cm.
+#'
+#' @return [data.frame] with columns CoreID, S.WC (organic carbon stock at the whole core), D.Max (maximum depth of the core), and Stock (organic carbon stock at the standardized depth)
+#' @export
+#'
+#' @examples
+
+estimate_stock <- function(df = NULL, Depth = 100) {
+
+  # class of the dataframe
+  if (is.data.frame(df)==FALSE) {stop("The data provided is not class data.frame, please chaeck data and transforme")}
+  if (is.numeric(Depth)==FALSE) {stop("The Depth provided is not class numeric, please chaeck data and transforme")}
+
+  # name of the columns
+  if ("CoreID" %in% colnames(df)==FALSE) {stop("There is not column named CoreID. Please, check necessary columns in functions documentation")}
+  if ("DMin" %in% colnames(df)==FALSE) {stop("There is not column named Min.D. Please, check necessary columns in functions documentation")}
+  if ("DMax" %in% colnames(df)==FALSE) {stop("There is not column named Max.D. Please, check necessary columns in functions documentation")}
+  if ("DBD" %in% colnames(df)==FALSE) {stop("There is not column named DBD. Please, check necessary columns in functions documentation")}
+  if ("fOC" %in% colnames(df)==FALSE) {stop("There is not column named fOC. Please, check necessary columns in functions documentation")}
+
+  # class of the columns
+  if (is.numeric(df$DMin)==FALSE) {stop("Minimum depth data is not class numeric, please check")}
+  if (is.numeric(df$DMax)==FALSE) {stop("Maximum depth data is not class numeric, please check")}
+  if (is.numeric(df$DBD)==FALSE) {stop("Dry Bulk Density data is not class numeric, please check")}
+  if (is.numeric(df$fOC)==FALSE) {stop("Organic carbon data is not class numeric, please check")}
+
+
+
+  df<-df[!is.na(df$fOC),]
+
+  # estimate thickness of the sample
+
+  dfh<-estimate_h(df)
+
+  # estimate stocks
+
+  X<-split(dfh, dfh$CoreID)
+
+  BCS <- data.frame(CoreID=character(),
+                    S.WC=numeric(),
+                    D.Max=numeric(),
+                    Stock=numeric())
+
+  for(i in 1:length(X)) {
+    BCS[i,1]<-names(X[i])
+    Data<-as.data.frame(X[i])
+    colnames(Data)<-colnames(dfh)
+
+    if(nrow(Data)<3) next
+
+    else{
+
+      #estimation of carbon g cm2 per sample, OCgcm2= carbon density (g cm3) by thickness (h)
+      Data <-Data |> dplyr::mutate (OCgcm2 = DBD*(fOC/100)*h)
+
+      #estimation of the OC stock in the whole core
+      BCS[i,2]<-sum(Data[,which(colnames(Data)=="OCgcm2")])
+      BCS[i,3]<-max(Data[,which(colnames(Data)=="EMax")])
+
+      #if core exactly the standarization depth, we keep the stock of the whole core
+      if(max(Data$EMax)==Depth) {BCS[i,4]<-sum(Data[,which(colnames(Data)=="OCgcm2")])}
+
+
+      else{
+
+        # if the core longer than the standardization depth we estimate the stock until that depth
+        if (max(Data$EMax)>=Depth)
+
+        {
+          Data<-Data[c(1:(length(which(Data$EMax <=Depth))+1)),]
+
+          if(nrow(Data)<3) next
+
+          else{
+
+          BCS[i,4]<-(((sum(Data[c(1:(nrow(Data)-1)),which(colnames(Data)=="OCgcm2")]))+
+                        (Data[nrow(Data),which(colnames(Data)=="OCgcm2")]/((max(Data$EMax)-Data[(nrow(Data)-1),which(colnames(Data)=="EMax")]))
+                         *(Depth-Data[(nrow(Data)-1),which(colnames(Data)=="EMax")]))))}}
+
+        #if core shorter than than the standardization depth we model the OC acumulated mass with depth and predict the stock at that depth
+        else {
+
+          Data <-Data |> dplyr::mutate (OCM = cumsum(OCgcm2))
+          model<-lm(OCM ~ EMax, data=Data)
+          BCS[i,4]<-coef(model)[1] + Depth*coef(model)[2]}}
+    }}
+  return(BCS)
+}