tricontour: Calculate contour curves on a triangulation

In finnlindgren/excursions: Excursion Sets and Contour Credibility Regions for Random Fields

Description

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 'inla.mesh'
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 'inla.mesh'
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", "inla.mesh.segment"),
  ...
)

Arguments

x	An object generated by a call to inla.mesh.2d or inla.mesh.create, 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 ≤q z < u_2, for - the regions fulfil u_1 < z ≤q 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 "inla.mesh.segment" (requires the INLA package).

Value

For tricontour, a list some of the same fields that inla.mesh.segment 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 inla.mesh.segment objects defining countour curves (level sets)
map	A list of sp or inla.mesh.segment objects enclosing regions between level sets

Author(s)

Finn Lindgren [e-mail]

Examples

## Not run: 
  if (require.nowarnings("INLA")) {
    ## 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 <- inla.mesh.lattice(x = x, y = x)
    mesh <- inla.mesh.create(lattice = lattice, extend = FALSE, refine = FALSE)
    spde <- inla.spde2.matern(mesh, alpha = 2)
    
    ## Generate an artificial sample
    sigma2.e <- 0.01
    n.obs <-1000
    obs.loc <- cbind(runif(n.obs) * diff(range(x)) + min(x),
                     runif(n.obs) * diff(range(x)) + min(x))
    Q <- inla.spde2.precision(spde, theta = c(log(sqrt(0.5)), log(sqrt(1))))
    x <- inla.qsample(Q = Q)
    A <- inla.spde.make.A(mesh = 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),
                               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)
    plot(setsc$map[idx.setsc], col = cols2[idx.setsc])
    par(mfrow = c(1,1))
  }

## End(Not run)

finnlindgren/excursions documentation built on Jan. 17, 2021, 12:20 a.m.