Interpolation in R

First, let’s load the data from the website. The data are stored as SpatialPointsDataFrame and SpatialPointsDataFrame objects. Most of the functions used in this exercise work off of these classes. The one exception is the direchlet function which requires a conversion to a ppp object.

You’ll note the line P@bbox <- W@bbox which forces the rectangular extent of the Texas map onto the point data object. This is an important step if the interpolation of the points are to cover the entire extent of Texas. Had this step been omitted, most of the interpolated layers would have been limited to the smallest rectangular extent enclosing the point object.

Thiessen polygons

The Thiessen polygons (or proximity interpolation) can be created using spatstat’s dirichlet function.

Many packages share the same function names. This can be a problem when these packages are loaded in a same R session. For example, the intersect function is available in the base, spatstat and raster packages–all of which are loaded in this current session. To ensure that the proper function is selected, it’s a good idea to preface the function name with the package name as in raster::intersect().

This tip will be used in the next chunk of code when calling the idw function which is available in both spatstat and gstat.

Note that the dirichlet function (like most functions in the spatsat package) require that the point object be in a ppp format hence the inline as.ppp(P) syntax.


The IDW output is a raster. This requires that we first create an empty raster grid, then interpolate the precipitation values to each unsampled grid cell. An IDW power value of 2 (idp=2.0) will be used.


Generate Kriged surface

Next, use the variogram model to generate a kriged interpolated surface. The krige function allows us to include the trend model thus saving us from having to de-trend the data, krige the residuals, then combine the two rasters. Instead, all we need to do is pass krige the trend formula f.1.

Generate the variance and confidence interval maps

The dat.krg object stores not just the interpolated values, but the variance values as well. These can be passed to the raster object for mapping as follows:

Are more readily interpretable map is the 95% confidence interval map which can be generated from the variance object as follows (the map values should be interpreted as the number of inches above and below the estimated rainfall amount).