kernelsmooth: function to apply a kernel to each cell in a grid or matrix...

In ethanplunkett/gridkernel: Package to apply kernels to raster data with NA and edge correction

function to apply a kernel to each cell in a grid or matrix with edge and NA correction

Description

These functions apply a kernel to a every cell in a grid or matrix. The returned object will be of the same type and dimension of the input but every cell will contain the sum of the product of the kernel and the corresponding cells of the input. If the kernel.fun argument is used than depending on the value of max.r and kernel.dim the grid method may perform upscaling before applying the kernel, this makes the algorythm much faster for large bandwidth smoothing operations.

Usage

kernelsmooth(x, ...)

## S3 method for class 'grid'
kernelsmooth(x, kernel, ...)
## S3 method for class 'matrix'
kernelsmooth(x, kernel, use.old=FALSE, ...)
## S3 method for class 'grid'
kernelsmooth(x, kernelfun, max.r, kernel.dim = 21, max.p.na = 1,
sample.points, weights, no.match = 0, use.old = FALSE, na.value = NA, ...)

Arguments

x	An object of class "grid" or "matrix" containing the raster data to be smoothed.
kernel	A kernel. It should be a square, odd dimensioned, numeric matrix that sums to 1.Use either this argument or kernelfun.
kernelfun	A function that generates a square, odd dimensioned, matrix that sums to 1. It must have max.r and cellsize arguments. Other arguments to kernelfun can be passed via .... Use either this argument or kernel, if you use both kernel will be ignored. See makegaussiankernel for one such function.
max.r	The maximum radius that weights are calculated for (in map units) - will be passed to kernelfun and used along with kernel.dim to determine how much upscaling is appropriate.
kernel.dim	The cell dimension of the kernel in the upscaled grid. It must be an odd number. This controls the amount of upscaling, which is adjusted such that a max.r radius kernel has this cell dimension. Smaller values result in more upscaling and faster processing at the expense of precision. The default value 21 is a good trade off between speed and precision for guassian kernels but may not be appropriate for other kernels. See
max.p.na	If the proportion of NA values (from the original resolution) in an upscaled cell exceeds max.p.na than it will receive an NA value during upscaling; otherwise the upscaled cell is given the mean value of all non-NA cells that it covers.
sample.points	(optional) a dataframe or matrix with columns "x" and "y" containing the map unit coordinates of points. If specified the character method will return the a vector of sample values (from the smoothed surface), and the grid method will return both the samples and a grid (in a list).
weights	(optional) a dataframe or matrix with columns "value" and "weight". If weights is provided than all the cells in the grid will be replaced with their associated weight; any cell whose value is not in the value column will be replaced with no.match, NA cells in x will be replaced with na.value. If weights does not contain columns named "value" and "weight" and contains only two columns it will be assumed that the first is "value" and the second is "weight".
no.match	If weights are supplied Non-NA values in x that don't match any value in original will be replaced with this value.
na.value	If weights are supplied this determines the value to use when cells in x are NA. Note if the weights argument contains a value for NA that will take precedence.
use.old	if FALSE then the newer Rcpp based implementation is used. Otherwise the old, pure R, version is used. This argument is here to facilitate testing; leave it at FALSE.
...	Arguments to be passed to other methods. If x is a character (path) these are passed on to kernelfun.

Details

If a gaussian kernel is used than the result is a smoothed version of the original grid or matrix. Other kernels could be used for other types of smoothing, or potentially be used to enhance edges or increase contrast.

This function adjusts for missing values in the matrix and at the edges.

The matrix version is equivalent to calling calckernel on each cell of x to calculate the corresponding value in the output.

If the kernel argument is used the grid version is the same but returns a grid object. If the kernelfun argument is used than the grid version may upscale x to a coarser resolution, make a kernel appropriate for the new cellsize, apply the kernel to the upscaled grid and then downscale (with interpolation) back to the original resolution.

Value

an object of the same class as x in which each cell is the result of applying a kernel to the correspononding focal cell in x, or if sample.points are provided a list with a grid and the sampled values.

Author(s)

Ethan Plunkett

See Also

calckernel, makegaussiankernel

Examples

# Matrix example
m <- matrix(1:110, 10, 11)
m[5, ] <- 80
k1 <- matrix(1, 3, 3)  # 3 x 3 kernel with more weight in center
k1[2,2] <- 3
k1 <- k1/sum(k1)
sm1 <- kernelsmooth(m, k1)

# An alternative to gaussiansmooth
k2 <- makegaussiankernel(sd = 1.5,  )
sm2 <- kernelsmooth(m, k2)



## Not run: 
# Visualize the results
# Note image rotates the matrix 90 CCW when plotting it.
cols <- heat.colors(n = 110)
image(m, col = cols)  # original
image(sm1, col = cols) # 3 x 3 with extra weight in center
image(sm2, col = cols)  # 1.5 cell

## End(Not run)

# Grid example
g  <- as.grid(m, xll = 0, yll = 0, cellsize = 10)
gs1  <- kernelsmooth(g, k1)
gs2 <- kernelsmooth(g, k2)


## Not run: 
# Visualize
plot(g)
plot(gs1)
plot(gs2)

## End(Not run)

# With upscaling
# Upscaling only comes into play with large kernel bandwidth and a
# kernel function. Here will make  a big example landscape and
#  use make.gaussian.kernel as the kernel function the result will be the same
# as if we had used gaussiansmooth()

dim  <- 1000  # cells
m <- matrix(runif(dim^2,  min = 0, max  = 100), dim, dim)
m[500:520, ] <- 0  # Low value horizontal band
m[ , 500:520] <- 100   # High value vertical band
m[1:100, 1:100] <- NA  # Patch of NA in upper left
g <- as.grid(m, xll = 0, yll = 0, cellsize = 30)

# Do smoothing  with  upscaling
sd <- 700  # meters not cells
a <- system.time({
  gs <- kernelsmooth(g, kernelfun = makegaussiankernel,  max.r = 3 * sd, sd = sd  )
})


## Not run: 
# Compare to results without upscaling
#  SLOW!!
b <- system.time({
  gs2 <- kernelsmooth(g, kernelfun = makegaussiankernel, kernel.dim = Inf,
                      max.r = 3 * sd, sd = sd  )
})
(a/b)[3]   #  0.007
cor(as.numeric(gs$m), as.numeric(gs2$m), use = "complete")  # 0.99969
mean( as.numeric(gs$m  - gs2$m), na.rm = TRUE)  # 0.0013
# Upscaling took 0.007 of the time and is almost  perfectly correlated to the
# non-upscaled version.
# However, upscaling did produce some artifacts in the interior of the NA region.

## End(Not run)

#  Pareto Kernel example
pk <- makeparetokernel(scale=10, shape=.8, max.r=250, cellsize=30)

# Plot a cross  section of the kernel
plot(k[ , ceiling(nrow(k)/2)])

## Not run: 
# By passing the kernel directly:
ps <- kernelsmooth(g, pk)

# With a kernel function
# Note, pareto  kernels are pointy, you might want to increase kernel.dim
ps2 <- kernelsmooth(g, kernelfun = makeparetokernel, max.r = 250, scale = 10, shape = 0.8 )

## End(Not run)

ethanplunkett/gridkernel documentation built on March 2, 2024, 9:06 p.m.