Spatial Operations

Load some data

This is data of police neighbourhoods in Cumbria. It is in lat-long coordinates, so we convert to Ordnance Survey coordinates for analysis. Each neighbourhood has a North/South region code.

require(rgeos)
neighs = readOGR("./Data/", "neighbourhoods")
## OGR data source with driver: ESRI Shapefile 
## Source: "./Data/", layer: "neighbourhoods"
## with 46 features and 4 fields
## Feature type: wkbPolygon with 2 dimensions
summary(neighs)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
##      min    max
## x -3.641 -2.159
## y 54.040 55.189
## Is projected: FALSE 
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
##        ID                 Name      Population      region  
##  GARS01 : 1   Alston Moor   : 1   Min.   : 1334   North:24  
##  GARS02 : 1   Appleby       : 1   1st Qu.: 7001   South:22  
##  GARS03 : 1   Aspatria      : 1   Median :10829             
##  GARS04 : 1   Barrow Central: 1   Mean   :10600             
##  GARS05 : 1   Barrow East   : 1   3rd Qu.:13752             
##  GARS06 : 1   Barrow West   : 1   Max.   :18980             
##  (Other):40   (Other)       :40
plot(neighs, col = ifelse(neighs$region == "North", "blue", "red"))
plot of chunk unnamed-chunk-1
neighs = spTransform(neighs, CRS("+init=epsg:27700"))

Dissolving

To merge all the northern and southern neighbourhoods we use gUnaryUnion from the rgeos package. This returns a SpatialPolygons object, so a little work is needed to turn this into a SpatialPolygonsDataFrame so we can store attributes for each polygon.

regional = gUnaryUnion(neighs, neighs$region)
regional = SpatialPolygonsDataFrame(regional, data.frame(region = c("North", 
    "South")), match.ID = FALSE)
plot(regional)
text(coordinates(regional), label = regional$region)
plot of chunk unnamed-chunk-2

Buffering

Use gBuffer to get an inner or outer buffer of features.

buffs = lapply(c(2000, 1000, -1000), function(x) {
    gBuffer(neighs[17:18, ], width = x)
})
plot(neighs[17:18, ], border = "red", lwd = 2)
plot(buffs[[1]], add = TRUE)
plot(buffs[[2]], add = TRUE)
plot(buffs[[3]], add = TRUE)
plot of chunk unnamed-chunk-3

Multiple buffers

Create separate buffers for features by using byid=TRUE

bz = gBuffer(neighs[17:18, ], byid = TRUE, width = 2000)
plot(bz)
plot(neighs[17:18, ], add = TRUE, border = "red", lwd = 2)
plot of chunk unnamed-chunk-4

Intersections

Now suppose we want to find the area of the border strip 2km wide (on both sides) of these regions. The intersection of the boundaries is almost what we want, but we also have to intersect it with the union of the regions crop out some extra bits.

bzu = gIntersection(bz[1, ], bz[2, ])
plot(bz, lwd = 1)
plot(neighs[17:18, ], add = TRUE, border = "red", lwd = 1)
edgeZone = gIntersection(bzu, gUnaryUnion(neighs[17:18, ]))
plot(edgeZone, add = TRUE, col = "#FF000080")
plot of chunk unnamed-chunk-5
gArea(edgeZone)/1e+06
## [1] 95.09

Simplification

gSimplify uses the well-known Douglas-Peucker algorithm to simplify polygons. We'll create a 5km simplified version of Cumbria.

plot(regional, lwd = 4)
outline = gSimplify(gUnaryUnion(regional), 5000)
plot(outline, add = TRUE, border = "red", lwd = 2)
plot of chunk unnamed-chunk-6

Adjacency

Functions in the spdep package are used to compute adjacency

require(spdep)
adj = poly2nb(neighs)
summary(adj)
## Neighbour list object:
## Number of regions: 46 
## Number of nonzero links: 218 
## Percentage nonzero weights: 10.3 
## Average number of links: 4.739 
## Link number distribution:
## 
##  2  3  4  5  6  7  8  9 10 
##  3 11  8 11  5  5  1  1  1 
## 3 least connected regions:
## 18 24 35 with 2 links
## 1 most connected region:
## 21 with 10 links

Adjacency Maps

You can plot adjacencies by plotting it with a set of coordinates - here we use the polygon centroids:

plot(neighs, col = "#f0e0c0")
plot(adj, coordinates(neighs), add = TRUE, col = "red", lty = 2)
plot of chunk unnamed-chunk-8

Graph Theory

You can convert adjacency lists to graph objects and use the igraph package to analyse them.

require(igraph)
## Loading required package: igraph
g = graph.adjlist(adj)
plot(g)
plot of chunk unnamed-chunk-9

You can then do things like compute shortest paths from one region to another via the adjacencies.

