julienmoeys / soiltexture

The Soil Texture Wizard (R package soiltexture)
27 stars 10 forks source link

Bug in TT.text.transf.X #6

Open trevor-baker opened 3 years ago

trevor-baker commented 3 years ago

I am using package version 1.5.1, installed 2020-11-30.

The TT.text.transf.X function has a bug when the dat.ps.lim contains one or more sizes that are smaller than the first non-zero base.ps.lim size. The first column of the output dataframe contains only the value from the smallest input size class and the sum of the row does not come out to 100.

reprex:

psd.in <- data.frame(PS_0.001= c(5,10),
                         PS_0.002 = c(20,25), 
                         PS_0.02 = c(10,5), 
                         PS_0.05 = c(20,15), 
                         PS_2 = c(40,40),
                         PS_999 = c(5,5))
in.size <- c(0, 0.001, 0.002, 0.02, 0.05, 2, 999)
out.size <- c(0, 0.002, 0.05, 0.1, 0.5, 1, 2, 999)
psd.out <- as.data.frame( soiltexture::TT.text.transf.X(tri.data = psd.in, 
                                                            base.ps.lim = out.size, 
                                                            dat.ps.lim = in.size, 
                                                            text.tol=1) )

#need revised function below to run this part. This gives correct result.   
#psd.out2 <- as.data.frame( psd.transform(tri.data = psd.in, 
#                                             base.ps.lim = out.size, 
#                                             dat.ps.lim = in.size, 
#                                             text.tol=1) )

names(psd.out) <- c("PS_0.002", "PS_0.05", "PS_0.1", "PS_0.5", "PS_1", "PS_2", "PS_999")
rowSums(psd.in) #sum is 100
rowSums(psd.out) #sum is 80 and 75, i.e. 100 minus the PS_002 values. The value in PS_0.002 column is skipped over
rowSums(psd.out2) #sum is 100. PS_0.002 column is included in smallest size output

And here is the revised function. Just has a small piece to replace X[1] in the apply() function with the last column that has a PSD size smaller than or equal to the first non-zero base.ps.lim size. I hope it is clear enough.

psd.transform <- function (tri.data,
                           base.ps.lim,
                           dat.ps.lim,
                           text.sum = NULL,
                           text.tol = NULL,
                           tri.sum.tst = NULL,
                           tri.pos.tst = NULL){

  soiltexture::TT.auto.set(set.par = FALSE)
  soiltexture::TT.data.test.X(tri.data = tri.data, text.sum = text.sum,
                 text.tol = text.tol, tri.sum.tst = tri.sum.tst, tri.pos.tst = tri.pos.tst)

  tri.data <- t(apply(X = tri.data, MARGIN = 1, FUN = function(X) {
    cumsum(X)
  }))

  ps.end <- dim(tri.data)[2] + 1

  soiltexture::TT.check.ps.lim(base.ps.lim = base.ps.lim, dat.ps.lim = dat.ps.lim,
                  ps.lim.length = c(length(base.ps.lim), ps.end))
  if (base.ps.lim[1] != 0) {
    tri.data <- cbind(ZERO = rep(0, dim(tri.data)[1]), tri.data)
    ps.start <- 1
  } else {
    ps.start <- 2
  }

  base.ps.lim2 <- soiltexture::TT.dia2phi(base.ps.lim)
  dat.ps.lim2 <- soiltexture::TT.dia2phi(dat.ps.lim)

  old.col.nm <- colnames(tri.data)

  #this is the section I edited.
  tri.data2 <- t(apply(X = tri.data, MARGIN = 1,
                      FUN = function(X,base.ps.lim2, dat.ps.lim2) {

                              #need to find the columns in dat with PSD sizes smaller than the first non-zero base size.
                              dat.smaller.than.base <- which(sapply(dat.ps.lim2[ps.start:ps.end], function(x){
                                  x >= base.ps.lim2[ps.start]
                                } ))

                                c(sum(X[max(dat.smaller.than.base)]), #used to be X[1] - didn't capture exact matches and any sizes smaller
                                  diff(approx(x = dat.ps.lim2[ps.start:ps.end],
                                              y = X,
                                              xout = base.ps.lim2[ps.start:length(base.ps.lim)],
                                              method = "linear",
                                              rule = 1,
                                              ties = function(...) {
                                                        stop("error in TT.text.transf: Unexpected ties in text.cum = f(phi)")
                                                      })$y))
                            },
                base.ps.lim2, dat.ps.lim2))

  if (base.ps.lim[1] != 0) {
    tri.data2 <- tri.data2[, -1]
  }

  tri.data2 <- as.data.frame(tri.data2)

  colnames(tri.data2) <- paste(sep = "", "C", 1:dim(tri.data2)[2])

  return(tri.data2)

} #end function

Thanks for all your work on this package. It is very useful!