sk: Make a snapKrig grid list object

Description

Constructs snapKrig ("sk") class list, representing a 2-dimensional regular spatial grid

Usage

sk(..., vals = TRUE)

Value

a "sk" class list object

Arguments

...: raster, matrix, numeric vector, or list of named arguments (see details)
vals: logical indicating to include the data vector gval in return list

Details

This function accepts 'RasterLayer' and 'RasterStack' inputs from the raster package, 'SpatRaster' objects from terra, as well as any non-complex matrix, or a set of arguments defining the vectorization of one. It returns a sk class list containing at least the following three elements:

gdim: vector of two positive integers, the number of grid lines (n = their product)
gres: vector of two positive scalars, the resolution (in distance between grid lines)
gyx: list of two numeric vectors (lengths matching gdim), the grid line intercepts

and optionally,

crs: character, the WKT representation of the CRS for the grid (optional)
gval: numeric vector or matrix, the grid data
is_obs: logical vector indicating non-NA entries in the grid data
idx_grid: length-n numeric vector mapping rows of gval to grid points

Supply some/all of these elements (including at least one of gdim or gyx) as named arguments to .... The function will fill in missing entries wherever possible.

If gres is missing, it is computed from the first two grid lines in gyx; If gyx is missing, it is assigned the sequence 1:n (scaled by gres, if available) for each n in gdim; and if gdim is missing, it is set to the number of grid lines specified in gyx. gyx should be sorted (ascending order), regularly spaced (with spacing gres), and have lengths matching gdim.

Scalar inputs to gdim, gres are duplicated for both dimensions. For example the call sk(gdim=c(5,5)) can be simplified to sk(gdim=5), or sk(5).

For convenience, arguments can also be supplied together in a named list passed to .... If a single unnamed argument is supplied (and it is not a list) the function expects it to be either a numeric vector (gdim), a matrix, or a raster object.

Alternatively, you can supply an sk object as (unnamed) first argument, followed by individual named arguments. This replaces the named elements in the sk object then does a validity check.

Empty grids - with all data NA - can be initialized by setting vals=FALSE, in which case gval will be absent from the returned list). Otherwise gval is the column-vectorized grid data, either as a numeric vector (single layer case only) or as a matrix with grid data in columns. gval is always accompanied by is_obs, which supplies an index of NA entries (or rows)

A sparse representation is used when gval is a matrix, where only the non-NA entries (or rows) are stored. idx_grid in this case contains NA's were is_obs is FALSE, and otherwise contains the integer index of the corresponding row in gval. In the matrix case it is assumed that each layer (ie column) has the same NA structure. idx_grid is only computed for the first layer. If a point is missing from one layer, it should be missing from all layers.

Examples

Run this code


# simple grid construction from dimensions
gdim = c(12, 10)
g = sk(gdim)
summary(g)

# pass result to sk and get the same thing back
identical(g, sk(g))

# supply grid lines as named argument instead and get the same result
all.equal(g, sk(gyx=lapply(gdim, function(x) seq(x)-1L)))

# display coordinates and grid line indices
plot(g)
plot(g, ij=TRUE)

# same dimensions, different resolution, affecting aspect ratio in plot
gres_new = c(3, 4)
plot(sk(gdim=gdim, gres=gres_new))

# single argument (unnamed) can be grid dimensions, with shorthand for square grids
all.equal(sk(gdim=c(2,2)), sk(c(2,2)))
all.equal(sk(2), sk(gdim=c(2,2)))

# example with matrix data
gdim = c(25, 25)
gyx = as.list(expand.grid(lapply(gdim, seq)))
eg_vec = as.numeric( gyx[[2]] %% gyx[[1]] )
eg_mat = matrix(eg_vec, gdim)
g = sk(eg_mat)
plot(g, ij=TRUE, zlab='j mod i')

# y is in descending order
plot(g, xlab='x = j', ylab='y = 26 - i', zlab='j mod i')

# this is R's default matrix vectorization order
all.equal(eg_vec, as.vector(eg_mat))
all.equal(g, sk(gdim=gdim, gval=as.vector(eg_mat)))

# multi-layer example from matrix
n_pt = prod(gdim)
n_layer = 3
mat_multi = matrix(stats::rnorm(n_pt*n_layer), n_pt, n_layer)
g_multi = sk(gdim=gdim, gval=mat_multi)
summary(g_multi)

# repeat with missing data (note all columns must have consistent NA structure)
mat_multi[sample.int(n_pt, 0.5*n_pt),] = NA
g_multi_miss = sk(gdim=gdim, gval=mat_multi)
summary(g_multi_miss)

# only observed data points are stored, idx_grid maps them to the full grid vector
max(abs( g_multi[['gval']] - g_multi_miss[['gval']][g_multi_miss[['idx_grid']],] ), na.rm=TRUE)

# single bracket indexing magic does the mapping automatically
max(abs( g_multi[] - g_multi_miss[] ), na.rm=TRUE)

# vals=FALSE drops multi-layer information
sk(gdim=gdim, gval=mat_multi, vals=FALSE)

# raster import examples skipped to keep CMD check time < 5s on slower machines
# \donttest{
if( requireNamespace('raster') ) {

# open example file as RasterLayer
r_path = system.file('external/rlogo.grd', package='raster')
r = raster::raster(r_path)

# convert to sk (notice only first layer was loaded by raster)
g = sk(r)
summary(g)
plot(g)

# open a RasterStack - gval becomes a matrix with layers in columns
r_multi = raster::stack(r_path)
g_multi = sk(r_multi)
summary(g_multi)
plot(g_multi, layer=1)
plot(g_multi, layer=2)
plot(g_multi, layer=3)

# repeat with terra
if( requireNamespace('terra') ) {

# open example file as SpatRaster (all layers loaded by default)
r_multi = terra::rast(r_path)
g_multi = sk(r_multi)
summary(g_multi)

# open first layer only
g = sk(r[[1]])
summary(g)

}
}
# }