Inverse distance weighted (IDW) interpolation using spatstat

So many packages so little time

You can do inverse-distance weighted interpolation in the spatstat package, but it’s a little bit fiddly. This post has code examples to show you how.
R
geospatial
tutorial
Author

David O’Sullivan

Published

October 22, 2021

Modified

November 21, 2024

One of the joys (ahem) of R spatial is moving data around between formats so you can use the best packages for particular jobs. Here’s an example using IDW interpolation in spatstat.

Libraries

Libraries are the usual suspects plus spatstat (duh) and maptools for some extra conversions. We also need terra for the data prep.

library(sf)
library(dplyr)
library(tmap)
library(spatstat)
library(terra)

Introduction

As is often the case there is useful functionality in a package that doesn’t play nice with the core R-spatial packages. spatstat is really great for lots of things, but does not support sf and even needs a bit of persuading to handle sp data. Its implementation of IDW interpolation is nice however, so it’s nice to know how to use it. Whether or not you should ever use IDW is another question altogether, but we can worry about that some other time.

Data

First we need a set of values to interpolate. I made a projected version of the R core dataset volcano which is a nice place to start.

maungawhau <-  rast("maungawhau.tif")

Some random control points and a study area

We can get a dataframe of random points on the surface using terra::spatSample. We’ll make this into a sf object as a starting point because that’s the most likely situation when you want to interpolate data (you will have an sf source).

pts <- maungawhau |>
  spatSample(500, xy = TRUE) |>
  st_as_sf(coords = c("x", "y")) |>
  st_set_crs(2193) |>
  st_jitter(5)

We also need a spatial extent for the interpolation, so let’s just make a convex hull of the points

spatial_extent <- pts |>
  st_union() |>
  st_convex_hull() |>
  st_sf()

And just to see where we are at

tm_shape(maungawhau) + 
  tm_raster(col.scale = tm_scale_continuous(values = "hcl.terrain2")) +
  tm_shape(pts) + 
  tm_dots() + 
  tm_shape(spatial_extent) + 
  tm_borders() + 
  tm_layout(legend.frame = FALSE)

Make the data into a spatstat point pattern

spatstat has its own format for point patterns, including coordinates, marks (the values) and a window or owin (the spatial extent). It’s best to make the window first and then we can make the whole thing all at once. spatstat prefers sp objects, so we go via ‘Spatial’ to get a spatstat::owin object. maptools provides the conversion to an owin.

W <- spatial_extent |>
  as.owin()

We also need the control point coordinates

xy <- pts |>
  st_coordinates() |>
  as_tibble()

Now we can make a spatstat::ppp point pattern

pp <- ppp(x = xy$X, xy$Y, marks = pts$maungawhau, window = W)
plot(pp)

Success!

A previous notebook showed an even quicker way to do this, but where the window will be formed from a bounding box (and where’s the fun in that?)

pts |>
  as("Spatial") |>
  as.ppp()

Interpolation

It’s easy from here. power is the inverse power applied to distances, and eps is the resolution in units of the coordinate system.

result <- idw(pp, power = 2, eps = 10)
plot(result)

This can be converted back to a terra raster for comparison with the original surface.

# stack the layers so we can 'facet' plot them
comparison_raster <- rast(
  list(maungawhau = maungawhau, 
       interpolation = rast(result) |> resample(maungawhau)
  ))

tm_shape(comparison_raster) +
  tm_raster(col.scale = tm_scale_continuous(values = "hcl.terrain2"),
            col.free = FALSE,
            col.legend = tm_legend(
              title = "Elevation", 
              position = tm_pos_out("right", "center"), 
              frame = FALSE))

Like I said, IDW is not necessarily a great interpolation method!

Geospatial Stuff