| Title: | Illustrate Various Methods of Interpolation for Irregular Data |
|---|---|
| Description: | Examples of interpolating irregular data, to illustrate the mechanics of various methods and some easy tools to run them. |
| Authors: | Michael Sumner [aut, cre], Marcelino de la Cruz [ctb] (wrote maybetin example function on R-Sig-Geo), Tomas Remenyi [ctb] |
| Maintainer: | Michael Sumner <[email protected]> |
| License: | GPL-3 |
| Version: | 0.4.0.9000 |
| Built: | 2026-10-14 08:35:51 UTC |
| Source: | https://github.com/hypertidy/guerrilla |
Take the dimension, extent and coordinate system off a RasterLayer or a
terra SpatRaster and return them as a grid_spec(). Values are not
carried over; this is for describing a target to interpolate onto.
as_grid(x)as_grid(x)
x |
a |
It exists so that a grid you already have still works as a target, and so that the claim these classes are all the same few numbers can be checked rather than taken on trust.
A guerrilla_grid.
if (requireNamespace("raster", quietly = TRUE)) { as_grid(raster::raster(raster::extent(0, 8, 0, 6), ncols = 4, nrows = 3)) }if (requireNamespace("raster", quietly = TRUE)) { as_grid(raster::raster(raster::extent(0, 8, 0, 6), ncols = 4, nrows = 3)) }
One row per cell, with the cell centre coordinates and the value.
## S3 method for class 'guerrilla_grid' as.data.frame(x, row.names = NULL, optional = FALSE, ...)## S3 method for class 'guerrilla_grid' as.data.frame(x, row.names = NULL, optional = FALSE, ...)
x |
a |
row.names, optional
|
ignored, for consistency with the generic |
... |
ignored |
A data frame of x, y and value.
g <- grid_spec(dimension = c(3, 2), extent = c(0, 3, 0, 2)) g$values <- 1:6 head(as.data.frame(g))g <- grid_spec(dimension = c(3, 2), extent = c(0, 3, 0, 2)) g$values <- 1:6 head(as.data.frame(g))
The values as a matrix, one row per grid row, with the top row of the grid
first. This is the arrangement image() and raster both think in, and it
is where most orientation bugs come from, so it gets its own function.
## S3 method for class 'guerrilla_grid' as.matrix(x, ...)## S3 method for class 'guerrilla_grid' as.matrix(x, ...)
x |
a |
... |
ignored |
A matrix with nrow rows and ncol columns.
g <- grid_spec(dimension = c(3, 2), extent = c(0, 3, 0, 2)) g$values <- 1:6 as.matrix(g)g <- grid_spec(dimension = c(3, 2), extent = c(0, 3, 0, 2)) g$values <- 1:6 as.matrix(g)
The weights that say how much of each corner a point is made of. They sum to one, and they are all non-negative exactly when the point is inside the triangle – which makes them a point-in-triangle test and an interpolation rule at the same time.
bary_weights(triangle, xy)bary_weights(triangle, xy)
triangle |
three corners, as a 3 by 2 matrix |
xy |
points to place within it, a two column matrix |
This is the whole of what geometry::tsearch(bary = TRUE) computes in C,
written out so it can be read. For a triangle with corners
, , and a point
, solve for the weights by Cramer's rule on
Estimating a value is then one line: the weighted sum of the corner values.
A matrix with one row per point and three columns, one weight per
corner of triangle.
find_triangle(), which uses this to locate points, and
grid_barycentric(), which uses it to interpolate.
tri <- cbind(c(0, 1, 0), c(0, 0, 1)) ## the corners themselves bary_weights(tri, tri) ## the centroid is one third of each bary_weights(tri, cbind(1/3, 1/3)) ## outside the triangle, a weight goes negative bary_weights(tri, cbind(1, 1))tri <- cbind(c(0, 1, 0), c(0, 0, 1)) ## the corners themselves bary_weights(tri, tri) ## the centroid is one third of each bary_weights(tri, cbind(1/3, 1/3)) ## outside the triangle, a weight goes negative bary_weights(tri, cbind(1, 1))
A gridded elevation layer covering the Southern Ocean sector surveyed by the BROKE-West voyage, used as a backdrop for the interpolation examples.
A RasterLayer of 300 rows by 1200 columns at 1/15 degree
resolution, in longitude/latitude on WGS84. Values are integer metres, from
-5574 to 3612, with 17280 of 360000 cells missing.
Derived from ETOPO2, cropped to 15 to 95 degrees east and 75 to 55 degrees
south, aggregated by a factor of 2, smoothed with a Gaussian focal weight,
and rounded to whole metres. Values are elevation relative to sea level, so
the ocean is negative and the Antarctic continent is positive. Cells around
the edge are NA where the focal window ran off the grid.
ETOPO2, via data-raw/bathy.R. See
inst/examples/afternoon_exercise.Rmd for how it was built.
Triangulation has nothing to say about two values at the same place, and Qhull refuses outright. Real data does this constantly: a station occupied twice, a transect that crosses its own track, coordinates rounded to a printed precision.
collapse_duplicates( xy, value, fun = mean, quiet = FALSE, fun_label = deparse(substitute(fun)) )collapse_duplicates( xy, value, fun = mean, quiet = FALSE, fun_label = deparse(substitute(fun)) )
xy |
coordinates, two columns |
value |
one value per row of |
fun |
function used to combine the values at a repeated coordinate |
quiet |
do not report how many coordinates were collapsed |
fun_label |
name of |
So decide explicitly. collapse_duplicates() reduces repeated coordinates to
one, combining their values with fun, and says how many it collapsed.
A list of xy and value, with no repeated coordinate.
xy <- cbind(c(1, 2, 2, 3), c(1, 2, 2, 3)) collapse_duplicates(xy, c(10, 20, 30, 40))xy <- cbind(c(1, 2, 2, 3), c(1, 2, 2, 3)) collapse_duplicates(xy, c(10, 20, 30, 40))
Converters to the raster classes people actually have installed. Each one is a handful of lines, because there is nothing in a grid that these classes do not also hold: the same dimension, the same extent, the same values in the same order.
as_raster(grid) as_terra(grid) as_gdalraster(grid, filename = NULL)as_raster(grid) as_terra(grid) as_gdalraster(grid, filename = NULL)
grid |
a |
filename |
file to write; defaults to a path in GDAL's |
as_gdalraster() has to write a file, since that is how GDAL works; by
default it writes to GDAL's in-memory filesystem. It returns the open
GDALRaster object, not the file name, and the caller has to $close() it.
The file name is $getFilename() on the returned object.
An object of the corresponding class: a RasterLayer, a
SpatRaster, or an open GDALRaster.
xy <- cbind(runif(50), runif(50)) g <- grid_barycentric(xy, xy[, 1]) if (requireNamespace("terra", quietly = TRUE)) as_terra(g)xy <- cbind(runif(50), runif(50)) g <- grid_barycentric(xy, xy[, 1]) if (requireNamespace("terra", quietly = TRUE)) as_terra(g)
For every point, the row of triangles whose triangle contains it, or NA
for a point in none of them.
find_triangle(xy, triangles, points)find_triangle(xy, triangles, points)
xy |
vertex coordinates, a two column matrix |
triangles |
one triangle per row, as three indices into |
points |
points to locate, a two column matrix |
This is what the old tri_pip() did by building one sp polygon per
triangle and calling over(). GEOS does the same job with a spatial index it
builds for you: geos::geos_intersects_matrix() puts the triangles in an
STRtree, so each point is only tested against the few whose bounding boxes it
falls in. The exact containment test is then bary_weights(), in R.
An integer vector with one element per row of points.
xy <- cbind(c(0, 1, 1, 0), c(0, 0, 1, 1)) triangles <- rbind(c(1, 2, 3), c(1, 3, 4)) find_triangle(xy, triangles, cbind(c(0.9, 0.1, 5), c(0.5, 0.9, 5)))xy <- cbind(c(0, 1, 1, 0), c(0, 0, 1, 1)) triangles <- rbind(c(1, 2, 3), c(1, 3, 4)) find_triangle(xy, triangles, cbind(c(0.9, 0.1, 5), c(0.5, 0.9, 5)))
Triangulate the input coordinates, then estimate a value at every cell of a target grid from the barycentric coordinates of the cell centre within the triangle that contains it.
grid_barycentric( x, value = NULL, grid = NULL, duplicates = mean, engine = c("geometry", "R"), ... )grid_barycentric( x, value = NULL, grid = NULL, duplicates = mean, engine = c("geometry", "R"), ... )
x |
coordinates, or coordinates carrying their value as z; see
|
value |
one value per coordinate, or |
grid |
a |
duplicates |
function used to combine values at repeated coordinates,
or |
engine |
|
... |
ignored |
Within a triangle the surface is the plane through its three corner values,
so the result is continuous, passes exactly through the input values, and
invents nothing beyond them. This is what MATLAB calls
griddata(method = "linear") and what GDAL's gdal_grid calls linear.
Cells outside the convex hull of the points are in no triangle, and are left
NA. Points sharing a coordinate are combined first, because a triangulation
cannot hold two values in one place; see collapse_duplicates().
A guerrilla_grid with values.
engine = "geometry" locates the cells and computes their weights in one
pass of geometry::tsearch(), in C. It is the one to use.
engine = "R" does the same work in the open: find_triangle() to locate
each cell, bary_weights() to weight it. It is slower, and it is the point
of this package – the two agree to floating point, and the second one can be
read.
xy <- cbind(runif(100), runif(100)) grid_barycentric(xy, xy[, 1] + xy[, 2]) ## the value can be the z of the coordinates grid_barycentric(cbind(xy, xy[, 1])) ## the two engines agree g1 <- grid_barycentric(xy, xy[, 1], grid_spec(xy, dimension = c(20, 20))) g2 <- grid_barycentric(xy, xy[, 1], grid_spec(xy, dimension = c(20, 20)), engine = "R") max(abs(g1$values - g2$values), na.rm = TRUE)xy <- cbind(runif(100), runif(100)) grid_barycentric(xy, xy[, 1] + xy[, 2]) ## the value can be the z of the coordinates grid_barycentric(cbind(xy, xy[, 1])) ## the two engines agree g1 <- grid_barycentric(xy, xy[, 1], grid_spec(xy, dimension = c(20, 20))) g2 <- grid_barycentric(xy, xy[, 1], grid_spec(xy, dimension = c(20, 20)), engine = "R") max(abs(g1$values - g2$values), na.rm = TRUE)
Assign each point to the cell that contains it, and combine the values that land in the same cell.
grid_bin(x, value = NULL, grid = NULL, fun = mean, ...)grid_bin(x, value = NULL, grid = NULL, fun = mean, ...)
x |
coordinates, or coordinates carrying their value as z; see
|
value |
one value per coordinate, or |
grid |
a |
fun |
function used to combine the values in one cell |
... |
ignored |
This is the simplest thing that can be called interpolation, and it is not
really interpolation at all: nothing is estimated between the points, cells
with no point in them stay NA, and the answer depends entirely on how big
the cells are. It is here because it is the honest baseline, and because the
question it raises – what should happen when two points land in one cell –
is a question every other method answers silently.
fun is that answer, made explicit. mean averages them, length counts
them, function(x) x[1] keeps the first, and looking at grid_bin(fun = length) next to any other method is the quickest way to see where a surface
is supported by data and where it is supported by the algorithm.
Points on the outer edge of the extent are outside the last cell and are
dropped, which is one reason to pad a grid; see grid_spec().
The whole implementation is vaster::cell_from_xy() and tapply(). No
raster class is involved, because none is needed.
A guerrilla_grid with values.
xy <- cbind(runif(300), runif(300)) g <- grid_spec(xy, dimension = c(20, 20)) plot(grid_bin(xy, xy[, 1], g)) ## how many points went into each cell plot(grid_bin(xy, xy[, 1], g, fun = length))xy <- cbind(runif(300), runif(300)) g <- grid_spec(xy, dimension = c(20, 20)) plot(grid_bin(xy, xy[, 1], g)) ## how many points went into each cell plot(grid_bin(xy, xy[, 1], g, fun = length))
Triangulate the points, fit value ~ x + y by least squares to the three
corners of each triangle, and evaluate that fit at the cells inside it.
grid_facet_lm(x, value = NULL, grid = NULL, duplicates = mean, ...)grid_facet_lm(x, value = NULL, grid = NULL, duplicates = mean, ...)
x |
coordinates, or coordinates carrying their value as z; see
|
value |
one value per coordinate, or |
grid |
a |
duplicates |
function used to combine values at repeated coordinates,
or |
... |
ignored |
This is the long way round to grid_barycentric(), and the two agree to
floating point. Three points determine a plane, least squares through three
points is that plane, and barycentric interpolation evaluates that same
plane – so the expensive per-triangle lm() and the one vectorised
tsearch() call compute the same surface. Worth seeing once.
It replaces the Delaunay half of facets(), which did this through
spatstat and returned a point pattern rather than a grid.
The Voronoi half of facets() is not here, because it was not doing what it
looked like it was doing. A Voronoi tile contains exactly one point by
construction, so lm(value ~ x + y) on it fits an intercept and nothing
else, and predicts that one point's value across the whole tile. That is
nearest neighbour, computed the slowest way available. grid_voronoi() is
the same answer under its right name.
A guerrilla_grid with values.
grid_barycentric() for the fast form of this, grid_voronoi()
for what the Voronoi case really was.
xy <- cbind(runif(40), runif(40)) v <- xy[, 1] + xy[, 2] plot(grid_facet_lm(xy, v))xy <- cbind(runif(40), runif(40)) v <- xy[, 1] + xy[, 2] plot(grid_facet_lm(xy, v))
Fit a smooth of x and y with mgcv::gam() and predict it at every cell.
grid_gam( x, value = NULL, grid = NULL, formula = value ~ s(x, y), statistic = c("prediction", "se"), ... )grid_gam( x, value = NULL, grid = NULL, formula = value ~ s(x, y), statistic = c("prediction", "se"), ... )
x |
coordinates, or coordinates carrying their value as z; see
|
value |
one value per coordinate, or |
grid |
a |
formula |
a model formula in |
statistic |
|
... |
passed to |
The default formula is value ~ s(x, y), one isotropic two dimensional
smooth. That matters: mgcv's default basis for a two dimensional s() is a
thin plate regression spline, so this is grid_tps() with the number of
basis functions reduced and the smoothing parameter chosen by REML. The two
give visibly similar surfaces on the same data, and seeing that is worth more
than either surface on its own.
value ~ s(x) + s(y) is the wrong model for a surface and a useful thing to
try: additive in x and y means every column of the grid has the same shape as
every other, so it cannot represent a feature that sits in one place. Pass it
as formula and look.
A guerrilla_grid with values.
xy <- cbind(runif(200), runif(200)) v <- exp(-20 * ((xy[, 1] - 0.5)^2 + (xy[, 2] - 0.5)^2)) if (requireNamespace("mgcv", quietly = TRUE)) { op <- par(mfrow = c(1, 2)) plot(grid_gam(xy, v), main = "s(x, y)") ## additive in x and y cannot put a bump in one place plot(grid_gam(xy, v, formula = value ~ s(x) + s(y)), main = "s(x) + s(y)") par(op) }xy <- cbind(runif(200), runif(200)) v <- exp(-20 * ((xy[, 1] - 0.5)^2 + (xy[, 2] - 0.5)^2)) if (requireNamespace("mgcv", quietly = TRUE)) { op <- par(mfrow = c(1, 2)) plot(grid_gam(xy, v), main = "s(x, y)") ## additive in x and y cannot put a bump in one place plot(grid_gam(xy, v, formula = value ~ s(x) + s(y)), main = "s(x) + s(y)") par(op) }
Hand the points to GDAL and let gdal_grid do the interpolation, returning
the result as a grid like everything else here.
grid_gdal( x, value = NULL, grid = NULL, algorithm = "linear", nodata = NaN, ... )grid_gdal( x, value = NULL, grid = NULL, algorithm = "linear", nodata = NaN, ... )
x |
coordinates, or coordinates carrying their value as z; see
|
value |
one value per coordinate, or |
grid |
a |
algorithm |
a |
nodata |
value GDAL should write for empty cells, read back as |
... |
further arguments appended to the |
Four of the methods in this package have a gdal_grid twin, and running
them side by side is the best available check that any of them is right.
They are independent implementations in different languages by different
people, so agreement is evidence and disagreement is a bug in one of them:
"linear:radius=0.0" |
grid_barycentric() |
"nearest" |
grid_voronoi() |
"invdist:power=2.0" |
grid_idw() |
"average:radius1=,radius2=" |
grid_bin(), with a search radius
|
A guerrilla_grid with values.
algorithm = "linear" on its own does not stop at the convex hull. GDAL's
default is radius=-1, an infinite search, so a cell in no triangle takes
the value of the nearest point instead of being left out. That is
extrapolation, it is silent, and it is why "linear:radius=0.0" rather than
"linear" is the one that matches grid_barycentric().
Second, gdal_grid fills cells it could not estimate with 0 and does not
record that 0 as the band's no data value, so nothing downstream can tell it
apart from a real measurement – and 0 is a real temperature, a real
elevation and a real anomaly. This function therefore appends nodata=nan
unless the algorithm string already sets one. That both fills the empty
cells with NaN and tags the band, and gdalraster turns the tagged
value into NA on the way back, so an explicit nodata=-9999 arrives here
as NA just the same.
gdal_grid is reachable from R only through sf::gdal_utils().
gdalraster wraps warp, translate and rasterize but not
GDALGrid, and the unified gdal command line added in GDAL 3.11 has no
grid subcommand either, so there is nothing for it to wrap yet.
grid_barycentric(), grid_voronoi(), grid_idw()
xy <- cbind(runif(50), runif(50)) if (requireNamespace("sf", quietly = TRUE) && requireNamespace("gdalraster", quietly = TRUE)) { g <- grid_spec(xy, dimension = c(40, 40)) gdal <- grid_gdal(xy, xy[, 1], g, algorithm = "linear:radius=0.0") ours <- grid_barycentric(xy, xy[, 1], g) max(abs(gdal$values - ours$values), na.rm = TRUE) }xy <- cbind(runif(50), runif(50)) if (requireNamespace("sf", quietly = TRUE) && requireNamespace("gdalraster", quietly = TRUE)) { g <- grid_spec(xy, dimension = c(40, 40)) gdal <- grid_gdal(xy, xy[, 1], g, algorithm = "linear:radius=0.0") ours <- grid_barycentric(xy, xy[, 1], g) max(abs(gdal$values - ours$values), na.rm = TRUE) }
Estimate each cell as a weighted average of the input values, weighted by
one over distance to the power idp.
grid_idw(x, value = NULL, grid = NULL, idp = 2, nmax = Inf, ...)grid_idw(x, value = NULL, grid = NULL, idp = 2, nmax = Inf, ...)
x |
coordinates, or coordinates carrying their value as z; see
|
value |
one value per coordinate, or |
grid |
a |
idp |
the inverse distance power |
nmax |
use at most this many of the nearest points per cell |
... |
passed to |
There is no model here, only a rule, and that is the point of having it next
to the others. Raising idp makes the nearest point dominate, so the
surface tends to nearest neighbour (grid_voronoi()); lowering it flattens
towards the overall mean. Neither end is more correct than the other, and
nothing in the data tells you where to sit between them.
The consequence is that this method cannot report a standard error, because
it never claimed to be estimating anything. grid_kriging() is the same
idea taken seriously: weights chosen from a fitted model of how the values
vary with distance, which then also says how uncertain each cell is.
A guerrilla_grid with values.
xy <- cbind(runif(50), runif(50)) if (requireNamespace("gstat", quietly = TRUE)) { op <- par(mfrow = c(1, 2)) plot(grid_idw(xy, xy[, 1], idp = 0.5), main = "idp = 0.5") plot(grid_idw(xy, xy[, 1], idp = 8), main = "idp = 8") par(op) }xy <- cbind(runif(50), runif(50)) if (requireNamespace("gstat", quietly = TRUE)) { op <- par(mfrow = c(1, 2)) plot(grid_idw(xy, xy[, 1], idp = 0.5), main = "idp = 0.5") plot(grid_idw(xy, xy[, 1], idp = 8), main = "idp = 8") par(op) }
Fit a variogram to the data, then use it to weight the input values for every cell.
grid_kriging( x, value = NULL, grid = NULL, model = NULL, statistic = c("prediction", "se"), ... )grid_kriging( x, value = NULL, grid = NULL, model = NULL, statistic = c("prediction", "se"), ... )
x |
coordinates, or coordinates carrying their value as z; see
|
value |
one value per coordinate, or |
grid |
a |
model |
a |
statistic |
|
... |
passed to |
The variogram is a picture of how much two values differ, as a function of
how far apart they are. Kriging fits a curve to that picture and uses it to
choose the weights, so unlike grid_idw() the weighting is estimated from
the data rather than chosen.
That also makes this the one method here that can refuse. If the values have no spatial structure at these distances there is no curve to fit, the fit comes back with a negative range, and this stops rather than kriging from nonsense. Every other method in the package will happily interpolate pure noise and hand you a picture of it.
Because there is a fitted model, there is also a variance: statistic = "se"
returns the square root of the kriging variance, which is small near the data
and rises to the sill away from it. That surface is usually more informative
than the prediction, and it is the reason to reach for this method at all.
model is passed to gstat::fit.variogram() as the starting point. The
default is a spherical model with the sill at the overall variance and the
range at a third of the distances the variogram covers. That is a starting
guess, not an analysis: a real one plots gstat::variogram() first and
chooses. If the fit does not converge, that is what it is telling you.
A guerrilla_grid with values.
xy <- cbind(runif(80), runif(80)) v <- xy[, 1] * 3 + rnorm(80, sd = 0.1) if (requireNamespace("gstat", quietly = TRUE)) { plot(grid_kriging(xy, v)) plot(grid_kriging(xy, v, statistic = "se")) }xy <- cbind(runif(80), runif(80)) v <- xy[, 1] * 3 + rnorm(80, sd = 0.1) if (requireNamespace("gstat", quietly = TRUE)) { plot(grid_kriging(xy, v)) plot(grid_kriging(xy, v, statistic = "se")) }
Number of cells in a grid
grid_ncell(grid)grid_ncell(grid)
grid |
a |
A single integer, the product of the two dimensions.
grid_ncell(grid_spec(dimension = c(3, 2), extent = c(0, 6, 0, 4)))grid_ncell(grid_spec(dimension = c(3, 2), extent = c(0, 6, 0, 4)))
Cell size in each direction
grid_res(grid)grid_res(grid)
grid |
a |
Two numbers, the x and y size of one cell.
grid_res(grid_spec(dimension = c(4, 2), extent = c(0, 8, 0, 4)))grid_res(grid_spec(dimension = c(4, 2), extent = c(0, 8, 0, 4)))
Bin the points to a grid, then smooth that grid with a Gaussian kernel,
using fields::image.smooth().
grid_smooth(x, value = NULL, grid = NULL, theta = NULL, ...)grid_smooth(x, value = NULL, grid = NULL, theta = NULL, ...)
x |
coordinates, or coordinates carrying their value as z; see
|
value |
one value per coordinate, or |
grid |
a |
theta |
kernel bandwidth, in the units of the coordinates |
... |
passed to |
This is the two step method: rasterize first (grid_bin()), estimate second.
Everything else here works from the points directly, so this one is worth
having as the contrast – the binning throws away where inside its cell each
point was, and no amount of smoothing afterwards gets that back.
theta is the kernel bandwidth in coordinate units, and it does all the
work. There is no standard error, because there is no model.
fields::image.smooth() wants a square grid, so this uses one of
max(dimension) cells a side and returns it; the result covers the same
extent but need not have the dimension you asked for.
A guerrilla_grid with values.
xy <- cbind(runif(200), runif(200)) if (requireNamespace("fields", quietly = TRUE)) { plot(grid_smooth(xy, xy[, 1] + xy[, 2], theta = 0.1)) }xy <- cbind(runif(200), runif(200)) if (requireNamespace("fields", quietly = TRUE)) { plot(grid_smooth(xy, xy[, 1] + xy[, 2], theta = 0.1)) }
Describe a regular grid: how many cells across and down, what region of the plane it covers, and optionally what coordinate system that plane is in.
grid_spec( x = NULL, dimension = c(60L, 50L), extent = NULL, crs = NULL, pad = 0 )grid_spec( x = NULL, dimension = c(60L, 50L), extent = NULL, crs = NULL, pad = 0 )
x |
coordinates to cover: anything with two or more columns, such as a
matrix or data frame. Ignored when |
dimension |
number of columns and rows, in that order |
extent |
outer edges |
crs |
coordinate reference system, or |
pad |
fraction of each axis range to add on every side, so the outermost
input coordinates do not sit exactly on the boundary. |
The value returned is a list. That is deliberate. Every raster in every package is some version of these few numbers, and a reader who can print the whole thing has actually been told what a raster is:
dimension, the number of columns and rows, in that order;
extent, the outer edges as xmin, xmax, ymin, ymax;
crs, a coordinate reference system, or NULL for none;
values, one number per cell, or NULL for a grid that is only a target.
dimension and extent are exactly the first two arguments of the
vaster functions, so grid$dimension and grid$extent can be handed
straight to vaster::xy_from_cell() and friends. Cells are numbered left to
right along the top row first, the same order raster and GDAL use.
Nothing here assumes the plane is the Earth. crs stays NULL unless you
say otherwise.
A guerrilla_grid, which is a list of dimension, extent, crs
and values.
xy <- cbind(c(0, 4, 10), c(0, 2, 5)) grid_spec(xy) ## it really is just a list str(unclass(grid_spec(xy))) ## give the points some room grid_spec(xy, pad = 0.05)$extent ## or skip the points entirely grid_spec(dimension = c(4, 3), extent = c(0, 8, 0, 6))xy <- cbind(c(0, 4, 10), c(0, 2, 5)) grid_spec(xy) ## it really is just a list str(unclass(grid_spec(xy))) ## give the points some room grid_spec(xy, pad = 0.05)$extent ## or skip the points entirely grid_spec(dimension = c(4, 3), extent = c(0, 8, 0, 6))
Fit a thin plate spline to the points with fields::Tps() and evaluate it
at every cell.
grid_tps( x, value = NULL, grid = NULL, statistic = c("prediction", "se"), lon.lat = FALSE, ... )grid_tps( x, value = NULL, grid = NULL, statistic = c("prediction", "se"), lon.lat = FALSE, ... )
x |
coordinates, or coordinates carrying their value as z; see
|
value |
one value per coordinate, or |
grid |
a |
statistic |
|
lon.lat |
the coordinates are longitude and latitude, so distances
should be great circle ones; passed to |
... |
passed to |
A thin plate spline is the surface that passes near the data while bending
as little as possible, with a smoothing parameter deciding the trade between
those two. It is the engine behind GDAL's -tps warping, so this is also
what happens when you georeference an image from ground control points.
Unlike grid_barycentric() it fills the whole grid, including cells far
outside the data, where the answer is extrapolation and should be read as
such. statistic = "se" is how it tells you that: the standard error grows
with distance from the data, so the two surfaces together say both what the
method thinks and where it is guessing.
lon.lat is the one place in this package where what the coordinates mean
changes the arithmetic. A spline bends in the plane of its coordinates, and
a degree of longitude is not a degree of latitude anywhere but the equator.
It is not guessed from the crs, because the guess would be wrong exactly
when it mattered; say what you have.
A guerrilla_grid with values.
xy <- cbind(runif(60), runif(60)) v <- sin(xy[, 1] * 6) + xy[, 2] if (requireNamespace("fields", quietly = TRUE)) { plot(grid_tps(xy, v)) ## and where it is guessing plot(grid_tps(xy, v, statistic = "se")) }xy <- cbind(runif(60), runif(60)) v <- sin(xy[, 1] * 6) + xy[, 2] if (requireNamespace("fields", quietly = TRUE)) { plot(grid_tps(xy, v)) ## and where it is guessing plot(grid_tps(xy, v, statistic = "se")) }
Give every cell the value of the closest input point.
grid_voronoi(x, value = NULL, grid = NULL, duplicates = mean, ...)grid_voronoi(x, value = NULL, grid = NULL, duplicates = mean, ...)
x |
coordinates, or coordinates carrying their value as z; see
|
value |
one value per coordinate, or |
grid |
a |
duplicates |
function used to combine values at repeated coordinates,
or |
... |
ignored |
The tessellation is the picture of that rule. A Voronoi tile is exactly the region closer to its own point than to any other, so assigning each cell to the tile it falls in is nearest neighbour, drawn rather than computed.
Unlike grid_barycentric(), this fills the whole grid. Voronoi tiles cover
the plane, so there is no convex hull to fall outside of – which is
convenient and is also the method's main way of lying to you, since a cell
far beyond the data gets a confident answer from one distant point.
A guerrilla_grid with values.
xy <- cbind(runif(30), runif(30)) plot(grid_voronoi(xy, xy[, 1])) points(xy, pch = 16, cex = 0.5)xy <- cbind(runif(30), runif(30)) plot(grid_voronoi(xy, xy[, 1])) points(xy, pch = 16, cex = 0.5)
The centre of every cell of grid, in cell order: left to right along the
top row, then the next row down.
grid_xy(grid)grid_xy(grid)
grid |
a |
This is one call to vaster::xy_from_cell() with the grid's own first two
elements. It is the only thing the interpolation functions need to know
about the grid, which is why the grid can be a list.
A two column matrix with one row per cell.
grid_xy(grid_spec(dimension = c(3, 2), extent = c(0, 6, 0, 4)))grid_xy(grid_spec(dimension = c(3, 2), extent = c(0, 6, 0, 4)))
Is this a guerrilla grid?
is_grid(x)is_grid(x)
x |
any object |
TRUE or FALSE
is_grid(grid_spec(cbind(1:3, 1:3))) is_grid(1:10)is_grid(grid_spec(cbind(1:3, 1:3))) is_grid(1:10)
Create a grid by interpolating across triangles.
mesh_raster(x, grid = NULL, n = 128, ...) ## S3 method for class 'mesh3d' mesh_raster(x, grid = NULL, n = 128, ...) ## S3 method for class 'matrix' mesh_raster(x, grid = NULL, n = 128, ...) ## S3 method for class 'data.frame' mesh_raster(x, grid = NULL, n = 128, ...)mesh_raster(x, grid = NULL, n = 128, ...) ## S3 method for class 'mesh3d' mesh_raster(x, grid = NULL, n = 128, ...) ## S3 method for class 'matrix' mesh_raster(x, grid = NULL, n = 128, ...) ## S3 method for class 'data.frame' mesh_raster(x, grid = NULL, n = 128, ...)
x |
a matrix or data frame of three columns, or a |
grid |
a |
n |
cells across, when |
... |
passed to methods, and on to |
For x-y-z input this is grid_barycentric() with the arguments arranged
differently:
triangulate, then estimate each cell from the triangle containing it. The
difference is that mesh_raster() also takes a mesh3d, where the triangles
already exist and no triangulation is needed.
That case is worth having because it runs the argument backwards. A grid is one particular mesh – a regular one, with the values at the vertices – so turning an arbitrary mesh into a grid is resampling one mesh onto another, and barycentric interpolation is how you do it either way.
A guerrilla_grid with values.
## interpolate from raw points xyz <- quakes[c("long", "lat", "depth")] xyz$depth <- -xyz$depth gx <- mesh_raster(xyz) rat <- 1/cos(mean(xyz[["lat"]]) * pi/180) plot(gx, asp = rat, col = hcl.colors(12, "YlOrRd")) points(xyz, pch = "+", cex = 0.3) ## add some dummy points (we aren't modelling the world) xex <- cbind(expand.grid(long = range(xyz$long), lat = range(xyz$lat)), depth = 0) g2 <- mesh_raster(rbind(xex, xyz)) plot(g2, asp = rat) points(xyz, pch = "+", cex = 0.3) if (requireNamespace("maps", quietly = TRUE)) { maps::map(add = TRUE) } ## a mesh of triangles, rather than raw points if (requireNamespace("Rvcg", quietly = TRUE)) { data("humface", package = "Rvcg") grid <- mesh_raster(humface, n = 256) plot(grid, col = grey.colors(24)) } ## for the other direction, a grid drawn as a 3D mesh, see ## vignette("triangulation")## interpolate from raw points xyz <- quakes[c("long", "lat", "depth")] xyz$depth <- -xyz$depth gx <- mesh_raster(xyz) rat <- 1/cos(mean(xyz[["lat"]]) * pi/180) plot(gx, asp = rat, col = hcl.colors(12, "YlOrRd")) points(xyz, pch = "+", cex = 0.3) ## add some dummy points (we aren't modelling the world) xex <- cbind(expand.grid(long = range(xyz$long), lat = range(xyz$lat)), depth = 0) g2 <- mesh_raster(rbind(xex, xyz)) plot(g2, asp = rat) points(xyz, pch = "+", cex = 0.3) if (requireNamespace("maps", quietly = TRUE)) { maps::map(add = TRUE) } ## a mesh of triangles, rather than raw points if (requireNamespace("Rvcg", quietly = TRUE)) { data("humface", package = "Rvcg") grid <- mesh_raster(humface, n = 256) plot(grid, col = grey.colors(24)) } ## for the other direction, a grid drawn as a 3D mesh, see ## vignette("triangulation")
Draw the grid with graphics::image(). No raster class is involved, which
is the point: a grid is a dimension, an extent and a vector of values, and
that is enough to draw it.
## S3 method for class 'guerrilla_grid' plot(x, col = grDevices::hcl.colors(24, "YlGnBu"), asp = 1, ...)## S3 method for class 'guerrilla_grid' plot(x, col = grDevices::hcl.colors(24, "YlGnBu"), asp = 1, ...)
x |
a |
col |
colours |
asp |
aspect ratio; |
... |
passed to |
x, invisibly.
xy <- cbind(runif(50), runif(50)) g <- grid_barycentric(xy, xy[, 1] + xy[, 2]) plot(g)xy <- cbind(runif(50), runif(50)) g <- grid_barycentric(xy, xy[, 1] + xy[, 2]) plot(g)
Every interpolation function here takes the same thing: places, and a value at each place. This is what turns the various ways of expressing that into one matrix of coordinates and one vector of values.
xyz_input(x, value = NULL)xyz_input(x, value = NULL)
x |
coordinates: a matrix, a data frame, or anything wk handles |
value |
one value per coordinate, or |
A value can arrive two ways, and they are the same thing said differently:
as a separate value argument, or as the z of three dimensional
coordinates. Barycentric interpolation really is treating the value as a
height above the plane, so wk::xyz() input is not an abuse of notation
here, it is the notation. For kriging or a GAM it would be, which is why
value exists too.
Objects from wk and sf go through wk::wk_coords(), which also
carries their coordinate reference system across. A plain matrix or data
frame is read by position: first column x, second y, third z if present.
wk::as_xy() is deliberately not used for those, because it refuses to
guess at column names, and Lon/Lat is exactly what real data is called.
A list of xy (a two column matrix), value, and crs.
xyz_input(cbind(1:3, 4:6), c(10, 20, 30)) ## z is the value xyz_input(cbind(1:3, 4:6, c(10, 20, 30))) ## and the crs comes along xyz_input(wk::xyz(1:3, 4:6, c(10, 20, 30), crs = "EPSG:4326"))$crsxyz_input(cbind(1:3, 4:6), c(10, 20, 30)) ## z is the value xyz_input(cbind(1:3, 4:6, c(10, 20, 30))) ## and the crs comes along xyz_input(wk::xyz(1:3, 4:6, c(10, 20, 30), crs = "EPSG:4326"))$crs