Depth levels and the slabs they stand for

A SpatVoxel samples the water column at a set of depth levels, such as the World Ocean Atlas standard depths. A SpatEnvelope stores a continuous [depth_min, depth_max] interval per cell. Converting between the two with envelope_to_voxel(), or combining them with terra::mask(), intersect_3d(), intersects_3d() or volume(), needs a rule for which levels an interval occupies. This article sets out that rule.

The examples use a two-cell footprint, so every result can be printed in full.

library(ocean3d)

fp <- terra::rast(nrows = 1, ncols = 2, xmin = 0, xmax = 2, ymin = 0, ymax = 1)
terra::values(fp) <- c(1, 1)

A depth level stands for a slab, not a knife edge

A cell occupies a level when its envelope overlaps the slab that level stands for by more than a shared edge. This matters for thin intervals. An envelope of [10, 20] contains none of c(0, 100, 200), but it is clearly within the water the 0 m level represents:

thin <- as_envelope(fp, depth_min = 10, depth_max = 20)
terra::values(envelope_to_voxel(thin, depths = c(0, 100, 200)))
#>      presence_depth=0 presence_depth=100 presence_depth=200
#> [1,]                1                 NA                 NA
#> [2,]                1                 NA                 NA

Touching is not overlapping

An envelope of [0, 100] ends exactly where the 100 m slab begins, so it occupies the 0 m level alone:

upper <- as_envelope(fp, depth_min = 0, depth_max = 100)
lower <- as_envelope(fp, depth_min = 100, depth_max = 200)

terra::values(envelope_to_voxel(upper, depths = c(0, 100, 200)))
#>      presence_depth=0 presence_depth=100 presence_depth=200
#> [1,]                1                 NA                 NA
#> [2,]                1                 NA                 NA

Two envelopes that meet at 100 m therefore do not intersect. intersects_3d() gives the same answer for the envelopes themselves and for their voxels, which land on disjoint levels:

terra::values(intersects_3d(upper, lower))
#>      intersects
#> [1,]      FALSE
#> [2,]      FALSE

terra::values(intersects_3d(
  envelope_to_voxel(upper, depths = c(0, 100, 200)),
  envelope_to_voxel(lower, depths = c(0, 100, 200))
))
#>      intersects
#> [1,]      FALSE
#> [2,]      FALSE

Where the slab edges fall: bounds

Where the slab edges fall is a property of the dataset the levels came from, so bounds selects the convention:

  • "top" (the default) runs each slab from its level down to the next.
  • "midpoint" puts the edges halfway between neighbours. This is what World Ocean Atlas levels mean, so use it with WOA data.

An interval of [60, 70] lands in a different level under each:

shallow <- as_envelope(fp, depth_min = 60, depth_max = 70)

terra::values(envelope_to_voxel(shallow, depths = c(0, 100, 200)))
#>      presence_depth=0 presence_depth=100 presence_depth=200
#> [1,]                1                 NA                 NA
#> [2,]                1                 NA                 NA
terra::values(envelope_to_voxel(shallow, depths = c(0, 100, 200),
                                bounds = "midpoint"))
#>      presence_depth=0 presence_depth=100 presence_depth=200
#> [1,]               NA                  1                 NA
#> [2,]               NA                  1                 NA

The same bounds argument is accepted by terra::mask(), intersect_3d(), intersects_3d(), volume() and calc_volume_overlap() wherever they place an envelope on a voxel’s depth levels. See ?envelope_to_voxel for the full rule, including how the shallowest and deepest levels are clamped.