Lancaster University
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"))
neighs = spTransform(neighs, CRS("+init=epsg:27700"))
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)
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)
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)
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")
gArea(edgeZone)/1e+06
## [1] 95.09
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)
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
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)
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)
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)
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)
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 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")
regional$otterCount = pointCount(otters, regional) plot(regional) text(coordinates(regional), label = regional$otterCount)
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)))
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