---
title: "Depth levels and the slabs they stand for"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Depth levels and the slabs they stand for}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = TRUE
)
```

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.

```{r setup}
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:

```{r thin-envelope}
thin <- as_envelope(fp, depth_min = 10, depth_max = 20)
terra::values(envelope_to_voxel(thin, depths = c(0, 100, 200)))
```

## 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:

```{r adjacent-upper}
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)))
```

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:

```{r adjacent-intersects}
terra::values(intersects_3d(upper, lower))

terra::values(intersects_3d(
  envelope_to_voxel(upper, depths = c(0, 100, 200)),
  envelope_to_voxel(lower, depths = c(0, 100, 200))
))
```

## 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:

```{r bounds}
shallow <- as_envelope(fp, depth_min = 60, depth_max = 70)

terra::values(envelope_to_voxel(shallow, depths = c(0, 100, 200)))
terra::values(envelope_to_voxel(shallow, depths = c(0, 100, 200),
                                bounds = "midpoint"))
```

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.