Add accumulateDepths()

[brownag]

A generic function for dealing with old-style O horizons and other changes in datum for soil surface.

I developed this routine originally at request of @jskovlin and refined it with ideas of @dylanbeaudette. It has been idling in the /misc/ folder. This is a PR to get some discussion going on the function. I added tests, docs, examples etc. but am happy to adjust this function in any way. @smroecker I know has worked on similar things. I would welcome any comments!

I anticipate it will be useful in these two related issues downstream of aqp, so would like to move forward with it

1. Extending OSD parsing tools SoilKnowledgeBase#25
2. https://github.com/ncss-tech/parse-osd/issues/12

---

Example

Linear Z-axis datum shift by profile ID

Running the function on a valid SPC allows you to do datum shifts up or down, either in bulk for whole SPC or by profile

library(aqp)
#> This is aqp 1.31
data(sp4)
depths(sp4) <- id ~ top + bottom
hz <- accumulateDepths(sp4,
    id = "id",
    hzdepths = c("top", "bottom"),
    hzname = "name",
    hzdatum = 5 * 1:length(sp4)
  )
plot(hz)

Treat reversed top depth/bottom depth as horizons above datum (e.g. old style O horizons)

Treat reversed top and bottom depths as depths above datum. This is a special interpretation of source data that could have become that way for a variety of data entry related reasons other than datum.

