pkern_grid: Make a pkern grid list object

In deankoch/pkern: Fast Kriging with Product Kernels

Make a pkern grid list object

Description

Constructs a list representing a 2-dimensional regular spatial grid of data

Usage

pkern_grid(g, vals = TRUE)

Arguments

g	raster, matrix, numeric vector, or list (see details)
vals	logical indicating to include the data vector 'gval' in return list

Details

This function constructs pkern grid list objects, accepting 'RasterLayer' and 'RasterStack' inputs from the raster package, 'SpatRaster' objects from terra, as well as any non-complex matrix, or a list containing the vectorization of one. Empty grids can be initialized by specifying dimensions (and or setting vals=FALSE)

The function returns a list with the following 3-6 elements:

• gval: the data (if any) in column-major order with y descending, x ascending

• crs: character, the WKT string (if available) describing coordinate reference system

• idx_grid: integer vector, mapping rows of matrix gval to points on the grid

• gyx: a list containing the coordinates of the y and x grid lines in vectors y and x

• gres: the (numeric) y and x distances between grid lines

• gdim: the (integer) number of y and x grid lines

The last three elements are required to define a valid pkern grid list object. Note that regular grids must have equally spaced grid lines in gyx.

The input g can itself be a list containing some/all of these elements (including at least one of 'gdim' or 'gyx'), and the function will fill in missing entries wherever possible: If 'gval' is missing in the single-layer case, the function sets NAs; If 'res' 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 'res', if available) for each n in 'gdim'; and if 'gdim' is missing, it is set to equal the number of grid lines specified in (each vector of) 'gyx'.

gval can be a vector, for single-layer grids, or a matrix whose columns are a set of grid layers. In the matrix case, gval stores the observed data only, with NAs indicating by the mapping idx_grid. This mapping is assumed to be the same in all layers, but is only computed for the first layer. If a point is missing from one layer, it should be missing from all layers.

idx_grid is a length prod(gdim) integer vector with NAs for unobserved points, and otherwise the row number (in gval) of the observed point. These non-NA values must comprise seq(nrow(gval)) (ie all rows must be mapped), but they can be in any order. If gval is a matrix but idx_grid is missing, it is computed automatically (from the first column); and if idx_grid is supplied, but gval is a vector, it coerced to a 1-column matrix.

Scalar inputs to 'gdim', 'gres' are duplicated for both dimensions, and for convenience 'gdim' can be specified directly in g to initialize a simple grid; For example the call pkern_grid(list(gdim=c(5,5))) can be simplified to pkern_grid(list(gdim=5)) or pkern_grid(5).

Value

named list containing at least 'gyx', 'gres', and 'gdim' (see details)

Examples

# simple grid construction from dimensions
gdim = c(12, 10)
g = pkern_grid(g=list(gdim=gdim))
str(g)
str(pkern_grid(gdim, vals=FALSE))

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

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

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

# gres argument is ignored if a non-conforming gyx is supplied
gres_new = c(3, 4)
pkern_plot(pkern_grid(g=list(gyx=lapply(gdim, seq), gres=gres_new)))
pkern_plot(pkern_grid(g=list(gdim=gdim, gres=gres_new)))

# shorthand for square grids
all.equal(pkern_grid(2), pkern_grid(g=c(2,2)))

# example with random data
gdim = c(25, 25)
yx = as.list(expand.grid(lapply(gdim, seq)))
eg_vec = as.numeric( yx[[2]] %% yx[[1]] )
eg_mat = matrix(eg_vec, gdim)
g = pkern_grid(eg_mat)
pkern_plot(g, ij=T, zlab='j mod i')

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

# data vectors should be in R's default matrix vectorization order
all.equal(eg_vec, as.vector(eg_mat))
all.equal(g, pkern_grid(list(gdim=gdim, gval=as.vector(eg_mat))))

# multi-layer example from matrix
n_pt = prod(gdim)
n_layer = 3
mat_multi = matrix(rnorm(n_pt*n_layer), n_pt, n_layer)
g_multi = pkern_grid(list(gdim=gdim, gval=mat_multi))
str(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 = pkern_grid(list(gdim=gdim, gval=mat_multi))
str(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)

# vals=FALSE drops multi-layer information
pkern_grid(g=list(gdim=gdim, gval=mat_multi), vals=FALSE)

if( requireNamespace('raster') ) {

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

# convert to pkern list (notice only first layer was loaded by raster)
g = pkern_grid(r)
str(g)
pkern_plot(g)

# open a RasterStack - gval becomes a matrix with layers in columns
r_multi = raster::stack(r_path)
g_multi = pkern_grid(r_multi)
str(g_multi)
pkern_plot(g_multi, layer=1)
pkern_plot(g_multi, layer=2)
pkern_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 = pkern_grid(r_multi)
str(g_multi)
pkern_plot(g_multi, layer=1)
pkern_plot(g_multi, layer=2)
pkern_plot(g_multi, layer=3)

# open first layer only
g = pkern_grid(r[[1]])
str(g)
pkern_plot(g)

}
}

deankoch/pkern documentation built on Oct. 26, 2023, 8:54 p.m.