horizon depths corruption by `accumulateDepths()` with NA

[dylanbeaudette]

I noticed some unusual output from accumulateDepths() when the source SPC contains NA in the top or bottom depths. Further inspection demonstrated unusual output when top == bottom, modification to depths in the absence of horizon names matching pattern, and when the bottom of an old-style O horizon is not directly adjacent to the top depth of the first mineral horizon.

It would be nice if this function could fail gracefully in the presence of "to little" context, such as when horizon depths are NA or top == bottom.

Test effects of NA:

library(aqp)
x <- quickSPC('01:AAAA|BBBBBBBB|CCCCCCC|CrCrCrCrCr')

# no change - expected
x
accumulateDepths(x)

# all top/bottom set to NA
y <- x
y$top[1] <- NA
accumulateDepths(y)

# unexpected changes to all depths
y <- x
y$bottom[1] <- NA
accumulateDepths(y)

Test effects of modifying either top or bottom depths

# unexpected changes to all depths
y <- x
y$bottom[1] <- 0
accumulateDepths(y)

# unexpected changes to all depths
y <- x
y$bottom[1] <- y$top[1]
accumulateDepths(y)

# unexpected changes to all depths
y <- x
y$bottom[2] <- y$top[2]
accumulateDepths(y)


# unexpected changes to all depths
y <- x
y$top[1] <- 10
y$bottom[1] <- 0
accumulateDepths(y)


# unexpected changes  to all depths
y <- x
y$top[1] <- 10
y$bottom[1] <- 0
y$name[1] <- 'O'
accumulateDepths(y)

Test specific problem, a single old-style O horizon with bottom depth matching top depth of next horizon.

# works as expected
y <- x
y$top[1] <- 10
y$bottom[1] <- 0
y$name[1] <- 'O'
y$top[2] <- 0
accumulateDepths(y)

[brownag]

I am not clear on what you are expecting to happen in these cases. Definitely happy to consider some additional heuristics for fixing things, but I think what we have in there now was an attempt at striking the balance of what is fixable vs. not following discussions on repairMissingHzDepths() fillHzGaps() and friends.

The premise of the function is based on summing up thicknesses. If you pass depths to the function that are incomplete or invalid then thicknesses cannot be calculated reliably.

Do you want an error message? Are you suggesting we include some more of these cases in default fix=TRUE?

Currently we "fix":

• Add 1 cm to bottom-most horizons with NA bottom depth
• Add 1 cm thickness to horizons with top and bottom depth equal
• Add 1 cm thickness to horizons with NA top depth and bottom depth 0

From what I can tell from your examples, these outcomes are expected because cumulative sum of horizon thicknesses propagates NA down the profile, similarly zeroing out a thickness of a horizon means everything underneath it becomes shallower. There is no attempt to preserve "gaps" in this approach (though, prior explicit filling of implicit gaps with fillHzGaps() would do the trick), the output is always a "continuous" profile (provided all depths are present or fixable).

This is another case of how much should we leave up to the user to QC their data prior to running analysis. I personally don't think all of our functions need to do fixing of gaps, repairing missing depths etc. internally. And I also don't think they should require perfect data.

[dylanbeaudette]

Good points; I remember the discussion and overall agreement that functionality in aqp should not attempt to "fix" all possible problems in noisy data. I also agree that some baseline QC would avoid many of the examples above. Based on tinkering this afternoon, I do have a couple of ideas for accumulateDepths():

• since the function relies on sorting via top depth, it should probably skip profiles with NA in any top depth
• skip profiles when fix=FALSE and no horizon name matches pattern
• additional documentation of what happens outside of the standard use-case of converting old-style O horizon depths to modern convention

Elaborating on a previous example, does the output of the following make sense despite the fact that the input data are invalid? In my mind, the function should either throw an error or skip such data.

# copy original
y <- x
# ruin a perfectly fine profile with a missing bottom depth in 1st horizon
y$bottom[1] <- NA
# do no perform any "fixing"
z <- accumulateDepths(y, fix = FALSE)

horizons(x)[, 1:3]
horizons(y)[, 1:3]
horizons(z)[, 1:3]

In a large collection, the following unusual conversion can go unnoticed.

  top bottom name
1   0      0    A
2   0     80    B
3  80    150    C
4 150    200   Cr

Here is another example, which I think demonstrates a bug with handling the pattern argument. Consider the original profile, with top two horizon depths arranged like an old-style O horizon.

y <- x
# old-style O horizon, but this is labeled 'A'
y$top[1] <- 10
y$bottom[1] <- 0
# make sure next horizon top depth starts at bottom of first
y$top[2] <- 0
z <- accumulateDepths(y, fix = FALSE, pattern = 'O')

horizons(x)[, 1:3]
horizons(y)[, 1:3]
horizons(z)[, 1:6]

The result doesn't make sense to me:

  top bottom name id hzID thk
1   0    120    B 01    2 120
2 120    110    A 01    1 -10
3 110    180    C 01    3  70
4 180    230   Cr 01    4  50

[brownag]

I agree with the concerns raised in your examples: The current approach of setting missing thickness to zero is risky--and should be optional (e.g., controlled by an na.rm argument) rather than automatic. Also, profiles with NA top depths are ususally problematic and should be skipped during processing.

Expecting accumulateDepths() to fix all data issues is unrealistic, but it was originally designed for "dirty" horizon data in a data.frame, not one that had gone through an SPC object initialization, so there was some fixing required; it would probably be helpful to have a better way to report problems and leave repairs to functions like repairMissingHzDepths(), but that all should probably be in a different function

There’s precedence for identifying a “surface” horizon, but skipping profiles with intractable issues (e.g., invalid depths, top == bottom, negative thickness) would be safer and more efficient. Any skipping logic should also consider datum shifts, not just the pattern matches.

The documentation could better highlight what happens outside standard use cases in the documentation of datum shift and old O--especially with messy or previously corrupted data.

Here is another example, which I think demonstrates a bug with handling the pattern argument. Consider the original profile, with top two horizon depths arranged like an old-style O horizon.

  top bottom name id hzID thk
1   0    120    B 01    2 120
2 120    110    A 01    1 -10
3 110    180    C 01    3  70
4 180    230   Cr 01    4  50

While this is very ugly, I think it is mostly an unfortunate byproduct of corrupting the SPC rather than a "bug"

When you set top depth of the A horizon to 10, you are breaking the internal order if you don't shift the position of the A horizon to where depths<- would have put it (same place accumulateDepths() puts it). This doesn't change the output in question, but it should change our expectation for the output. The proper position for modified A horizon is in the second slice based on the top depth. Since the A horizon doesn't match the pattern it does not get put at the top of the profile, and the negative depth reflection does not get applied to it. So, it gets put after the B horizon (which all along had a shallower top depth, as you replaced the original value). The A horizon, with a calculated thickness of -10, retains the same depth logic error it started with (top > bottom) underneath the B, and is moved to the correct slice index. Indeed, a bizarre situation, but it would work fine if either your pattern was e.g. "O|A" (i.e. you knew ahead of time you had some A horizons that used the O horizon style) or, better yet, your A horizon didn't have a depth logic error to begin with.

I can draft a PR so that the function:

• skips and flags profiles with NA top depths or invalid thicknesses
• avoid automatic replacement of NA thickness with zero
• provide clear messaging for data it cannot process safely
• leave major repairs to dedicated QC functions, possibly replacing the current fix functionality with calls to one or more existing or to-be-created functions

[dylanbeaudette]

Thanks. There is no rush--I can "fix" my immediate concerns by doing some baseline QC on the source SPC. I'm happy to help with testing and documentation along the way.