generate.nonstationary.grf.object: Create object used for approximate computations with a...

In mathiasisaksen/GRFics: Efficient generation of noise functions from approximate Gaussian random fields

Description

Creates a list containing the elements necessary for doing essential computations with the non-stationary GRF approximation.

Usage

generate.nonstationary.grf.object(x.lower = -1, x.upper = 1,
  y.lower = -1, y.upper = 1, resolution.x = 100,
  resolution.y = 100, range.function = 1, scale.function = 1,
  vector.field = c(0, 0), initial.seed = 0, params = NULL,
  num.functions = NULL, function.names = NULL)

Arguments

x.lower, x.upper, y.lower, y.upper	Specifies the boundaries of the rectangular region that the GRF is defined on.
resolution.x, resolution.y	The number of grid cells used for the regular grid, in the x- and y-direction, respectively.
range.function	A strictly positive function that specifies the range of the GRF in each location. The function must be vectorized: It takes in an n x 2 matrix of locations and returns a vector of length n, containing the value of the range in each location.
scale.function	A strictly positive function that is used for controlling the scale (standard deviation) of the GRF in each location. The function must be vectorized: It takes in an n x 2 matrix of locations and returns a vector of length n, containing the value of the scale in each location.
vector.field	A vector field that controls nature of the non-stationarity in each location. For a given location, the length and the direction of the corresponding vector plays the same role as 'strength.parameter' and 'direction.parameter' in generate.grf.object. More precisely: Let X be a location. If 'r = range.function(X)', and 's' and 'theta' are the length and direction of 'vector.field(X)', then the longest range is going to be approximately 'r*(1+s)', along the direction of 'theta'. The shortest range is approximately 'r', along the direction perpendicular to 'theta'.
initial.seed	The initial seed value used when generating functions.
params	Optional, a list containing all of the above parameters.
num.functions	The number of functions to initialize the GRF object with.
function.names	The names of the functions to initialize the GRF object with.

Details

If the GRF object is to be initialized with functions, either num.functions or function.names must be specified.

Value

An object (technically, a list) containing:

model.components	A list containing the components defining the GRF: The precision matrix Q, and the Cholesky decomposition of Q, chol.object. The latter is computed and added when a function is added to the GRF object.
grid	A (resolution.x*resolution.y) x 2 matrix containing the center coordinates of the regular grid that the approximation is defined on.
initial.seed	The initial seed value used when generating functions.
params	A list of the remaining input parameters listed above.

Note

For an explanation of how 'range.function', 'scale.function' and 'vector.field' affect the GRF, see the corresponding parameters 'range.parameter', 'scale.parameter', 'strength.parameter' and 'direction.parameter' in generate.grf.object. While the parameter interpretations for the stationary GRF only hold as an approximation in the non-stationary case, it still provides useful intuition.

Author(s)

Mathias Isaksen [e-mail]

Examples

# Example of non-stationary GRF on the square [-1, 1]^2
library(GRFics)

# A range function that increases linearly with distance.
# The range is 0.01 in the origin (0, 0), and 0.5 in the corners (±1, ±1).
range.function = function(X) {
  r = sqrt(X[, 1]^2+X[, 2]^2)
  return(0.01 + (0.5-0.01)*r/sqrt(2))
}

# A spiral-like vector field, where the length increases with distance from origin.
vector.field = function(X) {
  r = sqrt(X[, 1]^2+X[, 2]^2)
  theta = atan2(X[, 2], X[, 1]) + pi/4
  return(data.frame(vx = 3*r*cos(theta), vy = 3*r*sin(theta)))
}

library(ggplot2)
# Plots of range.function and vector.field
plot.grid = generate.grid.centers(-1, 1, -1, 1, 20, 20)

range.df = cbind(plot.grid, range = range.function((plot.grid)))
# Length of vector.field is scaled by 0.02 for plotting
vector.df = cbind(plot.grid, 0.02*vector.field(plot.grid))
vector.length = sqrt(vector.df$vx^2+vector.df$vy^2)
# range.function is shown using color, while vector.field is shown using arrows
ggplot()+
  geom_raster(data = range.df, aes(x = x, y = y, fill = range))+
  geom_segment(data = vector.df, aes(x = x-vx/2, xend = x+vx/2,
                                     y = y-vy/2, yend = y+vy/2),
               arrow = arrow(length = unit(0.2*vector.length, "npc")))+
  scale_fill_gradient(low = "blue", high = "white")+
  coord_fixed()

# Prepare GRF object
grf.object = generate.nonstationary.grf.object(-1, 1, -1, 1, 100, 100,
                                 range.function = range.function,
                                 vector.field = vector.field,
                                 function.names = "Original",
                                 initial.seed = 1000)

# Extract data frame containing the 100 x 100 grid coordinates and
# the value of the GRF in each grid cell.
original.df = get.function.df(grf.object)
# Interpolate GRF on finer grid for comparison
interp.locations = generate.grid.centers(-1, 1, -1, 1, 500, 500)
interp.values = evaluate.grf(interp.locations, grf.object, rescale.method = "none")
interp.df = data.frame(x = interp.locations[, 1],
                       y = interp.locations[, 2],
                       z = interp.values,
                       name = "Interpolation")

# Plot of original and interpolated function
ggplot()+
  geom_raster(data = rbind(original.df, interp.df),
              aes(x = x, y = y, fill = z))+
  coord_fixed()+
  facet_grid(cols = vars(name))+
  scale_fill_gradient(low = "black", high = "white")

# Create new vector field, where the direction is given by the non-stationary GRF
vector.grid = grf.object$grid # Use same 100 x 100 grid as the GRF is defined on
# Use uniform.transform = T to ensure that values are between 0 and 1 and all directions are equally respresented,
# and multiply by 2pi to get directions between 0 and 2pi (in radians)
vector.direction = 2*pi*evaluate.grf(vector.grid, grf.object, rescale.method = "uniform")
nonstat.vector.df = data.frame(
  x = vector.grid$x,
  y = vector.grid$y,
  vx = 0.02*cos(vector.direction),
  vy = 0.02*sin(vector.direction))
# Plot of random vector field (should be viewed in a large window)
ggplot()+
  geom_segment(data = nonstat.vector.df, aes(x = x-vx/2, xend = x+vx/2,
                                     y = y-vy/2, yend = y+vy/2),
               arrow = arrow(length = unit(0.003, "npc")))+
  coord_fixed()

mathiasisaksen/GRFics documentation built on May 20, 2021, 5:55 a.m.