Contributors: Clark S. Rushing
This activity will introduce you to spatial interpolation in R.
R Skill Level: Intermediate - this activity assumes you have a working knowledge of
This activity incorporates the skills learned in the Points, Raster & Projections activities.
Download R script Last modified: 2019-09-20 18:26:28
This is still very much a work in progress
One basic principle of geography is that variables that are close in space typically have similar values. We can use this principle to make predictions and interpolate a smooth surface from data sampled at points. There are a few ways to do this. Generally, when people refer to ‘kriging’ they are referring to interpolation.
The following examples will use the raw number of species detected along breed bird survey (BBS) routes within Arizona. The data were summarized for brevity but the data represent the number of species detected along each route through time. We didn’t do any formal analysis that accounts for imperfect detection. We just wanted some real data to illustrate interpolation.
library(raster) library(sp) library(gstat) library(spatstat)
Now that we have some data let’s interpolate (predict values for locations were we have no data) species richness across Arizona.
Let’s make an empty raster so that we can predict species richness across AZ. We’ll use this several times so we’ll go ahead and make it here.
# Spatial grid over the extent of sppStateRoute # approximately 5000 random points # set.seed() makes sure the random points are the same each time # increases reproducability set.seed(12345) samplegrid <- raster::raster(sppStateRoute,res = c(0.025,0.025)) # define CRS raster::crs(samplegrid) <- raster::crs(sppStateRoute) <- sp::CRS("+init=epsg:4326")
We’ll use the
gstat package to interpolate species richness using a few different methods.
Inverse distance weighting
idw.model <- gstat(formula=sppStateRoute$SpeciesDetected~1, locations=sppStateRoute) idw.spp <- raster::interpolate(samplegrid, idw.model)
## [inverse distance weighted interpolation]
## class : RasterLayer ## dimensions : 221, 222, 49062 (nrow, ncol, ncell) ## resolution : 0.025, 0.025 (x, y) ## extent : -114.6683, -109.1183, 31.38823, 36.91323 (xmin, xmax, ymin, ymax) ## crs : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ## source : memory ## names : var1.pred ## values : 11.42869, 75.10614 (min, max)
Let’s use the
mask function to clip the raster to Arizona.
idw.spp.az <- raster::mask(idw.spp, States[States$NAME_1=="Arizona",]) par(bty="n",mar = c(0,0,0,0)) raster::plot(idw.spp.az, axes = FALSE, legend.args = list("Species\nRichness")) plot(States, add = TRUE)
# remove duplicate locations because otherwise Kriging returns NA sppStateRoute <- sppStateRoute[-sp::zerodist(sppStateRoute)[,1],] # create variogram vario <- variogram(sppStateRoute$SpeciesDetected~1, sppStateRoute)
vario.fit <- fit.variogram(vario, vgm(psill=250, # approx asymptote range = 150, model="Exp", nugget=50)) predict.model <- gstat(g=NULL, formula = sppStateRoute$SpeciesDetected~1, locations = sppStateRoute, model=vario.fit) # make an empty raster to predict on # r <- raster(ext = extent(sppStateRoute), res = c(0.025, 0.025), # same res as above crs = crs(sppStateRoute)@projargs) # using ordinary kriging # SpeciesRichness <- raster::interpolate(object = r, model = predict.model)
## [using ordinary kriging]
If you receive warnings that look similar to below, you have duplicate locations in your data. See above for how to remove duplicate locations.
Covariance matrix singular at location [-114.543,36.7882,0]: skipping...