Title: | Vector Fields from Spatial Time Series of Population Abundance |
---|---|
Description: | Functions for converting time series of spatial abundance or density data in raster format to vector fields of population movement using the digital image correlation technique. More specifically, the functions in the package compute cross-covariance using discrete fast Fourier transforms for computational efficiency. Vectors in vector fields point in the direction of highest two dimensional cross-covariance. The package has a novel implementation of the digital image correlation algorithm that is designed to detect persistent directional movement when image time series extend beyond a sequence of two raster images. |
Authors: | Devin Goodsman [aut, cre] |
Maintainer: | Devin Goodsman <[email protected]> |
License: | GPL (>= 3) |
Version: | 0.1.2 |
Built: | 2024-10-27 04:35:04 UTC |
Source: | https://github.com/goodsman/icvectorfields |
Calculates a displacement field based on the cross-covariance of two input rasters presumably representing spatial population abundance or density at two different instances of time.
DispField(inputrast1, inputrast2, factv1, facth1, restricted = FALSE)
DispField(inputrast1, inputrast2, factv1, facth1, restricted = FALSE)
inputrast1 |
a raster as produced by terra::rast |
inputrast2 |
a raster of equivalent dimension to inputrast1 as produced by terra::rast |
factv1 |
an odd integer for the vertical dimension of sub-grids |
facth1 |
an odd integer for the horizontal dimension of sub-grids |
restricted |
logical (TRUE or FALSE) |
The input rasters are first converted to equivalent matrices. The function then divides the domain up into sub-grids of size factv1 X facth1, which are vertical and horizontal sub-grid dimensions.
If restricted is set to FALSE (the default), the function computes cross-covariance between each sub-grid of the first input raster and the entirety of the second input raster and then uses the location of maximum cross-covariance to estimate displacement in the vertical and horizontal directions from the centre of each sub-grid.
If restricted is set to TRUE, the function uses cross-covariance between each sub-grid in the first input raster and the equivalent sub-grid in the second input raster to estimate vertical and horizontal displacement.
Reference coordinates and cell size are extracted from the first input raster such that the locations from whence displacement is estimated as well as displacement estimates can be expressed in the units of the projected coordinates.
The coordinates are assumed to increase vertically and horizontally from the lower left corner of the two-dimensional domain.
Caution is warranted when defining the sub-grid dimensions because the function can produce erroneous results when sub-grids are too small.
A data frame is returned with the following column names: rowcent, colcent, frowmin, frowmax, fcolmin, fcolmax, centx, centy, dispx, and dispy. The rowcent and colcent column names are the row and column indices for the center of each sub-grid; frowmin and frowmax are the sub-grid minimum and maximum row indices; fcolmin and fcolmax are the sub-grid minimum and maximum column indices; centx and centy are the projected coordinates of the centre of the subgrid derived from the raster input files; dispx and dispy are the displacement in the horizontal and vertical directions in the same units as the projected coordinates of the raster input files.
DispFieldbb
for a similar function using a bounding
box to define a focal region, DispFieldST
for a version
designed to quantify persistent directional movement when the time series
features more than two time instances, DispFieldSTall
for a
version designed to quantify persistent directional movement when velocity
is variable in space, and Xcov2D
for demonstration of how
two-dimensional cross-covariance is used to determine displacement (see
examples of Xcov2D function documentation).
(Mat1 <- matrix(rep(c(1:5, 0, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat2 <- matrix(rep(c(0, 1:5, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) (VFdf1 <- DispField(rast1, rast2, factv1 = 9, facth1 = 9)) # The second raster is shifted right by 1 unit relative to the first raster # dispx = 1
(Mat1 <- matrix(rep(c(1:5, 0, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat2 <- matrix(rep(c(0, 1:5, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) (VFdf1 <- DispField(rast1, rast2, factv1 = 9, facth1 = 9)) # The second raster is shifted right by 1 unit relative to the first raster # dispx = 1
Calculates a displacement field based on the cross-covariance of two input
rasters presumably representing spatial population abundance or density at
two different instances of time. This version differs from
DispField
in that the user defines a bounding box that
determines a single sub-grid. The center of the bounding box is the location
from whence displacement is estimated.
DispFieldbb( inputrast1, inputrast2, rowmn, rowmx, colmn, colmx, restricted = FALSE )
DispFieldbb( inputrast1, inputrast2, rowmn, rowmx, colmn, colmx, restricted = FALSE )
inputrast1 |
a raster as produced by terra::rast |
inputrast2 |
a raster of equivalent dimension to inputrast1 as produced by terra::rast |
rowmn |
an integer denoting the minimum row index of the sub-grid |
rowmx |
an integer denoting the maximum row index of the sub-grid |
colmn |
an integer denoting the minimum column index of the sub-grid |
colmx |
an integer denoting the maximum column index of the sub-grid |
restricted |
logical (TRUE or FALSE) |
The input rasters are first converted to equivalent matrices. If restricted is set to FALSE (the default), the function computes cross-covariance between the sub-grid of the first input raster and the entirety of the second input raster and then uses the location of maximum cross-covariance to estimate displacement in the vertical and horizontal directions from the centre of the sub-grid.
If restricted is set to TRUE, the function uses cross-covariance between the sub-grid of the first input raster and the equivalent sub-grid of the second input raster to estimate vertical and horizontal displacement.
Reference coordinates and cell size are extracted from the first input raster such that the locations from whence displacement is estimated as well as displacement estimates can be expressed in the units of the projected coordinates.
The coordinates are assumed to increase vertically and horizontally from the lower left corner of the two-dimensional domain.
Caution is warranted when defining the bounding box because the function can produce erroneous results when the bounding box is too small.
A data frame is returned with the following column names: rowcent, colcent, frowmin, frowmax, fcolmin, fcolmax, centx, centy, dispx, and dispy. The rowcent and colcent column names are the row and column indices for the center of the sub-grid; frowmin and frowmax are the sub-grid minimum and maximum row indices; fcolmin and fcolmax are the sub-grid minimum and maximum column indices; centx and centy are the projected coordinates of the centre of the subgrid derived from the raster input files; dispx and dispy are the displacement in the horizontal and vertical directions in the same units as the projected coordinates of the raster input files.
DispField
for a similar function with a grid of focal
regions, DispFieldSTbb
for a version designed to quantify
persistent directional movement when the time series features more than two
time instances, DispFieldSTbball
for a version designed to
quantify persistent directional movement when velocity is variable in
space, and Xcov2D
for demonstration of how two-dimensional
cross-covariance is used to determine displacement (see examples of Xcov2D
function documentation).
rseq <- stats::runif(72) Mat1 <- matrix(rep(0, 81), nrow = 9) Mat2 <- Mat1 Mat1[1:9, 1:8] <- rseq Mat1 Mat2[1:9, 2:9] <- rseq Mat2 # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) (VFdf1 <- DispFieldbb(rast1, rast2, 2, 8, 2, 8)) # The second raster is shifted right by 1 unit relative to the first raster # dispx = 1
rseq <- stats::runif(72) Mat1 <- matrix(rep(0, 81), nrow = 9) Mat2 <- Mat1 Mat1[1:9, 1:8] <- rseq Mat1 Mat2[1:9, 2:9] <- rseq Mat2 # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) (VFdf1 <- DispFieldbb(rast1, rast2, 2, 8, 2, 8)) # The second raster is shifted right by 1 unit relative to the first raster # dispx = 1
This is an implementation of a novel algorithm that differs from more
traditional digital image correlation implementations that are applied in the
DispField
and DispFieldbb
functions. The function
calculates a displacement field representing persistent movement based on the
cross-covariance in a raster stack (in this case a sequential series of
rasters) presumably representing spatial population abundance or density at
more than two different instances of time. If analysis is restricted to only
two time instances, DispField
is more appropriate.
DispFieldST(inputstack1, lag1, factv1, facth1, restricted = FALSE)
DispFieldST(inputstack1, lag1, factv1, facth1, restricted = FALSE)
inputstack1 |
a raster stack with each raster layer representing an instance of time. The raster stack should be organized such that the first raster in the stack is the first observed spatial dataset and time progresses forward with the third dimension index of the raster stack. The raster stack should contain only numeric values. Any NA value will be converted to a zero |
lag1 |
an integer time lag |
factv1 |
an odd integer for the vertical dimension of subgrids |
facth1 |
an odd integer for the horizontal dimension of subgrids |
restricted |
logical (TRUE or FALSE) |
The input rasters in the raster stack are first converted to equivalent matrices, which together represent a three-dimensional array with two spatial dimensions and one time dimension. The prescribed lag is applied to the three dimensional array derived from the raster stack by first producing two equivalent arrays and then removing appropriate numbers of layers from the top of one and the bottom of the other. These are referred to as unlagged and lagged spatiotemporal arrays in the description that follows.
Prior to computing displacement based on direction of maximum cross-covariance, the function divides the spatial domain up into sub-grids of size factv1 X facth1, which are vertical and horizontal sub-grid spatial dimensions.
The function converts three dimensional lagged and unlagged spatiotemporal arrays to two-dimensional lagged and unlagged spatiotemporal matrices by averaging along one of the spatial dimensions (either rows or columns) to obtain two pairs of two-dimensional matrices in which one dimension is spatial (either rows or columns) and one dimension is temporal. One of each pair corresponds to the unlagged spatiotemporal array and the other corresponds to the lagged spatiotemporal array. Displacement in the vertical direction is computed using unlagged and lagged matrices that have been averaged along rows and displacement in the horizontal direction is computed using unlagged and lagged matrices that have been averaged along columns.
If restricted is set to FALSE (the default), the function computes cross-covariance between each sub-grid of the unlagged row-averaged spatiotemporal matrix and the whole row-averaged lagged spatiotemporal matrix and between each sub-grid of the unlagged column-averaged spatiotemporal matrix and the entirety corresponding lagged matrix.
If restricted is set to TRUE, the function uses cross-covariance between lagged and unlagged version of row-averaged and column averaged spatiotemporal matrices that have all been either row or column-averaged within sub-grids to estimate vertical and horizontal displacement.
Regardless of whether restricted is set to TRUE or FALSE, for each sub-grid, displacement in the x and y direction is divided by the shift in the time dimension to produce orthogonal velocity vetors. Note that for this reason, the lag1 argument of the function does not necessarily determine the time lag that is used to produce each orthoganal velocity vector.
Reference coordinates and cell size are extracted from the first raster stack such that the locations from whence displacement is estimated as well as displacement (or velocity) estimates can be expressed in the units of the projected coordinates.
The coordinates are assumed to increase vertically and horizontally from the lower left corner of the two-dimensional domain.
Caution is warranted when defining the sub-grid dimensions because the function can produce erroneous results when sub-grids are too small.
In addition, results can be quite sensitive to specification of the time lag.
If velocities are highly variable in space or over time, avoid specifying a
single time lag by calling the related DispFieldSTall
function.
A data frame is returned with the following column names: rowcent, colcent, frowmin, frowmax, fcolmin, fcolmax, centx, centy, dispx, and dispy. The rowcent and colcent column names are the row and column indices for the center of each sub-grid; frowmin and frowmax are the sub-grid minimum and maximum row indices; fcolmin and fcolmax are the sub-grid minimum and maximum column indices; centx and centy are the projected coordinates of the centre of the subgrid derived from the raster input files; dispx and dispy are the orthoganal velocity vectors in units of space per timestep in the horizontal and vertical directions in the same spatial units as the projected coordinates of the raster input files.
DispField
for a similar function with a grid of focal
regions for only two time instances, DispFieldSTbb
for a
version designed to quantify persistent directional movement when the time
series features more than two time instances but using a bounding pox to
define a focal region, see DispFieldSTall
for a version
designed to quantify persistent directional movement when velocity is
variable in space, and Xcov2D
for demonstration of how
two-dimensional cross-covariance is used to determine displacement (see
examples of Xcov2D function documentation).
(Mat1 <- matrix(rep(c(1:5, 0, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat2 <- matrix(rep(c(0, 1:5, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat3 <- matrix(rep(c(0, 0, 1:5, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat4 <- matrix(rep(c(0, 0, 0, 1:5, 0), 9), nrow = 9, byrow = TRUE)) # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) rast3 <- terra::rast(Mat3) terra::plot(rast3) rast4 <- terra::rast(Mat4) terra::plot(rast4) teststack1 <- c(rast1, rast2, rast3, rast4) (VFdf2 <- DispFieldST(teststack1, lag1 = 1, factv1 = 9, facth1 = 9)) # block is moving rightward at a speed of 1 unit of space per unit of time # dispx = 1
(Mat1 <- matrix(rep(c(1:5, 0, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat2 <- matrix(rep(c(0, 1:5, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat3 <- matrix(rep(c(0, 0, 1:5, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat4 <- matrix(rep(c(0, 0, 0, 1:5, 0), 9), nrow = 9, byrow = TRUE)) # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) rast3 <- terra::rast(Mat3) terra::plot(rast3) rast4 <- terra::rast(Mat4) terra::plot(rast4) teststack1 <- c(rast1, rast2, rast3, rast4) (VFdf2 <- DispFieldST(teststack1, lag1 = 1, factv1 = 9, facth1 = 9)) # block is moving rightward at a speed of 1 unit of space per unit of time # dispx = 1
This is an implementation of a novel algorithm that differs from more
traditional digital image correlation implementations that are applied in the
DispField
and DispFieldbb
functions. This version
is similar to the DispFieldST
function except that it does not
require a specific time lag. Instead the user specifies a maximum time lag
and the function computes displacement vectors using the time lag that
produces the maximum speed (magnitude of displacement divided by time lag).
The function calculates a displacement field representing persistent movement
based on the cross-covariance in a raster stack (in this case a sequential
series of rasters) presumably representing spatial population abundance or
density at more than two different instances of time. If analysis is
restricted to only two time instances, DispField
is more
appropriate.
DispFieldSTall(inputstack1, lagmax, factv1, facth1, restricted = FALSE)
DispFieldSTall(inputstack1, lagmax, factv1, facth1, restricted = FALSE)
inputstack1 |
a raster stack with each raster layer representing an instance of time. The raster stack should be organized such that the first raster in the stack is the first observed spatial dataset and time progresses forward with the third dimension index of the raster stack. The raster stack should contain only numeric values. Any NA value will be converted to a zero |
lagmax |
an integer representing the maximum time lag |
factv1 |
an odd integer for the vertical dimension of subgrids |
facth1 |
an odd integer for the horizontal dimension of subgrids |
restricted |
logical (TRUE or FALSE) |
The DispFieldSTall function has the same inner workings as the
DispFieldST
function except that instead of specifying a
specific time lag, the user specifies a maximum time lag. The function then
cycles through all lags up to the maximum time lag and chooses the for each
location the maximum speed. The DispFieldSTall function is more appropriate
than DispFieldST
when velocity is variable in space.
Caution is warranted when defining the sub-grid dimensions because the function can produce erroneous results when sub-grids are too small.
A data frame is returned with the following column names: rowcent, colcent, frowmin, frowmax, fcolmin, fcolmax, centx, centy, dispx, and dispy. The rowcent and colcent column names are the row and column indices for the center of each sub-grid; frowmin and frowmax are the sub-grid minimum and maximum row indices; fcolmin and fcolmax are the sub-grid minimum and maximum column indices; centx and centy are the projected coordinates of the centre of the subgrid derived from the raster input files; dispx and dispy are the orthoganal velocity vectors in units of space per timestep in the horizontal and vertical directions in the same spatial units as the projected coordinates of the raster input files.
DispField
for a similar function with a grid of focal
regions for only two time instances, DispFieldST
for a
version designed to quantify persistent directional movement when the time
series features more than two time instances and the velocity is constant
in space, DispFieldSTbball
for a version designed to quantify
persistent directional movement when velocity is variable in space and the
focal region is defined using a bounding box, and Xcov2D
for
demonstration of how two-dimensional cross-covariance is used to determine
displacement (see examples of Xcov2D function documentation).
(Mat1 <- matrix(rep(c(1:5, 0, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat2 <- matrix(rep(c(0, 1:5, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat3 <- matrix(rep(c(0, 0, 1:5, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat4 <- matrix(rep(c(0, 0, 0, 1:5, 0), 9), nrow = 9, byrow = TRUE)) # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) rast3 <- terra::rast(Mat3) terra::plot(rast3) rast4 <- terra::rast(Mat4) terra::plot(rast4) teststack1 <- c(rast1, rast2, rast3, rast4) (VFdf4 <- DispFieldSTall(teststack1, lagmax = 2, factv1 = 9, facth1 = 9)) # block is moving rightward at a speed of 1 unit of space per unit of time # dispx = 1
(Mat1 <- matrix(rep(c(1:5, 0, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat2 <- matrix(rep(c(0, 1:5, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat3 <- matrix(rep(c(0, 0, 1:5, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat4 <- matrix(rep(c(0, 0, 0, 1:5, 0), 9), nrow = 9, byrow = TRUE)) # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) rast3 <- terra::rast(Mat3) terra::plot(rast3) rast4 <- terra::rast(Mat4) terra::plot(rast4) teststack1 <- c(rast1, rast2, rast3, rast4) (VFdf4 <- DispFieldSTall(teststack1, lagmax = 2, factv1 = 9, facth1 = 9)) # block is moving rightward at a speed of 1 unit of space per unit of time # dispx = 1
This is an implementation of a novel algorithm that differs from more
traditional digital image correlation implementations that are applied in the
DispField
and DispFieldbb
functions. The function
calculates a displacement field representing persistent movement based on the
cross-covariance in a raster stack (in this case a sequential series of
rasters) presumably representing spatial population abundance or density at
more than two different instances of time. If analysis is restricted to only
two time instances, DispFieldbb
is more appropriate.
DispFieldSTbb( inputstack1, lag1, rowmn, rowmx, colmn, colmx, restricted = FALSE )
DispFieldSTbb( inputstack1, lag1, rowmn, rowmx, colmn, colmx, restricted = FALSE )
inputstack1 |
a raster stack with each raster layer representing an instance of time. The raster stack should be organized such that the first raster in the stack is the first observed spatial dataset and time progresses forward with the third dimension index of the raster stack. The raster stack should contain only numeric values. Any NA value will be converted to a zero |
lag1 |
an integer time lag |
rowmn |
an integer for the minimum row index of the bounding box |
rowmx |
an integer for the maximum row index of the bounding box |
colmn |
an integer for the minimum column index of the bounding box |
colmx |
an integer for the maximum column index of the bounding box |
restricted |
logical (TRUE or FALSE) |
The input rasters in the raster stack are first converted to equivalent matrices, which together represent a three-dimensional array with two spatial dimensions and one time dimension. The prescribed lag is applied to the three dimensional array derived from the raster stack by first producing two equivalent arrays and then removing appropriate numbers of layers from the top of one and the bottom of the other. These are referred to as unlagged and lagged spatiotemporal arrays in the description that follows.
The function converts three dimensional lagged and unlagged spatiotemporal arrays to two-dimensional lagged and unlagged spatiotemporal matrices by averaging along one of the spatial dimensions (either rows or columns) to obtain two pairs of two-dimensional matrices in which one dimension is spatial (either rows or columns) and one dimension is temporal. One of each pair corresponds to the unlagged spatiotemporal array and the other corresponds to the lagged spatiotemporal array. Displacement in the vertical direction is computed using unlagged and lagged matrices that have been averaged along rows and displacement in the horizontal direction is computed using unlagged and lagged matrices that have been averaged along columns.
If restricted is set to FALSE (the default), the function computes cross-covariance between the values within the bounding box of the unlagged row-averaged spatiotemporal matrix and the whole row-averaged lagged spatiotemporal matrix and between the values within the bounding box of the unlagged column-averaged spatiotemporal matrix and the entirety corresponding lagged matrix.
If restricted is set to TRUE, the function uses cross-covariance between lagged and unlagged version of row-averaged and column averaged spatiotemporal matrices that have all been either row or column-averaged within the bounding box to estimate vertical and horizontal displacement.
Regardless of whether restricted is set to TRUE or FALSE, for each sub-grid, displacement in the x and y direction is divided by the shift in the time dimension to produce orthogonal velocity vetors. Note that for this reason, the lag1 argument of the function does not necessarily determine the time lag that is used to produce each orthoganal velocity vector.
Reference coordinates and cell size are extracted from the first raster stack such that the locations from whence displacement is estimated as well as displacement (or velocity) estimates can be expressed in the units of the projected coordinates.
The coordinates are assumed to increase vertically and horizontally from the lower left corner of the two-dimensional domain.
Caution is warranted when defining the sub-grid dimensions because the function can produce erroneous results when sub-grids are too small.
#' In addition, results can be quite sensitive to specification of the time
lag. If velocities are highly variable in space or over time, avoid
specifying a single time lag by calling the related
DispFieldSTbball
function.
A data frame is returned with the following column names: rowcent, colcent, frowmin, frowmax, fcolmin, fcolmax, centx, centy, dispx, and dispy. The rowcent and colcent column names are the row and column indices for the center of each sub-grid; frowmin and frowmax are the sub-grid minimum and maximum row indices; fcolmin and fcolmax are the sub-grid minimum and maximum column indices; centx and centy are the projected coordinates of the centre of the subgrid derived from the raster input files; dispx and dispy are the orthoganal velocity vectors in units of space per timestep in the horizontal and vertical directions in the same spatial units as the projected coordinates of the raster input files.
DispField
for a similar function with a grid of focal
regions for only two time instances, DispFieldST
for a
version designed to quantify persistent directional movement when the time
series features more than two time instances but using a grid to define
focal regions, see DispFieldSTbball
for a version designed to
quantify persistent directional movement when velocity is variable in
space, and Xcov2D
for demonstration of how two-dimensional
cross-covariance is used to determine displacement (see examples of Xcov2D
function documentation).
rseq <- stats::runif(54) Mat1 <- matrix(rep(0, 9*9), nrow = 9) Mat2 <- Mat1; Mat3 <- Mat1; Mat4 <- Mat1 Mat1[1:9, 1:6] <- rseq Mat1 Mat2[1:9, 2:7] <- rseq Mat2 Mat3[1:9, 3:8] <- rseq Mat3 Mat4[1:9, 4:9] <- rseq Mat4 # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) rast3 <- terra::rast(Mat3) terra::plot(rast3) rast4 <- terra::rast(Mat4) terra::plot(rast4) teststack1 <- c(rast1, rast2, rast3, rast4) (VFdf3 <- DispFieldSTbb(teststack1, lag1 = 1, 2, 8, 2, 8)) # block is moving rightward at a speed of 1 unit of space per unit of time # dispx = 1
rseq <- stats::runif(54) Mat1 <- matrix(rep(0, 9*9), nrow = 9) Mat2 <- Mat1; Mat3 <- Mat1; Mat4 <- Mat1 Mat1[1:9, 1:6] <- rseq Mat1 Mat2[1:9, 2:7] <- rseq Mat2 Mat3[1:9, 3:8] <- rseq Mat3 Mat4[1:9, 4:9] <- rseq Mat4 # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) rast3 <- terra::rast(Mat3) terra::plot(rast3) rast4 <- terra::rast(Mat4) terra::plot(rast4) teststack1 <- c(rast1, rast2, rast3, rast4) (VFdf3 <- DispFieldSTbb(teststack1, lag1 = 1, 2, 8, 2, 8)) # block is moving rightward at a speed of 1 unit of space per unit of time # dispx = 1
This is an implementation of a novel algorithm that differs from more
traditional digital image correlation implementations that are applied in the
DispField
and DispFieldbb
functions. This version
is similar to the DispFieldSTbb
function except that it does
not require a specific time lag. Instead the user specifies a maximum time
lag and the function computes displacement vectors using the time lag that
produces the maximum speed (magnitude of displacement divided by time lag).
The function calculates a displacement field representing persistent movement
based on the cross-covariance in a raster stack (in this case a sequential
series of rasters) presumably representing spatial population abundance or
density at more than two different instances of time. If analysis is
restricted to only two time instances, DispFieldbb
is more
appropriate.
DispFieldSTbball( inputstack1, lagmax, rowmn, rowmx, colmn, colmx, restricted = FALSE )
DispFieldSTbball( inputstack1, lagmax, rowmn, rowmx, colmn, colmx, restricted = FALSE )
inputstack1 |
a raster stack with each raster layer representing an instance of time. The raster stack should be organized such that the first raster in the stack is the first observed spatial dataset and time progresses forward with the third dimension index of the raster stack. The raster stack should contain only numeric values. Any NA value will be converted to a zero |
lagmax |
an integer representing the maximum time lag |
rowmn |
an integer denoting the minimum row index of the sub-grid |
rowmx |
an integer denoting the maximum row index of the sub-grid |
colmn |
an integer denoting the minimum column index of the sub-grid |
colmx |
an integer denoting the maximum column index of the sub-grid |
restricted |
logical (TRUE or FALSE) |
The DispFieldSTbball function has the same inner workings as the
DispFieldSTbb
function except that instead of specifying a
specific time lag, the user specifies a maximum time lag. The function then
cycles through all lags up to the maximum time lag and chooses the for each
location the maximum speed. The DispFieldSTbball function is more appropriate
than DispFieldSTbb
when velocity is variable in space.
Caution is warranted when defining the bounding box dimensions because the function can produce erroneous results when the bounding box is too small.
A data frame is returned with the following column names: rowcent, colcent, frowmin, frowmax, fcolmin, fcolmax, centx, centy, dispx, and dispy. The rowcent and colcent column names are the row and column indices for the center of each sub-grid; frowmin and frowmax are the sub-grid minimum and maximum row indices; fcolmin and fcolmax are the sub-grid minimum and maximum column indices; centx and centy are the projected coordinates of the centre of the subgrid derived from the raster input files; dispx and dispy are the orthoganal velocity vectors in units of space per timestep in the horizontal and vertical directions in the same spatial units as the projected coordinates of the raster input files.
DispFieldbb
for a similar function with focal region
defined using a bounding box for only two time instances,
DispFieldSTbb
for a version designed to quantify persistent
directional movement when velocity is constant in space and the focal
region is defined using a bounding box, see DispFieldSTall
for a version designed to quantify persistent directional movement when
velocity is variable in space and focal regions are defined based on a
grid, and Xcov2D
for demonstration of how two-dimensional
cross-covariance is used to determine displacement (see examples of Xcov2D
function documentation).
(Mat1 <- matrix(rep(c(1:5, 0, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat2 <- matrix(rep(c(0, 1:5, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat3 <- matrix(rep(c(0, 0, 1:5, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat4 <- matrix(rep(c(0, 0, 0, 1:5, 0), 9), nrow = 9, byrow = TRUE)) # Rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) rast3 <- terra::rast(Mat3) terra::plot(rast3) rast4 <- terra::rast(Mat4) terra::plot(rast4) teststack1 <- c(rast1, rast2, rast3, rast4) (VFdf5 <- DispFieldSTbball(teststack1, lagmax = 2, 1, 9, 1, 9)) # block is moving rightward at a speed of 1 unit of space per unit of time # dispx = 1
(Mat1 <- matrix(rep(c(1:5, 0, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat2 <- matrix(rep(c(0, 1:5, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat3 <- matrix(rep(c(0, 0, 1:5, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat4 <- matrix(rep(c(0, 0, 0, 1:5, 0), 9), nrow = 9, byrow = TRUE)) # Rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) rast3 <- terra::rast(Mat3) terra::plot(rast3) rast4 <- terra::rast(Mat4) terra::plot(rast4) teststack1 <- c(rast1, rast2, rast3, rast4) (VFdf5 <- DispFieldSTbball(teststack1, lagmax = 2, 1, 9, 1, 9)) # block is moving rightward at a speed of 1 unit of space per unit of time # dispx = 1
Functions for computing the statistics which may be driving variables of
movement that has been quantified using the DispField
or
DispFieldbb
functions. The same raster data as were supplied to
the aforementioned functions must be supplied to these in addition to a
raster layer for which statistics are sought. Then for each region of
interest defined when DispField
or DispFieldbb
were called, these functions compute statistics for presumed source (sourceloc =
TRUE) locations or presumed sink locations (sourceloc = FALSE). Note that in the
DispMornasI function, defining radius using distance means that a radius of
one corresponds to the rook's neighbourhood.
DispMoransI(inputrast1, inputrast2, statrast, vfdf, sourceloc = TRUE, rad1) DispStats( inputrast1, inputrast2, statrast, vfdf, sourceloc = TRUE, statistic = "var" )
DispMoransI(inputrast1, inputrast2, statrast, vfdf, sourceloc = TRUE, rad1) DispStats( inputrast1, inputrast2, statrast, vfdf, sourceloc = TRUE, statistic = "var" )
inputrast1 |
a raster as produced by terra::rast |
inputrast2 |
a raster of equivalent dimension to inputrast1 as produced by terra::rast |
statrast |
a raster of equivalent dimension to inputrast1 as produced by terra::rast which contains the variable that will be used to compute statistics |
vfdf |
a data frame returned by the |
sourceloc |
logical (TRUE or FALSE) indicating whether statistics are to be returned at source or sink locations |
rad1 |
an ingeger indicating the neighbourhood radius for Moran's I statistic calculations in rows/columns. Any cell within a distance of rad1 cells of the focal cell is considered to be in its neighbourhood. |
statistic |
desired output statistic: It should be one of "mean", "var", or "sum". Default setting is var. |
A data frame is returned with all of the same columns as the vfdf input data frame plus an additional column containing the computed statistic in each region of interest defined in vfdf.
# Illustrating use of DispMoransI: (Mat1 <- matrix(c(0.1,1,0.1,0,0,0,0,0,0, 1,0.1,1,0,0,0,0,0,0, 0.1,1,0.1,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0), nrow = 9)) (Mat2 <- matrix(c(0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0.1,1,0.1,0,0,0,0,0,0, 1,0.1,1,0,0,0,0,0,0, 0.1,1,0.1,0,0,0,0,0,0), nrow = 9)) # Note that rasterizing a matrix causes it to be rotated 90 degrees. # Therefore, any shift in the x direction is in fact now a shift in the y direction rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) (VFdf1 <- DispField(rast1, rast2, factv1 = 3, facth1 = 3)) # The second raster is shifted down by -0.6666667 units relative to the first raster # dispy = -0.6666667 (the width of each box is 0.1111111). # Now to compute the statistics at the source: the Moran's I of the original values # in each region of interest (should be minus one in first row) (VFdf2 <- DispMoransI(rast1, rast2, rast1, VFdf1, sourceloc = TRUE, rad1 = 1)) # Illustrating use of DispStats: (Mat1 <- matrix(c(1,1,1,0,0,0,0,0,0, 1,1,1,0,0,0,0,0,0, 1,1,1,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0), nrow = 9)) (Mat2 <- matrix(c(0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 1,1,1,0,0,0,0,0,0, 1,1,1,0,0,0,0,0,0, 1,1,1,0,0,0,0,0,0), nrow = 9)) # Note that rasterizing a matrix causes it to be rotated 90 degrees. # Therefore, any shift in the x direction is in fact now a shift in the y direction rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) (VFdf1 <- DispField(rast1, rast2, factv1 = 3, facth1 = 3)) # The second raster is shifted down by -0.6666667 units relative to the first raster # dispy = -0.6666667 (the width of each box is 0.1111111). # Now to compute the statistics at the source: the mean of the original values # in each region of interest (should be one in first row) (VFdf2 <- DispStats(rast1, rast2, rast1, VFdf1, sourceloc = TRUE, statistic = "mean")) # sum in each region of interest (should be nine in first row) (VFdf3 <- DispStats(rast1, rast2, rast1, VFdf1, sourceloc = TRUE, statistic = "sum")) # variance in each region of interest (should be zero in all rows) (VFdf4 <- DispStats(rast1, rast2, rast1, VFdf1, sourceloc = TRUE, statistic = "var"))
# Illustrating use of DispMoransI: (Mat1 <- matrix(c(0.1,1,0.1,0,0,0,0,0,0, 1,0.1,1,0,0,0,0,0,0, 0.1,1,0.1,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0), nrow = 9)) (Mat2 <- matrix(c(0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0.1,1,0.1,0,0,0,0,0,0, 1,0.1,1,0,0,0,0,0,0, 0.1,1,0.1,0,0,0,0,0,0), nrow = 9)) # Note that rasterizing a matrix causes it to be rotated 90 degrees. # Therefore, any shift in the x direction is in fact now a shift in the y direction rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) (VFdf1 <- DispField(rast1, rast2, factv1 = 3, facth1 = 3)) # The second raster is shifted down by -0.6666667 units relative to the first raster # dispy = -0.6666667 (the width of each box is 0.1111111). # Now to compute the statistics at the source: the Moran's I of the original values # in each region of interest (should be minus one in first row) (VFdf2 <- DispMoransI(rast1, rast2, rast1, VFdf1, sourceloc = TRUE, rad1 = 1)) # Illustrating use of DispStats: (Mat1 <- matrix(c(1,1,1,0,0,0,0,0,0, 1,1,1,0,0,0,0,0,0, 1,1,1,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0), nrow = 9)) (Mat2 <- matrix(c(0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 1,1,1,0,0,0,0,0,0, 1,1,1,0,0,0,0,0,0, 1,1,1,0,0,0,0,0,0), nrow = 9)) # Note that rasterizing a matrix causes it to be rotated 90 degrees. # Therefore, any shift in the x direction is in fact now a shift in the y direction rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) (VFdf1 <- DispField(rast1, rast2, factv1 = 3, facth1 = 3)) # The second raster is shifted down by -0.6666667 units relative to the first raster # dispy = -0.6666667 (the width of each box is 0.1111111). # Now to compute the statistics at the source: the mean of the original values # in each region of interest (should be one in first row) (VFdf2 <- DispStats(rast1, rast2, rast1, VFdf1, sourceloc = TRUE, statistic = "mean")) # sum in each region of interest (should be nine in first row) (VFdf3 <- DispStats(rast1, rast2, rast1, VFdf1, sourceloc = TRUE, statistic = "sum")) # variance in each region of interest (should be zero in all rows) (VFdf4 <- DispStats(rast1, rast2, rast1, VFdf1, sourceloc = TRUE, statistic = "var"))
Here is a function that will find the row and column indices of a matrix that are associated with a vector index.
GetRowCol(Index, dim1, dim2)
GetRowCol(Index, dim1, dim2)
Index |
an integer vector index |
dim1 |
integer row dimension of the matrix from which the row and column indices are to be extracted |
dim2 |
integer column dimension of the matrix from which the row and column indices are to be extracted |
Often when applying functions like the R function which.max(matrix) to a matrix, a vector index is returned when the coder would prefer to have a row and column indices. This function converts the vector index to row and column indices.
The function assumes that the elements of the matrix are filled by column (byrow = FALSE), which is the default R matrix behaviour.
a numeric vector of length two with two integers indicating row and column respectively
GetRowCol(6, dim1 = 3, dim2 = 3) # should return c(3, 2)
GetRowCol(6, dim1 = 3, dim2 = 3) # should return c(3, 2)
Compute Moran's I for a matrix. A fast implementation of Moran's I for gridded data, with neighbours defined based on a radial distance. Note that when using radius to define the neighbourhood, a radius of one corresponds to the rook's neibhourhood. There is currently no equivalent to queen's neighbourhood.
MoransI(mat1, r1)
MoransI(mat1, r1)
mat1 |
a matrix of values; NA/Inf values must be coded as NA and are ignored |
r1 |
an integer representing the distance (radius), within which nearby cells are considered neighbours in units of rows/columns |
a single numeric value for Moran's I
(TestMat <- matrix(c(1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1), nrow =5)) # the code below should return -1 MoransI(TestMat, r1 = 1)
(TestMat <- matrix(c(1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1), nrow =5)) # the code below should return -1 MoransI(TestMat, r1 = 1)
Detect patterns in vector fields represented on a grid by looking in the rook's neighbourhood of each grid cell. Four patterns are detected: convergences occur when the vectors in the four adjacent cells in the rook's neighbourhood point towards the focal cell; divergences occur when the vectors in the four adjacent cells point away from the focal cell; Partial convergences occur when three of the four vectors point towards the focal cell and the final vector points neither towards nor away from the focal cell; Partial divergences occur when three of the four vectors point away the focal grid cell and the final vector points neither towards nor away from the focal grid. For all of the patterns above a sub-pattern is specified if all arrows point clockwise or counter-clockwise.
PatternDetect(vfdf)
PatternDetect(vfdf)
vfdf |
A data frame as returned by |
A data frame as returned by DispField
,
DispFieldST
, or DispFieldSTall
, with three
additional columns. The first additional column is called Pattern in which
the patterns around each focal cell are categorized as convergence,
divergence, partial convergence, partial divergence, or NA. The second
additional column, called SubPattern, indicates whether all arrows point
clockwise or counter-clockwise. The third additional column is called
PatternCt, which contains a one if all four neighbourhood grid cells
contain displacement estimates, and a NA otherwise.
# creating convergence/divergence patterns Mat1 <- matrix(rep(0,9*9), nrow = 9) Mat1[3, c(4, 6)] <- 1 Mat1[7, c(4, 6)] <- 1 Mat1[c(4, 6), 3] <- 1 Mat1[c(4, 6), 7] <- 1 Mat1 Mat2 <- matrix(rep(0,9*9), nrow = 9) Mat2[2, c(4, 6)] <- 1 Mat2[8, c(4, 6)] <- 1 Mat2[c(4, 6), 2] <- 1 Mat2[c(4, 6), 8] <- 1 Mat2 # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) # Detecting a divergence (VFdf1 <- DispField(rast1, rast2, factv1 = 3, facth1 = 3, restricted = TRUE)) (patdf1 <- PatternDetect(VFdf1)) # Detecting a convergence (VFdf2 <- DispField(rast2, rast1, factv1 = 3, facth1 = 3, restricted = TRUE)) (patdf2 <- PatternDetect(VFdf2))
# creating convergence/divergence patterns Mat1 <- matrix(rep(0,9*9), nrow = 9) Mat1[3, c(4, 6)] <- 1 Mat1[7, c(4, 6)] <- 1 Mat1[c(4, 6), 3] <- 1 Mat1[c(4, 6), 7] <- 1 Mat1 Mat2 <- matrix(rep(0,9*9), nrow = 9) Mat2[2, c(4, 6)] <- 1 Mat2[8, c(4, 6)] <- 1 Mat2[c(4, 6), 2] <- 1 Mat2[c(4, 6), 8] <- 1 Mat2 # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) # Detecting a divergence (VFdf1 <- DispField(rast1, rast2, factv1 = 3, facth1 = 3, restricted = TRUE)) (patdf1 <- PatternDetect(VFdf1)) # Detecting a convergence (VFdf2 <- DispField(rast2, rast1, factv1 = 3, facth1 = 3, restricted = TRUE)) (patdf2 <- PatternDetect(VFdf2))
In order to produce reliable results, the cross-covariance approach
implemented in DispField
. DispFieldST
, and
DispFieldSTall
needs a certain minimum number of non-zero or
non-NA valued pixels in pairs of images or pairs of arrays derived from a
raster stack that it uses to compute cross-covariance. The user may define a
threshold such as ten percent of the pixels within each sub-grid. This
function allows the user to assess whether the minimum threshold number of
non-zero pixels per sub-grid are surpassed by returning the number of
non-zero pixels within each sub-grid over all of the time instances in the
user-supplied raster stack. The user may choose to disregard or mistrust
displacement or velocity estimates derived from sub-grids with insufficient
numbers of non-zero pixels. This function is designed to be called before or
after one of the functions referenced above in order to enable the user to
quantify their confidence in displacement or velocity estimates.
PixelCt(inputstack1, factv1, facth1)
PixelCt(inputstack1, factv1, facth1)
inputstack1 |
a raster stack with each raster layer representing an instance of time. The raster stack should be organized such that the first raster in the stack is the first observed spatial dataset and time progresses forward with the third dimension index of the raster stack. The raster stack should contain only numeric values. Any NA value will be converted to a zero |
factv1 |
an odd integer for the vertical dimension of subgrids |
facth1 |
an odd integer for the horizontal dimension of subgrids |
A data frame is returned with the following column names: rowcent, colcent, frowmin, frowmax, fcolmin, fcolmax, and PixelCt. The rowcent and colcent column names are the row and column indices for the center of each sub-grid; frowmin and frowmax are the sub-grid minimum and maximum row indices; fcolmin and fcolmax are the sub-grid minimum and maximum column indices; pixelct is the count of non-zero pixels in the sub-grid over the entire time period covered by the input raster stack.
# below the example in the DispField function documentation is reproduced (Mat1 <- matrix(rep(c(1:5, 0, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat2 <- matrix(rep(c(0, 1:5, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) # Note that rasterizing a matrix causes it to be rotated 90 degrees. # Therefore, any shift in the x direction is in fact now a shift in the y direction rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) (Confdf1 <- PixelCt(c(rast1, rast2), factv1 = 9, facth1 = 9)) # This should return a pixel count of 54: This is the number # of pixels that were occupied in either the first or second # time instance.
# below the example in the DispField function documentation is reproduced (Mat1 <- matrix(rep(c(1:5, 0, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) (Mat2 <- matrix(rep(c(0, 1:5, 0, 0, 0), 9), nrow = 9, byrow = TRUE)) # Note that rasterizing a matrix causes it to be rotated 90 degrees. # Therefore, any shift in the x direction is in fact now a shift in the y direction rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) (Confdf1 <- PixelCt(c(rast1, rast2), factv1 = 9, facth1 = 9)) # This should return a pixel count of 54: This is the number # of pixels that were occupied in either the first or second # time instance.
This function converts the data that accompany the ICvectorfields R package
to a raster stack. The raster stack is the only accepted data input format
for the following ICvectorfields functions: DispFieldST
,
DispFieldSTbb
, DispFieldSTall
,
DispFieldSTbball
.
RastStackData(inputdf)
RastStackData(inputdf)
inputdf |
a data frame object in which the first column is longitude (or x coordinate), the second column is latitude (or y coordinate), and all of the subsequent columns represent a measure of population abundance or density at a unique instance of time. Each row of the input data frame, therefore, represents a unique spatial location, which should be on an evenly spaced grid. Note, however, that not all grid locations need to have observations; some grid locations can have values of NA or can be missing entirely. |
Once a raster stack has been created, individual layers can be subsetted using rasterstack[[index]], where index is an integer index for the third dimension of the raster stack.
The function returns a raster stack constructed using inputdf. Each layer in the stack corresponds to a column of the input dataset (after the first two columns, which are longitude and latitude). The extent of all of the rasters in the stack is constructed using the minimum and maximum longitudes and latitudes.
# creating random data in the correct data format xyzdf <- expand.grid(x = c(1:3), y = c(1:3)) xyzdf$z1 <- runif(9) xyzdf$z2 <- runif(9) xyzdf$z3 <- runif(9) zstack <- RastStackData(xyzdf) dim(zstack) terra::plot(zstack[[1]]) terra::plot(zstack[[2]]) terra::plot(zstack[[3]])
# creating random data in the correct data format xyzdf <- expand.grid(x = c(1:3), y = c(1:3)) xyzdf$z1 <- runif(9) xyzdf$z2 <- runif(9) xyzdf$z3 <- runif(9) zstack <- RastStackData(xyzdf) dim(zstack) terra::plot(zstack[[1]]) terra::plot(zstack[[2]]) terra::plot(zstack[[3]])
The movement of populations into adjacent cells may sometimes be influenced
by the gradient of some predictive variable. This function enables the
calculation of a simple gradient statistic in the rook's neighbourhood of
each cell in a dataset. The statistic must first be computed for each grid
cell using SubgridStats
. Then for each grid, the RooksGradient
function computes the arithmetic average of the difference between the
statistic at the focal grid cell and the statistic in the four (or fewer)
adjacent neighbours in its Rook's neibourhood. This arithmetically averaged
difference is then returned under the column header of 'Gradient'. A negative
gradient estimate indicates that the statistic in the central cell is higher
than that in neighbouring cells whereas a positive gradient estimate
indicates the opposite.
RooksGradient(vfdf, statistic = "mean")
RooksGradient(vfdf, statistic = "mean")
vfdf |
A data frame as returned by |
statistic |
desired output statistic: It should be one of "mean", "var", or "sum". Default setting is mean. |
A data frame similar to vfdf except that it includes an additional column called Gradient as described above.
# creating pattern patterns Mat1 <- matrix(rep(0,9*9), nrow = 9) Mat1[c(4:6), c(4:6)] <- 2 Mat1[c(4:6), c(1:3)] <- 1 Mat1[c(1:3), c(4:6)] <- 1 Mat1[c(7:9), c(4:6)] <- 1 Mat1[c(4:6), c(7:9)] <- 1 Mat1 Rast1 <- terra::rast(Mat1) terra::plot(Rast1) # calculating the mean in 9 subgrids (statsdf1 <- SubgridStats(Rast1, factv1 = 3, facth1 = 3, statistic = "mean")) # computing the gradient statistic on the mean (graddf1 <- RooksGradient(statsdf1, statistic = "mean")) # the Gradient statistic in the central grid in row 5 should # be equal to negative one
# creating pattern patterns Mat1 <- matrix(rep(0,9*9), nrow = 9) Mat1[c(4:6), c(4:6)] <- 2 Mat1[c(4:6), c(1:3)] <- 1 Mat1[c(1:3), c(4:6)] <- 1 Mat1[c(7:9), c(4:6)] <- 1 Mat1[c(4:6), c(7:9)] <- 1 Mat1 Rast1 <- terra::rast(Mat1) terra::plot(Rast1) # calculating the mean in 9 subgrids (statsdf1 <- SubgridStats(Rast1, factv1 = 3, facth1 = 3, statistic = "mean")) # computing the gradient statistic on the mean (graddf1 <- RooksGradient(statsdf1, statistic = "mean")) # the Gradient statistic in the central grid in row 5 should # be equal to negative one
This function prunes the data frame returned by the
PatternDetect
function such that it includes only rook's
neighborhoods that do not overlap. Pruning is done by sequential removal of
observations that are too near one another and as a result are overlapping.
Locations that are the most highly connected are removed first.
RooksNeighCt(vfdf)
RooksNeighCt(vfdf)
vfdf |
A data frame as returned by |
The reason for this function's existence is to facilitate probabilistic calculations regarding whether certain patterns are occurring more or less often that would be expected by chance. If rook's neighborhoods in which patterns are observed overlap, then the assumption of probabilistic independence is necessarily incorrect. Thus, overlap invalidates any calculation of the probability of occurrence of a particular pattern if that calculation assumes independence. The pruning actions of this function enable the user to more safely assume probabilistic independence.
A data frame similar to vfdf except that it includes only grid locations with speed estimates in all four adjacent grid locations in their rook's neighborhood. An additional column called IndPatternCt is appended which contains NA values for locations that are overlapping other locations and ones for all non-overlapping locations that have speed estimates in all four adjacent cells.
# creating convergence/divergence patterns Mat1 <- matrix(rep(0,9*9), nrow = 9) Mat1[3, c(4, 6)] <- 1 Mat1[7, c(4, 6)] <- 1 Mat1[c(4, 6), 3] <- 1 Mat1[c(4, 6), 7] <- 1 Mat1 Mat2 <- matrix(rep(0,9*9), nrow = 9) Mat2[2, c(4, 6)] <- 1 Mat2[8, c(4, 6)] <- 1 Mat2[c(4, 6), 2] <- 1 Mat2[c(4, 6), 8] <- 1 Mat2 Mat1 <- cbind(Mat1, Mat1) Mat1 <- rbind(Mat1, Mat1) Mat2 <- cbind(Mat2, Mat2) Mat2 <- rbind(Mat2, Mat2) # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) # Detecting a divergence (VFdf1 <- DispField(rast1, rast2, factv1 = 3, facth1 = 3, restricted = TRUE)) (patdf1 <- PatternDetect(VFdf1)) (subdf1 <- RooksNeighCt(patdf1)) # The last call should print a data table with four rows, each with a one under # the column header of IndPatternCt
# creating convergence/divergence patterns Mat1 <- matrix(rep(0,9*9), nrow = 9) Mat1[3, c(4, 6)] <- 1 Mat1[7, c(4, 6)] <- 1 Mat1[c(4, 6), 3] <- 1 Mat1[c(4, 6), 7] <- 1 Mat1 Mat2 <- matrix(rep(0,9*9), nrow = 9) Mat2[2, c(4, 6)] <- 1 Mat2[8, c(4, 6)] <- 1 Mat2[c(4, 6), 2] <- 1 Mat2[c(4, 6), 8] <- 1 Mat2 Mat1 <- cbind(Mat1, Mat1) Mat1 <- rbind(Mat1, Mat1) Mat2 <- cbind(Mat2, Mat2) Mat2 <- rbind(Mat2, Mat2) # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) # Detecting a divergence (VFdf1 <- DispField(rast1, rast2, factv1 = 3, facth1 = 3, restricted = TRUE)) (patdf1 <- PatternDetect(VFdf1)) (subdf1 <- RooksNeighCt(patdf1)) # The last call should print a data table with four rows, each with a one under # the column header of IndPatternCt
After running the PatternDetect
function, this function enables
classification of neighbour cells. Because the pattern detect function
classifies central cells according to the patterns of vector direction in
their Rook's neighbourhood, the neighbouring grid locations that comprise the
pattern are not labeled. This is remedied by the RooksNeighFind function
which classifies neighbour cells around focal grids classified with one of
the four patterns that the PatternDetect
function is able to
recognize. The function returns a data frame similar to the input data frame
with a column appended. In the appended column, neighbours surrounding focal
cells labeled with a particular pattern will be labeled as follows:
neighbours of the divergence pattern are labeled with a one, neighbours of
the convergence pattern are labeled with a two, neighbours of the partial
divergence pattern are labeled with a three, and neighbours of the partial
convergence pattern are labeled with a four. In cases where neighbours are
shared, the priority order from lowest to highest is four for partial
convergence to one for divergence. Thus, a neighbour that is shared between a
focal grid classified as a convergence and a nearby focal grid classified as
a divergence will be labeled with a one instead of a two.
RooksNeighFind(vfdf)
RooksNeighFind(vfdf)
vfdf |
A data frame as returned by |
A data frame similar to vfdf except that it includes an additional column called NeighType as described above.
# creating convergence/divergence patterns Mat1 <- matrix(rep(0,9*9), nrow = 9) Mat1[3, c(4, 6)] <- 1 Mat1[7, c(4, 6)] <- 1 Mat1[c(4, 6), 3] <- 1 Mat1[c(4, 6), 7] <- 1 Mat1 Mat2 <- matrix(rep(0,9*9), nrow = 9) Mat2[2, c(4, 6)] <- 1 Mat2[8, c(4, 6)] <- 1 Mat2[c(4, 6), 2] <- 1 Mat2[c(4, 6), 8] <- 1 Mat2 # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) # Detecting a divergence (VFdf1 <- DispField(rast1, rast2, factv1 = 3, facth1 = 3, restricted = TRUE)) (patdf1 <- PatternDetect(VFdf1)) (neighdf1 <- RooksNeighFind(patdf1)) # Rook's neighbour grids are labeled with a one. # Detecting a convergence (VFdf2 <- DispField(rast2, rast1, factv1 = 3, facth1 = 3, restricted = TRUE)) (patdf2 <- PatternDetect(VFdf2)) (neighdf2 <- RooksNeighFind(patdf2)) # Rook's neighbour grids are labeled with a two.
# creating convergence/divergence patterns Mat1 <- matrix(rep(0,9*9), nrow = 9) Mat1[3, c(4, 6)] <- 1 Mat1[7, c(4, 6)] <- 1 Mat1[c(4, 6), 3] <- 1 Mat1[c(4, 6), 7] <- 1 Mat1 Mat2 <- matrix(rep(0,9*9), nrow = 9) Mat2[2, c(4, 6)] <- 1 Mat2[8, c(4, 6)] <- 1 Mat2[c(4, 6), 2] <- 1 Mat2[c(4, 6), 8] <- 1 Mat2 # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) # Detecting a divergence (VFdf1 <- DispField(rast1, rast2, factv1 = 3, facth1 = 3, restricted = TRUE)) (patdf1 <- PatternDetect(VFdf1)) (neighdf1 <- RooksNeighFind(patdf1)) # Rook's neighbour grids are labeled with a one. # Detecting a convergence (VFdf2 <- DispField(rast2, rast1, factv1 = 3, facth1 = 3, restricted = TRUE)) (patdf2 <- PatternDetect(VFdf2)) (neighdf2 <- RooksNeighFind(patdf2)) # Rook's neighbour grids are labeled with a two.
Detect patterns in vector fields represented on a grid by looking in the
rook's neighbourhood of each grid cell. This function is analogous to
PatternDetect
, except that it detects rotational patterns. Four
patterns are detected: clockwise rotation when rotation in all four
neighbbour grids appears clockwise, counter clockwise rotation when rotation
in all four neighbour grids appears counter-clockwise, and partial clockwise
and counter-clockwise rotation, when all but one of the four adjacent
neighbour cells has vectors that indicate rotation. For all of the patterns
above a sub-pattern is specified as convergence when all of the vectors in
the four adjacent grids point towards the focal cell or a divergence when all
of the vectors in the four adjacent grids point away from the focal cell.
RotationDetect(vfdf)
RotationDetect(vfdf)
vfdf |
A data frame as returned by |
A data frame as returned by DispField
,
DispFieldST
, or DispFieldSTall
, with three
additional columns. The first additional column is called Pattern in which
the patterns around each focal cell are categorized as clockwise,
counter-clockwise, partial clockwise, partial counter-clockwise, or NA. The
second additional column, called SubPattern, indicates whether all arrows
point towards (convergence) or away (divergence) from the focal cell. The
third additional column is called PatternCt, which contains a one if all
four neighbourhood grid cells contain displacement estimates, and a NA
otherwise.
# creating rotation patterns Mat1 <- matrix(rep(0,9*9), nrow = 9) Mat1[c(1:3), 4] <- 1 Mat1[c(7:9), 6] <- 1 Mat1[4, c(7:9)] <- 1 Mat1[6, c(1:3)] <- 1 Mat1 Mat2 <- matrix(rep(0,9*9), nrow = 9) Mat2[c(1:3), 5] <- 1 Mat2[c(7:9), 5] <- 1 Mat2[5, c(7:9)] <- 1 Mat2[5, c(1:3)] <- 1 Mat2 # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) # Detecting clockwise rotation (VFdf1 <- DispField(rast1, rast2, factv1 = 3, facth1 = 3, restricted = TRUE)) (patdf1 <- RotationDetect(VFdf1)) # Detecting counter-clockwise rotation (VFdf2 <- DispField(rast2, rast1, factv1 = 3, facth1 = 3, restricted = TRUE)) (patdf2 <- RotationDetect(VFdf2))
# creating rotation patterns Mat1 <- matrix(rep(0,9*9), nrow = 9) Mat1[c(1:3), 4] <- 1 Mat1[c(7:9), 6] <- 1 Mat1[4, c(7:9)] <- 1 Mat1[6, c(1:3)] <- 1 Mat1 Mat2 <- matrix(rep(0,9*9), nrow = 9) Mat2[c(1:3), 5] <- 1 Mat2[c(7:9), 5] <- 1 Mat2[5, c(7:9)] <- 1 Mat2[5, c(1:3)] <- 1 Mat2 # rasterizing rast1 <- terra::rast(Mat1) terra::plot(rast1) rast2 <- terra::rast(Mat2) terra::plot(rast2) # Detecting clockwise rotation (VFdf1 <- DispField(rast1, rast2, factv1 = 3, facth1 = 3, restricted = TRUE)) (patdf1 <- RotationDetect(VFdf1)) # Detecting counter-clockwise rotation (VFdf2 <- DispField(rast2, rast1, factv1 = 3, facth1 = 3, restricted = TRUE)) (patdf2 <- RotationDetect(VFdf2))
Data based on a partial differential equation were simulated using the ReacTran R package (see details).
data(SimData)
data(SimData)
A data-frame with 40804 rows and 8 columns.
in arbitrary units
in arbitrary units
concentration in arbitrary units at t = 1
concentration in arbitrary units at t = 2
concentration in arbitrary units at t = 3
concentration in arbitrary units at t = 4
concentration in arbitrary units at t = 5
concentration in arbitrary units at t = 6
The simulation algorithm uses a finite differencing scheme with backwards differencing. The model used for simulation is a reaction diffusion-advection equation in which the advection term is variable in space but diffusion and reactions are constant in space see convection-diffusion equation for an example.
The parameters used in the general partial differntial equation in the link above are
D = 0.01 per squared spatial unit
R = 0.5 per unit time
v (advection is variable in space): in the upper left quadrant of the square domain v = (0.2, 0); in the upper right quadrant v = (0, -0.2); in the lower right quadrant v = (-0.2, 0); in the lower right quadrant v = (0, 0.2). Obviously v is discontinous at the quadrant boundaries, which causes some interesting model behaviour that is limited by considering only the first six time steps such that the bulk of the concentration in each quadrant does not cross a quadrant boundary.
The intial condition at time = 0 is a concentration of one unit per arbitrary unit of volume in the central cell of each quadrant.
External boundary conditions are zero-gradient (reflecting).
The data are formatted such that they can easily be converted to a raster stack using ICvectorfields::RastStackData(SimData).
Functions that facilitate calculation of statistics at the sub-grid level.
These may be useful for drivers of movement speed or direction if used in
tandem with DispField
, DispFieldST
, or
DispFieldSTall
.
SubgridMoransI(inputrast1, factv1, facth1, rad1 = 1) SubgridStats(inputrast1, factv1, facth1, statistic = "var")
SubgridMoransI(inputrast1, factv1, facth1, rad1 = 1) SubgridStats(inputrast1, factv1, facth1, statistic = "var")
inputrast1 |
a raster as produced by terra::rast |
factv1 |
an odd integer for the vertical dimension of sub-grids |
facth1 |
an odd integer for the horizontal dimension of sub-grids |
rad1 |
an integer indicating the neighbourhood radius for Moran's I statistic calculations in rows/columns. Any cell within a distance of rad1 cells of the focal cell is considered to be in its neighbourhood. |
statistic |
desired output statistic: It should be one of "mean", "var", or "sum". Default setting is var. |
Note that when using radius to define the neighbourhood in Moran's I calculations, a radius of one corresponds to the rook's neibhourhood. Values that are NA or Inf are not included in calculations of the Moran's I statistic nor in any of the other statistics that can be computed.
A data frame is returned with the following column names: rowcent, colcent, frowmin, frowmax, fcolmin, fcolmax, and a column for the output statistic.
DispStats
and DispMoransI
for functions
that compute statistics at presumed source or sink locations in each region
of interest.
(TestMat <- matrix(c(1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1), nrow = 5)) TestRast <- terra::rast(TestMat) terra::plot(TestRast) SubgridMoransI(TestRast, factv1 = 5, facth1 = 5, rad1 = 1) # using rad1 = 1 is equivalent to using the rooks neighbourhood # and so the output should be -1. (TestMat <- matrix(c(1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1), nrow = 5)) TestRast <- terra::rast(TestMat) terra::plot(TestRast) SubgridStats(TestRast, factv1 = 5, facth1 = 5, statistic = "mean") SubgridStats(TestRast, factv1 = 5, facth1 = 5, statistic = "var") SubgridStats(TestRast, factv1 = 5, facth1 = 5, statistic = "sum")
(TestMat <- matrix(c(1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1), nrow = 5)) TestRast <- terra::rast(TestMat) terra::plot(TestRast) SubgridMoransI(TestRast, factv1 = 5, facth1 = 5, rad1 = 1) # using rad1 = 1 is equivalent to using the rooks neighbourhood # and so the output should be -1. (TestMat <- matrix(c(1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1), nrow = 5)) TestRast <- terra::rast(TestMat) terra::plot(TestRast) SubgridStats(TestRast, factv1 = 5, facth1 = 5, statistic = "mean") SubgridStats(TestRast, factv1 = 5, facth1 = 5, statistic = "var") SubgridStats(TestRast, factv1 = 5, facth1 = 5, statistic = "sum")
This function efficiently computes two dimensional cross-covariance of two equal dimensioned matrices of real numbers using efficient discrete fast Fourier trasforms.
Xcov2D(mat1, mat2)
Xcov2D(mat1, mat2)
mat1 |
a real valued matrix |
mat2 |
a real valued matrix of equal dimension to mat1 |
The algorithm first pads each matrix with zeros so that the outer edges of the matrices do not interact with one another due to the circular nature of the discrete fast Fourier transform. Cross-covariance calculations require computation of the complex conjugate of one of the two imput matrices. Assuming all of it's elements are real, computing the complex conjugate is equivalent to flipping the matrix in the horizontal and vertical directions. Then to compute cross-covariance, the first matrix is convolved with the flipped second matrix as described in the convolution theorem.
This function is called by the main functions that compute displacement fields and vector fields and is included here primarily for demonstration purposes. Specifically, the method for computing the magnitude and direction of shifts is demonstrated in the examples.
The shift that produces the maximum cross-covariance between the two input matrices can be obtained by finding the row and column indices associated with the maximum cross-covariance. The shift in each direction is obtained by subracting one plus the half the dimension of the output matrix (the same for rows and columns) from the row and column values that are associated with the maximum cross-covariance as demonstrated in the examples below. Note that shifts to the right and up are denoted with positive numbers and shifts to the left and down are denoted by negative numbers. This is contrary to some conventions but efficient for producing vector fields. For more details on cross-covariance see cross-correlation.
a real valued matrix showing cross-covariance in each direction
matrix(c(1:6, rep(0, 3)), nrow = 3); matrix(c(rep(0, 3), 1:6), nrow = 3) dim(Xcov2D(matrix(c(1:6, rep(0, 3)), nrow = 3), matrix(c(rep(0, 3), 1:6), nrow = 3))) ICvectorfields::GetRowCol( which.max(Xcov2D(matrix(c(1:6, rep(0, 3)), nrow = 3), matrix(c(rep(0, 3), 1:6), nrow = 3))), dim1 = dim(Xcov2D(matrix(c(1:6, rep(0, 3)), nrow = 3), matrix(c(rep(0, 3), 1:6), nrow = 3)))[1], dim2 = dim(Xcov2D(matrix(c(1:6, rep(0, 3)), nrow = 3), matrix(c(rep(0, 3), 1:6), nrow = 3)))[2] ) # This implies that the shift is 6 - (10/2 + 1) in the vertical # direction and 7 - (10/2 + 1) in the horizonatal direction.
matrix(c(1:6, rep(0, 3)), nrow = 3); matrix(c(rep(0, 3), 1:6), nrow = 3) dim(Xcov2D(matrix(c(1:6, rep(0, 3)), nrow = 3), matrix(c(rep(0, 3), 1:6), nrow = 3))) ICvectorfields::GetRowCol( which.max(Xcov2D(matrix(c(1:6, rep(0, 3)), nrow = 3), matrix(c(rep(0, 3), 1:6), nrow = 3))), dim1 = dim(Xcov2D(matrix(c(1:6, rep(0, 3)), nrow = 3), matrix(c(rep(0, 3), 1:6), nrow = 3)))[1], dim2 = dim(Xcov2D(matrix(c(1:6, rep(0, 3)), nrow = 3), matrix(c(rep(0, 3), 1:6), nrow = 3)))[2] ) # This implies that the shift is 6 - (10/2 + 1) in the vertical # direction and 7 - (10/2 + 1) in the horizonatal direction.