A common task in GIS is comparing the spatial extent of one layer with another.
Say you have a load of points, some of which overlay a polygon layer. You are only interested
in the points that intersect with the points. What do you do? Also, how can you aggregated-up
the values contained in the points to correspond with the polygons.
These are complex computational problems. In this post
we will see how recent updates to R's
sp package make
the solution surprisingly intuitive and incredibly terse.
Loading the data
All of the data (and more) for this can be downloaded from the
tutorial page on GitHub.
To make this tutorial reproducible on any computer, we will download each input dataset
from within R using
library(sp) download.file("http://robinlovelace.net/data/lnd.RData", destfile = "lnd.RData") download.file("http://robinlovelace.net/data/stations.RData", destfile = "stations.RData") load("lnd.RData") load("stations.RData") plot(stations[sample(1:nrow(stations), 500), ]) plot(lnd, add = T, col = "red")
Spatial subsetting (clipping)
As the plot demonstrates, the stations are far more exentsive than polygons of
central London. We must therefore clip them. Doing this manually would take
much time – we'd have to interrogate the coordinates of each point to see
whether or not it falls within one of the polygon boundaries.
Fortunately with have the
over function from the
sp package to make this
operation more terse:
sel <- over(stations, lnd) stations <- stations[!is.na(sel[, 1]), ]
As if that weren't enough, the developers of
sp have integrated
spatial subsetting into R's main index system with the square brackets.
Because this is a common procedure it is actually possible
to perform it with a single line of code:
stations <- stations[lnd, ] plot(stations)
As the figure shows, only stations within the London borroughs are now shown.
All that was needed was to place another spatial object in the row
index of the points (
[lnd, ]) and R automatically understood that a
subset based on location should be produced. This line of code is an example
of R's 'terseness' – only a single line of code is needed to perform what
is in fact quite a complex operation.
The third way to acheive the
same result uses the
This is more complex and not included in this tutorial
(interested readers can see a vignette of this, to accomany the tutorial
The next section demonstrates
spatial aggregation, a more advanced version of spatial subsetting.
As with R's very terse code for spatial subsetting, the base function
aggregate (which provides summaries of variables based on some grouping variable)
also behaves differently when the inputs are spatial objects.
stations.c <- aggregate(stations, lnd, length) stations.c@data[, 1]
##  48 22 43 18 12 13 25 24 12 46 18 20 28 32 38 19 30 25 31 7 10 38 12 ##  16 28 17 16 28 4 6 14 26 5
The above code performs a number of steps in just one line:
lndpolygon (borrough) each
stationis located in and groups them accordingly
- it counts the number of stations in each borrough
- a new spatial object is created and assigned the name
stations.c, the count of stations
As shown below, the spatial implementation of
aggregate can provide summary statistics of variables.
In this case we take the variable
NUMBER and find its mean value for the stations in each ward.
stations.m <- aggregate(stations[c("NUMBER")], by = lnd, FUN = mean)
For an optional advanced task, let us analyse and plot the result.
q <- cut(stations.m$NUMBER, breaks = c(quantile(stations.m$NUMBER)), include.lowest = T) summary(q)
## [1.82e+04,1.94e+04] (1.94e+04,1.99e+04] (1.99e+04,2.05e+04] ## 9 8 8 ## (2.05e+04,2.1e+04] ## 8
clr <- as.character(factor(q, labels = paste0("grey", seq(20, 80, 20)))) plot(stations.m, col = clr) legend(legend = paste0("q", 1:4), fill = paste0("grey", seq(20, 80, 20)), "topright")
areas <- sapply(stations.m@polygons, function(x) x@area)
This results in a simple choropleth map and a new vector containing the area of each
borrough. As an additional step, try comparing the mean area of each borrough with the
mean value of stations within it:
If you'd like to learn more about R's rapidly improving spatial functionality,
you can download the complete tutorial, in .pdf or .Rmd form, from