Extract raster values to points with bilinear interpolation

Posted on by

A student recently asked me how exactly the R terra::extract() function worked when using method="bilinear" to get raster values for points. The help file rightly says that ‘With “bilinear” the returned values are interpolated from the values of the four nearest raster cells‘, but this wasn’t immediately clear without a visual example. So, I put together some code to illustrate it:

library(terra)

# create a simple raster map:
r <- rast(nrows = 3, ncols = 3, xmin = 0, xmax = 3, ymin = 0, ymax = 3)
values(r) <- 1:ncell(r)
plot(r)

# get the pixel centroids:
centr <- as.points(r)
plot(centr, add = TRUE)

# show pixel values on the map:
plot(r, main = "Pixel values")
text(centr, values(r))

With terra::extract(), when a point falls exactly at the centre of a pixel, “bilinear” gives the same result as the “simple” method, i.e. the value of the pixel where the point exactly falls:

# get the middle pixel centroid:
pnt <- centr[5, ]
plot(pnt, col = "blue", pch = 4, cex = 2, add = TRUE)

extract(r, pnt, ID = FALSE, method = "simple")
extract(r, pnt, ID = FALSE, method = "bilinear")
# both 5

Now here’s a function to check what values are extracted for points at other parts of this raster map:

checkPoint <- function(dx = 0, dy = 0, col = "black") {
  pnt_shift <- shift(pnt, dx = dx, dy = dy)
  plot(pnt_shift, col = col, add = TRUE)
  val <- extract(r, pnt_shift, ID = FALSE, method = "bilinear")
  text(pnt_shift, val, pos = 3, col = col, cex = 0.8, halo = TRUE)
}

Let’s use this function to see what values we get from this raster (with method="bilinear") for different points:

plot(r)
text(centr, values(r))

checkPoint(dx = 0.2, col = "blue")
checkPoint(dx = 0.5, col = "purple")
checkPoint(dx = 0.8, dy = 0.2, col = "steelblue")
checkPoint(dx = 0.5, dy = -0.5, col = "magenta")
checkPoint(dx = -0.3, dy = 0.3, col = "darkgreen")
checkPoint(dx = -1.3, dy = 0, col = "darkturquoise")
checkPoint(dx = -1.3, dy = 0.3, col = "red3")
checkPoint(dx = -1.3, dy = 1.2, col = "cadetblue")

So that’s how this works! Pretty cool, isn’t it? Consider using “bilinear” rather than the default “simple” method when extracting raster values to points with ‘terra‘. And thanks to Robert J. Hijmans for such a great package!

Comment