obrl-soil / mpspline2 Goto Github PK

So GSIF::mpspline() and ithir::ea_spline() both iterate over 1:n when making 1cm spline predictions, so that 1 = 0-1cm. mpspline2 and splinetool.exe iterate over 0:n-1, so that 0 = 0-1cm. The index, as k, is fed into the 1cm fit equations - see lines 246-285 of mpspline and lines 193-207 of ea_spline.

I changed my indexing to match splinetool.exe's outputs on purpose, but now I'm having second thoughts. Not sure which approach is actually correct. With 0-indexing, the 0-1cm prediction

alfa[h] + b0[h] * (k - s[[2]][h]) + gamma[h] * (k - s[[2]][h])^2

simplifies out to alfa(h1) + gamma(h1) as k and s[[2]]h are both 0. Not sure if that's a problem.

Turns out the Queensland Maths B curriculum ca. 2001 leaves me somewhat ill-equipped to take a critical look at the TMSE estimator 🙃, so I'm asking for some context/explanation here:

The TMSE code in existing implementations depends on a variance figure s2 calculated from the full input dataset (all supplied sites). In GSIF::mpspline() it appears to be 5% of the dataset standard deviation, squared. In ithir::ea_spline() it appears to be 5% of the mean, squared.

a) assuming use of mean is a bug, yes?
b) why take 5%?
c) doesn't this make the TMSE heavily dependant on input dataset size? Is that a problem?

There are also a couple of lines in the TMSE code that don't go anywhere, and I'm not sure what to do with them. @brendo1001, are you able to help out please?

Howdy. This is minor, but worth considering since the internal structure of SPC objects is subject to change. This is more important as aqp 2.0 looms (asymptotically?) closer.

The following code works fine but can / will break when internals change (sorry).

mpspline_conv.SoilProfileCollection <- function(obj = NULL) {
  dc <- obj@depthcols
  ic <- obj@idcol
  ac <- names(obj@horizons)[-which(names(obj@horizons) %in% c(dc, ic))]
  data.frame(c(obj@horizons[ic],
               obj@horizons[dc],
               obj@horizons[ac]), stringsAsFactors = FALSE)
}

Safely get everything you need like this (please). Comments / whitespace added for clarity.

# profile ID field name, set at SPC init time
ic <- idname(obj)

# horizon depth field names, set at SPC init time
dc <- horizonDepths(obj)

# horizon level attributes, includes IDs and depths
h <- horizons(obj)

# horizon level field names, minus IDs and depths
ac <- names(h)[-which(names(h) %in% c(ic, dc))]

# re-order to expected format: ID, top, bottom, {everything else}
h <- h[, c(ic, dc, ac)]

# done, stringsAsFactors logic mirrors source data
return(h)

I think that this accomplishes what is intended. A couple of things gained:

• safe getters used vs slot names
• single extraction of horizon level attributes: this may become more expensive in the future
• horizon level attributes are re-ordered without altering strings/factors set prior to function call

Let me know if I missed something.

Interested in integrating the wrapper functions from here?

• mpsplineSPC()
• singleMPS2SPC()

https://gist.github.com/dylanbeaudette/14aa4b4d43093f90f85d2379c2ae48cc

I am not partial to the names and would be happy to write docs / tests. I think that there might be a lot of interest in SPC --- MPS ---> SPC type work-flow.

If not, I'll likely add to aqp and add mpspline2 to the SUGGESTS field. Either way is OK with me.

Embedding d into the new style output object would be very helpful.