# example using old-style O horizons
hz <- read.table(text = "peiidref hzdept hzdepb hzname seqnum phiid
                    1        11      0      5      A      2   295
                    2        11      1      0     Oe      1   294
                    3        11      5     13     C1      3   296
                    4        11     13     58     C2      4   297
                    5        11     58    152     C3      5   298
                    6        13      0      5      A      2   303
                    7        13      1      0     Oe      1   302
                    8        13      5     25     Bw      3   304
                    9        13     25     61      C      4   305
                    10       13     61     NA      R      5   306
                    11      136      0     13     A1      3   695
                    12      136      1      0     Oe      2   694
                    13      136      2      1     Oi      1   693
                    14      136     13     61     C1      4   696
                    15      136     61     76     C2      5   697")

depths(hz) <- peiidref ~ hzdept + hzdepb
#> converting profile IDs from integer to character
#> Warning: Horizon bottom depths contain NA! Check depth logic with
#> aqp::checkHzDepthLogic()

hz_fixed <- accumulateDepths(hz,
                             id = "peiidref",
                             hzdepths = c("hzdept", "hzdepb"),
                             hzname = "hzname")
is_valid <- checkHzDepthLogic(hz_fixed)$valid

test0 <- subset(hz_fixed, !is_valid)
test1 <- subset(hz_fixed, is_valid)

origO <- subset(hz, grepl("O", hzname))
fixedO <- subset(hz_fixed, grepl("O", hzname))

par(mfrow=c(2,1), mar=c(0,0,3,2))
plotSPC(origO, max.depth = 25)
plotSPC(fixedO, max.depth = 25)

accumulateDepths()

Accumulate horizon depths, and reflect reversed depths, relative to new datum

Fix old-style organic horizon depths, or depths with a non-standard datum, by the "depth accumulation" method.

The "depth accumulation" method calculates thicknesses of individual horizons and then cumulative sums them after putting them in id + top depth order. The routine tries to determine context based on hzname and pattern. The main transformation is if a top depth is deeper than the bottom depth, the depths are reflected on the Z-axis (made negative). The data are then id + top depth sorted again, the thickness calculated and accumulated to replace the old depths.

This function uses several heuristics to adjust data before transformation and thickness calculation:

Regex matching of horizon designation patterns and similar

• matches of pattern where both top and bottom depth NA -> [0,1] [top,bottom] depth
• REMOVE horizons that do not match pattern where both top and bottom depths NA

Over-ride hzname handling with the sequence column argument seqnum

• if seqnum column specified "first record with NA hzname" is considered a pattern match if seqnum == 1

Trigger "fixing" with the fix argument:

• 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

Parameters:

• param x -- A data.frame or SoilProfileCollection
• param id -- unique profile ID. Default: NULL, if x is a SoilProfileCollection idname(x)
• param hzdepths -- character vector containing horizon top and bottom depth column names. Default: NULL, if x is a SoilProfileCollection horizonDepths(x)
• param hzname -- character vector containing horizon designation or other label column names. Default: NULL, if x is a SoilProfileCollection hzdesgnname(x)
• param hzdatum -- a numeric constant to add to accumulated depths. Default: 0
• param seqnum -- Optional: character vector containing record "sequence number" column name; used in-lieu of hzname (when NA) to identify "first" record in a profile
• param pattern -- pattern to search for in hzname to identify matching horizons to append the profile to
• param fix -- apply adjustments to missing (NA) depths and expand 0-thickness horizons? Default: TRUE

Returns: A horizon-level data.frame, suitable for promoting to SPC with depths<-, or a SoilProfileCollection, depending on the class of x.

[dylanbeaudette]

Thanks for pushing this forward, this is critical functionality that has been a long time coming. The "depth accumulation" method makes sense and is intuitive.

Overall, I like the approach but have some minor concerns about the function name and bundling of functionality. There appear to be three distinct things happening here: 1. vertical shifts, 2. flipping + shifting, 3. "fixing" problematic depths.

Vertical Shifting of Profiles

The basic interface for datum / z-axis shifts is intuitive and nicely vectorized. I'd suggest different terminology like "depth" for the manual pages, but that is minor. I'm on the fence about the use of datum, even though it is probably the best fit here. For this specific task, the function name is a little hard to guess / find.

Fixing / Flipping Old-Style O horizons

Excellent: the pattern-matching on horizon name with a seqnum-style backup is a good idea. I do think that we should encourage more strict-use of the hzdesgnname() metadata-setter in our examples and documentation. It is very easy to supply clobber the SPC if that isn't set or specified by argument correctly.

For the purposes of flipping old-style O horizons, the function name may be hard to find.

Fixing Horizon Depths (NA or equal)

I can see how this is related to the problems encountered with or caused by old-style O horizons. However, including this functionality reduces the modularity of the code. Also, 2/3 of this functionality exists in repairMissingHzDepths. Now that function is currently limited to the bottom-most horizons, so would miss anything with top == bottom elsewhere in the profile.

What do you think about the following?

• Vertical shifting is implemented as a stand-alone function, but still builds on the "accumulation of depths" algorithm, this becomes something like shiftZ, shiftVerticalDatum, or the like.
• Flipping old-style O horizons builds on the same algorithm but is called something different: this is a highly specialized task that should probably do that "one thing" as well / simply as possible. Perhaps flipHzAboveZero (not my favorite), flipHzAboveDatum, or similar.
• Consolidate all "fixing" of horizon depths (NA and equal) into a single function. repairMissingHzDepths is a starting place and is open to an overhaul and / or name change. It builds on the data.table optimizations you put in place earlier this year. It could be called fixHzDepths for clarity as it applies to both NA and equal depth "problems".

[brownag]

Thanks for reviewing this and all of the thoughtful suggestions.

Vertical Shifting of Profiles

The basic interface for datum / z-axis shifts is intuitive and nicely vectorized. I'd suggest different terminology like "depth" for the manual pages, but that is minor. I'm on the fence about the use of datum, even though it is probably the best fit here. For this specific task, the function name is a little hard to guess / find.

I don't think that "datum" is the wrong word to use here.

Fixing / Flipping Old-Style O horizons

Excellent: the pattern-matching on horizon name with a seqnum-style backup is a good idea.

@jskovlin can be thanked for that! I would not have come up with it on my own.

I do think that we should encourage more strict-use of the hzdesgnname() metadata-setter in our examples and documentation. It is very easy to supply clobber the SPC if that isn't set or specified by argument correctly.

This should now be addressed/addressable with the changes that came from #237

For the purposes of flipping old-style O horizons, the function name may be hard to find.

Fixing Horizon Depths (NA or equal)

I can see how this is related to the problems encountered with or caused by old-style O horizons. However, including this functionality reduces the modularity of the code.
Also, 2/3 of this functionality exists in repairMissingHzDepths. Now that function is currently limited to the bottom-most horizons, so would miss anything with top == bottom elsewhere in the profile.

This original routine that became this function preceded the repairMissingHzDepths() being available in aqp by two days. I think that we certainly could make repair missing hz depths more general, or have another method that fixes other types of logic errors (0-thickness, reversed depths).

• Vertical shifting is implemented as a stand-alone function, but still builds on the "accumulation of depths" algorithm, this becomes something like shiftZ, shiftVerticalDatum, or the like.

I think it would be easy to define one of these as a wrapper/alias of accumulateDepths() but am not sure how I would implement it as a stand alone function. In this routine the vertical shifting (or ability to shift vertically) comes as a side effect of the fact that all of the depths are re-accumulated for each profile. It just provides the option for each profile to have a starting value greater or less than zero.

• Flipping old-style O horizons builds on the same algorithm but is called something different: this is a highly specialized task that should probably do that "one thing" as well / simply as possible. Perhaps flipHzAboveZero (not my favorite), flipHzAboveDatum, or similar.

Again I think these could be an alias. I like "reflect" over "flip" but don't really care what the alias would be called.

• Consolidate all "fixing" of horizon depths (NA and equal) into a single function. repairMissingHzDepths is a starting place and is open to an overhaul and / or name change.

This would be fine. Is this something that must be implemented in order to merge the accumulateDepths() algorithm?

[dylanbeaudette]

Good to merge.