p1 = get.shortest.paths(g, 3, 23)[[1]]
plot(neighs)
plot(neighs[p1, ], col = "#0080FF", add = TRUE)
plot(neighs[c(3, 23), ], col = c("green", "red"), add = TRUE)
plot of chunk unnamed-chunk-10

Spatial Points

We get some otter report data (downloaded from GBIF), transform to Ordnance Survey grid, and plot.

otters = readOGR("./Data", "otters")
## OGR data source with driver: ESRI Shapefile 
## Source: "./Data", layer: "otters"
## with 3298 features and 2 fields
## Feature type: wkbPoint with 2 dimensions
otters = spTransform(otters, CRS("+init=epsg:27700"))
plot(neighs)
plot(otters, add = TRUE)
plot of chunk unnamed-chunk-11

Overlay

over is the main function for point/polygon operations.

Care must be taken since the output from over is very dependent on its exact input types.

Usage:

     over(x, y, returnList = FALSE, fn = NULL, ...)

Methods:
     x = "SpatialPoints", y = "SpatialPolygons" returns a numeric
          vector of length equal to the number of points; the number is
          the index (number) of the polygon of ‘y’ in which a point
          falls; NA denotes the point does not fall in a polygon; if a
          point falls in multiple polygons, the last polygon is
          recorded.

     x = "SpatialPointsDataFrame", y = "SpatialPolygons" equal to the
          previous method, except that an argument ‘fn=xxx’ is allowed,
          e.g. ‘fn = mean’ which will then report a data.frame with the
          mean attribute values of the ‘x’ points falling in each
          polygon (set) of ‘y’

     x = "SpatialPoints", y = "SpatialPolygonsDataFrame" returns a
          data.frame of the second argument with row entries
          corresponding to the first argument

over(points,polys)

Overlaying an SPointsDF with a SPolysDF...

otterN = over(otters, neighs)
dim(otterN)
## [1] 3298    4
summary(otterN)
##        ID                   Name        Population      region    
##  GARY07 : 649   Penrith Rural : 649   Min.   : 1334   North:1958  
##  GARY04 : 301   Longtown      : 301   1st Qu.: 8624   South:1340  
##  GARS08 : 227   High Furness  : 227   Median : 9585               
##  GARY06 : 217   Lakes And Shap: 217   Mean   :10139               
##  GARY02 : 185   Brampton      : 185   3rd Qu.:10622               
##  GARY03 : 181   Dalston       : 181   Max.   :18980               
##  (Other):1538   (Other)       :1538

...returns a plain data frame, one row per otter, with the attributes from the area. We can tabulate this to get the area with the most otter reports.

head(rev(sort(table(otterN$Name))), 10)
## 
##     Penrith Rural          Longtown      High Furness    Lakes And Shap 
##               649               301               227               217 
##          Brampton           Dalston            Wigton South Westmorland 
##               185               181               152               150 
##           Silloth           Keswick 
##               144               131

over(polys,points)

We can compute the number of otter reports in each area this way. The returnList argument returns a list, one element for each region, with the otter row numbers. We can then use sapply to get the length of each element of the list.

We'll wrap this in a function for re-use.

pointCount = function(pts, poly) {
    sapply(over(poly, geometry(pts), returnList = TRUE), length)
}
neighs$otterCount = pointCount(otters, neighs)
spplot(neighs, "otterCount")
plot of chunk unnamed-chunk-14
regional$otterCount = pointCount(otters, regional)
plot(regional)
text(coordinates(regional), label = regional$otterCount)
plot of chunk unnamed-chunk-14

Raster ops

require(raster)
landuse = raster("./Data/cumbriaLandUse.tif")
lucodes = read.csv("./Data/legend.csv", row.names = "GRID_CODE", as.is = TRUE)
lucodes$LABEL3[21] = "Agriculture/Natural"
otters = spTransform(otters, CRS(projection(landuse)))

Raster Sampling

We can now lookup the land use category at the otter report locations

otterLU = extract(landuse, otters)
otterTable = table(otterLU)
names(otterTable) = lucodes[names(otterTable), "LABEL3"]
sort(otterTable)
##                       Airports                   UNCLASSIFIED 
##                              1                              1 
##                      Peat bogs Industrial or commercial units 
##                              2                              6 
##            Agriculture/Natural      Non-irrigated arable land 
##                             16                             22 
##              Coniferous forest                   Mixed forest 
##                             23                             26 
##               Intertidal flats   Sport and leisure facilities 
##                             29                             34 
##                   Water bodies     Discontinuous urban fabric 
##                             49                             58 
##   Complex cultivation patterns            Moors and heathland 
##                             60                             73 
##            Broad-leaved forest             Natural grasslands 
##                            139                            150 
##                       Pastures 
##                           2123