Calculate contour curves on a triangulation

Calculates contour curves and/or regions between them, for functions defined on a triangulation

Usage

tricontour(
  x,
  z,
  nlevels = 10,
  levels = pretty(range(z, na.rm = TRUE), nlevels),
  ...
)

# S3 method for class 'fm_mesh_2d'
tricontour(
  x,
  z,
  nlevels = 10,
  levels = pretty(range(z, na.rm = TRUE), nlevels),
  ...
)

# S3 method for class 'matrix'
tricontour(
  x,
  z,
  nlevels = 10,
  levels = pretty(range(z, na.rm = TRUE), nlevels),
  loc,
  ...
)

# S3 method for class 'list'
tricontour(
  x,
  z,
  nlevels = 10,
  levels = pretty(range(z, na.rm = TRUE), nlevels),
  loc,
  type = c("+", "-"),
  tol = 1e-07,
  ...
)

tricontourmap(
  x,
  z,
  nlevels = 10,
  levels = pretty(range(z, na.rm = TRUE), nlevels),
  ...
)

# S3 method for class 'fm_mesh_2d'
tricontourmap(
  x,
  z,
  nlevels = 10,
  levels = pretty(range(z, na.rm = TRUE), nlevels),
  ...
)

# S3 method for class 'matrix'
tricontourmap(
  x,
  z,
  nlevels = 10,
  levels = pretty(range(z, na.rm = TRUE), nlevels),
  loc,
  ...
)

# S3 method for class 'list'
tricontourmap(
  x,
  z,
  nlevels = 10,
  levels = pretty(range(z, na.rm = TRUE), nlevels),
  loc,
  type = c("+", "-"),
  tol = 1e-07,
  output = c("sp", "fm", "inla.mesh.segment"),
  ...
)

Arguments

x: An object generated by a call to fm_mesh_2d or fm_rcdt_2d, a triangle-vertex index matrix, or a list of triangulation information, list(loc, graph=list(tv)).
z: A vector containing the values to be contoured (NAs are allowed).
nlevels: Number of contour levels desired, if and only if levels is not supplied.
levels: Numeric vector of levels at which to calculate contour lines.
...: Additional arguments passed to the other methods.
loc: coordinate matrix, to be supplied when x is given as a triangle-vertex index matrix only.
type: "+" or "-", indicating positive or negative association. For +, the generated contours enclose regions where \(u_1 \leq z < u_2\), for - the regions fulfil \(u_1 < z \leq u_2\).
tol: tolerance for determining if the value at a vertex lies on a level.
output: The format of the generated output. Implemented options are "sp" (default) and "fm" (and deprecated "inla.mesh.segment", converted to "fm").

Value

For tricontour, a list with some of the fields that fm_segm objects have:

loc: A coordinate matrix
idx: Contour segment indices, as a 2-column matrix, each row indexing a single segment
grp: A vector of group labels. Each segment has a label, in 1,...,nlevels*2+1, where even labels indicate interior on-level contour segments, and odd labels indicate boundary segments between levels.

For tricontourmap, a list:

contour: A list of sp or fm_segm objects defining countour curves (level sets)
map: A list of sp or fm_segm objects enclosing regions between level sets

Author

Finn Lindgren finn.lindgren@gmail.com

Examples

if (FALSE) { # \dontrun{
if (require("fmesher") &&
  require("sp")) {
  ## Generate mesh and SPDE model
  n.lattice <- 20 # increase for more interesting, but slower, examples
  x <- seq(from = 0, to = 10, length.out = n.lattice)
  lattice <- fm_lattice_2d(x = x, y = x)
  mesh <- fm_rcdt_2d_inla(lattice = lattice, extend = FALSE, refine = FALSE)

  ## Generate an artificial sample
  sigma2.e <- 0.1
  n.obs <- 1000
  obs.loc <- cbind(
    runif(n.obs) * diff(range(x)) + min(x),
    runif(n.obs) * diff(range(x)) + min(x)
  )
  Q <- fm_matern_precision(mesh, alpha = 2, rho = 3, sigma = 1)
  x <- fm_sample(n = 1, Q = Q)
  A <- fm_basis(mesh, loc = obs.loc)
  Y <- as.vector(A %*% x + rnorm(n.obs) * sqrt(sigma2.e))

  ## Calculate posterior
  Q.post <- (Q + (t(A) %*% A) / sigma2.e)
  mu.post <- as.vector(solve(Q.post, (t(A) %*% Y) / sigma2.e))

  ## Calculate continuous contours
  tric <- tricontour(mesh,
    z = mu.post,
    levels = as.vector(quantile(x, c(0.25, 0.75)))
  )

  ## Discrete domain contours
  map <- contourmap(
    n.levels = 2, mu = mu.post, Q = Q.post,
    alpha = 0.1, compute = list(F = FALSE), max.threads = 1
  )

  ## Calculate continuous contour map
  setsc <- tricontourmap(mesh,
    z = mu.post,
    levels = as.vector(quantile(x, c(0.25, 0.75)))
  )

  ## Plot the results
  reo <- mesh$idx$lattice
  idx.setsc <- setdiff(names(setsc$map), "-1")
  cols2 <- contourmap.colors(map,
    col = heat.colors(100, 0.5, rev = TRUE),
    credible.col = grey(0.5, 0)
  )
  names(cols2) <- as.character(-1:2)

  par(mfrow = c(1, 2))
  image(matrix(mu.post[reo], n.lattice, n.lattice),
    main = "mean", axes = FALSE, asp = 1
  )
  plot(setsc$map[idx.setsc], col = cols2[idx.setsc])
  par(mfrow = c(1, 1))
}
} # }