{"id":841,"date":"2015-04-01T10:20:00","date_gmt":"2015-04-01T17:20:00","guid":{"rendered":"http:\/\/homepages.uc.edu\/~yaozo\/wordpress\/?p=841"},"modified":"2015-04-01T10:20:00","modified_gmt":"2015-04-01T17:20:00","slug":"crime-data-mapping","status":"publish","type":"post","link":"https:\/\/zhuoyao.net\/index.php\/2015\/04\/01\/crime-data-mapping\/","title":{"rendered":"Crime Data Mapping"},"content":{"rendered":"
For this workshop session we’ll read in some data and do some basic mapping.<\/p>\n

We need a number of packages for this.<\/p>\n

\n
\n
\n
library<\/a>(<\/span>\"sp\"<\/span>)<\/span>\nlibrary<\/a>(<\/span>\"raster\"<\/span>)<\/span>\nlibrary<\/a>(<\/span>\"rgdal\"<\/span>)<\/span>\nlibrary<\/a>(<\/span>\"rgeos\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n
All our data files are in a folder called Data<\/code> at the same level as the current folder:<\/p>\n
\n
\n
\n
Datadir<\/span> =<\/span> \"..\/Data\"<\/span>\nDatafile<\/span> =<\/span> function<\/span>(<\/span>f<\/span>)<\/span> {<\/span>\n    file.path<\/a>(<\/span>Datadir<\/span>,<\/span> f<\/span>)<\/span>\n}<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n
You should make sure Datadir<\/code> is set to the place you’ve put the data files. If you have the tcltk<\/code> package, then you can use as.character(tcltk::tkchooseDirectory())<\/code>. Don’t put a trailing slash or backslash at the end of Datadir<\/code> or files might not be found.<\/p>\n<\/section>\n
\n
Read The Street-level Crime Data<\/h1>\n
The crime data at the street level is supplied in a CSV file each month. We’ll read in one month and plot it. The coordinates are Ordnance Survey of Great Britain grid references, which is EPSG code 27700<\/a>.<\/p>\n
\n
\n
\n
streetCrime<\/span> =<\/span> read.csv<\/a>(<\/span>Datafile<\/a>(<\/span>\"2012-03-cumbria-street.csv\"<\/span>)<\/span>)<\/span>\nhead<\/a>(<\/span>streetCrime<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
##     Month          Reported.by         Falls.within Easting Northing                          Location Crime.type Context\n## 1 2012-03 Cumbria Constabulary Cumbria Constabulary  351534   492353          On or near Dowker'S Lane   Burglary      NA\n## 2 2012-03 Cumbria Constabulary Cumbria Constabulary  299210   528515            On or near Supermarket   Burglary      NA\n## 3 2012-03 Cumbria Constabulary Cumbria Constabulary  351183   492741        On or near Serpentine Road   Burglary      NA\n## 4 2012-03 Cumbria Constabulary Cumbria Constabulary  301277   524712          On or near Hallwood Road   Burglary      NA\n## 5 2012-03 Cumbria Constabulary Cumbria Constabulary  299416   527779 On or near Sports\/Recreation Area   Burglary      NA\n## 6 2012-03 Cumbria Constabulary Cumbria Constabulary  361096   478662            On or near Chapel Lane   Burglary      NA\n<\/pre>\n<\/div>\n
\n
coordinates<\/a>(<\/span>streetCrime<\/span>)<\/span> =<\/span> ~<\/span>Easting<\/span> +<\/span> Northing<\/span>\nproj4string<\/a>(<\/span>streetCrime<\/span>)<\/span> =<\/span> CRS<\/a>(<\/span>\"+init=epsg:27700\"<\/span>)<\/span>\nplot<\/a>(<\/span>streetCrime<\/span>)<\/span>\ntitle<\/a>(<\/span>\"Street-level crime locations\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
OSM Background<\/h2>\n
Can we do a nicer map using the OpenStreetMap<\/code> package? This gives us a function to download and print map tiles from Open Street Map<\/a>, and other tile servers. These tiles use the Google Mercator coordinate system, so we have to transform our points to that. The osm()<\/code> function returns the CRS string for the tiles.<\/p>\n
\n
\n
\n
library<\/a>(<\/span>\"OpenStreetMap\"<\/span>)<\/span>\nlatLongBox<\/span> =<\/span> bbox<\/a>(<\/span>spTransform<\/a>(<\/span>streetCrime<\/span>,<\/span> CRS<\/a>(<\/span>\"+init=epsg:4326\"<\/span>)<\/span>)<\/span>)<\/span>\nmap<\/span> =<\/span> openmap<\/a>(<\/span>c<\/a>(<\/span>latLongBox<\/span>[<\/span>2<\/span>,<\/span> 2<\/span>]<\/span>,<\/span> latLongBox<\/span>[<\/span>1<\/span>,<\/span> 1<\/span>]<\/span>)<\/span>,<\/span> c<\/a>(<\/span>latLongBox<\/span>[<\/span>2<\/span>,<\/span> 1<\/span>]<\/span>,<\/span> latLongBox<\/span>[<\/span>1<\/span>,<\/span> 2<\/span>]<\/span>)<\/span>,<\/span>\n    type<\/span> =<\/span> \"osm\"<\/span>)<\/span>\nplot<\/a>(<\/span>map<\/span>)<\/span>\npoints<\/a>(<\/span>spTransform<\/a>(<\/span>streetCrime<\/span>,<\/span> osm<\/a>(<\/span>)<\/span>)<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
Using ggplot<\/code><\/h2>\n
If we add the coordinates back to the data frame we can use ggplot<\/code> to produce some nice plots:<\/p>\n
\n
\n
\n
library<\/a>(<\/span>\"ggplot2\"<\/span>)<\/span>\nstreetCrime<\/span>@<\/span>data<\/span> =<\/span> cbind<\/a>(<\/span>streetCrime<\/span>@<\/span>data<\/span>,<\/span> coordinates<\/a>(<\/span>streetCrime<\/span>)<\/span>)<\/span>\nggplot<\/a>(<\/span>streetCrime<\/span>@<\/span>data<\/span>)<\/span> +<\/span> geom_point<\/a>(<\/span>aes<\/a>(<\/span>x<\/span> =<\/span> Easting<\/span>,<\/span> y<\/span> =<\/span> Northing<\/span>,<\/span> colour<\/span> =<\/span> Crime.type<\/span>)<\/span>)<\/span> +<\/span>\n    coord_equal<\/a>(<\/span>)<\/span> +<\/span> opts<\/a>(<\/span>title<\/span> =<\/span> \"Street Crime By Type\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n<\/section>\n
\n
Relative Rates<\/h1>\n
It is clear from the maps that there are more crimes in urban areas, because that’s where people live. What we really want to do is look at rates of crime relative to some baseline. We could use population if we had a good measure of baseline population, or we can compare one crime rate against all crimes.<\/p>\n
The sparr<\/code> package provides functions for computing estimates of densities of point patterns using kernel smoothing. It has several methods for computing the kernel bandwidth, including adaptive methods. To use it we have to convert our points into an ppp<\/code> object from the spatstat<\/code> package. For the window containing the points we’ll use the convex hull. We’ll write a function that does the smoothing for a particular crime type. The function will convert the output from bivariate.density<\/code> into a raster<\/code> object. We’ll also read in some town locations and get the largest of them for some context when mapping.<\/p>\n
\n
\n
\n
library<\/a>(<\/span>\"spatstat\"<\/span>)<\/span>\nlibrary<\/a>(<\/span>\"sparr\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
## Warning message: RGL: unable to open X11 display\n<\/pre>\n<\/div>\n
\n
## Warning message: error in rgl_init\n<\/pre>\n<\/div>\n
\n
library<\/a>(<\/span>\"raster\"<\/span>)<\/span>\n\ncoords<\/span> =<\/span> coordinates<\/a>(<\/span>streetCrime<\/span>)<\/span>\ncHull<\/span> =<\/span> as.matrix<\/a>(<\/span>coords<\/span>[<\/span>rev<\/a>(<\/span>chull<\/a>(<\/span>coords<\/span>)<\/span>)<\/span>,<\/span> ]<\/span>)<\/span>\ndata<\/span> =<\/span> ppp<\/a>(<\/span>x<\/span> =<\/span> coords<\/span>[<\/span>,<\/span> 1<\/span>]<\/span>,<\/span> y<\/span> =<\/span> coords<\/span>[<\/span>,<\/span> 2<\/span>]<\/span>,<\/span> window<\/span> =<\/span> owin<\/a>(<\/span>poly<\/span> =<\/span> cHull<\/span>)<\/span>,<\/span> marks<\/span> =<\/span> streetCrime<\/span>$<\/span>Crime.type<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
## Warning message: data contain duplicated points\n<\/pre>\n<\/div>\n
\n
makeIntensity<\/span> <-<\/span> function<\/span>(<\/span>data<\/span>,<\/span> type<\/span>)<\/span> {<\/span>\n    d<\/span> =<\/span> bivariate.density<\/a>(<\/span>data<\/span> =<\/span> data<\/span>,<\/span> ID<\/span> =<\/span> type<\/span>,<\/span> intensity<\/span> =<\/span> TRUE<\/span>,<\/span> comment<\/span> =<\/span> FALSE<\/span>,<\/span> pilotH<\/span> =<\/span> diff<\/a>(<\/span>range<\/a>(<\/span>data<\/span>$<\/span>x<\/span>)<\/span>)<\/span>\/<\/span>10<\/span>)<\/span>\n    raster<\/a>(<\/span>list<\/a>(<\/span>x<\/span> =<\/span> d<\/span>$<\/span>X<\/span>,<\/span> y<\/span> =<\/span> d<\/span>$<\/span>Y<\/span>,<\/span> z<\/span> =<\/span> d<\/span>$<\/span>Zm<\/span>)<\/span>)<\/span>\n}<\/span>\n\nallCrime<\/span> =<\/span> makeIntensity<\/a>(<\/span>data<\/span>,<\/span> NULL)<\/span>\nasbos<\/span> =<\/span> makeIntensity<\/a>(<\/span>data<\/span>,<\/span> \"Anti-social behaviour\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
## Warning message: more than two distinct ID values detected in 'data'\n<\/pre>\n<\/div>\n
\n
drugs<\/span> =<\/span> makeIntensity<\/a>(<\/span>data<\/span>,<\/span> \"Drugs\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
## Warning message: more than two distinct ID values detected in 'data'\n<\/pre>\n<\/div>\n
\n
ukgaz<\/span> =<\/span> readOGR<\/a>(<\/span>Datadir<\/span>,<\/span> \"ukgaz\"<\/span>,<\/span> verbose<\/span> =<\/span> FALSE<\/span>)<\/span>\ntopTowns<\/span> =<\/span> ukgaz<\/span>[<\/span>order<\/a>(<\/span>ukgaz<\/span>$<\/span>pop<\/span>,<\/span> decreasing<\/span> =<\/span> TRUE<\/span>)<\/span>[<\/span>1<\/span>:<\/span>8<\/span>]<\/span>,<\/span> ]<\/span>\n\ncolours<\/span> =<\/span> colorRampPalette<\/a>(<\/span>c<\/a>(<\/span>\"blue\"<\/span>,<\/span> \"gray\"<\/span>,<\/span> \"red\"<\/span>)<\/span>)<\/span>(<\/span>20<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n
ASBO ratio<\/h2>\n
\n
\n
\n
plot<\/a>(<\/span>asbos<\/span>\/<\/span>allCrime<\/span>,<\/span> col<\/span> =<\/span> colours<\/span>)<\/span>\ntext<\/a>(<\/span>topTowns<\/span>,<\/span> topTowns<\/span>$<\/span>Name<\/span>)<\/span>\ntitle<\/a>(<\/span>\"ASBOs as fraction of all crime\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
The regional rate of ASBO crime to all crime is 0.5248<\/code>. We can see that Kendal seems higher and Workington lower than this average value.<\/p>\n
Drugs Ratio<\/h2>\n
\n
\n
\n
plot<\/a>(<\/span>drugs<\/span>\/<\/span>allCrime<\/span>,<\/span> col<\/span> =<\/span> colours<\/span>)<\/span>\ntext<\/a>(<\/span>topTowns<\/span>,<\/span> topTowns<\/span>$<\/span>Name<\/span>)<\/span>\ntitle<\/a>(<\/span>\"Drug offences as fraction of all crime\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
This shows almost the opposite pattern. The average drug crime ratio is 0.0333<\/code>, with Workington over double this, and Kendal much lower.<\/p>\n
These kind of ‘heat maps’ are aften used to show crime maps, house fires, or other point-based data, but without reference to some kind of background population as a denominator they can’t display meaningful risk. The sparr<\/code> package has functions for computing relative risk with error calculations. Read the package documentation for more information.<\/p>\n<\/section>\n
\n
Police Neighbourhood Areas<\/h1>\n
UK crime data is organised by areas known as neighbourhoods. The boundaries of the neighbourhoods are available as a shapefile. The coordinates are lat-long, so we convert to the same CRS as the street crime data. In this section we will aggregate the cases to neighbourhood level.<\/p>\n
\n
\n
\n
neighbourhoods<\/span> =<\/span> spTransform<\/a>(<\/span>readOGR<\/a>(<\/span>Datadir<\/span>,<\/span> \"neighbourhoods\"<\/span>)<\/span>,<\/span> CRS<\/a>(<\/span>proj4string<\/a>(<\/span>streetCrime<\/span>)<\/span>)<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
## OGR data source with driver: ESRI Shapefile \n## Source: \"..\/Data\", layer: \"neighbourhoods\"\n## with 46 features and 3 fields\n## Feature type: wkbPolygon with 2 dimensions\n<\/pre>\n<\/div>\n
\n
plot<\/a>(<\/span>neighbourhoods<\/span>)<\/span>\nplot<\/a>(<\/span>streetCrime<\/span>,<\/span> add<\/span> =<\/span> TRUE<\/span>)<\/span>\ntitle<\/a>(<\/span>\"Neighbourhood polygons\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
Topological Colouring<\/h2>\n
Plain regions like in the previous map never look too interesting, but colouring everything one colour isn’t much better. Can we colour regions so that two regions with a common boundary don’t have the same colour? The theory tells us that we only need at most four colours to do this for a planar map. However, the algorithm to do that is very hard, and our map might not be strictly planar, with islands and holes here and there.<\/p>\n
The following function uses poly2nb<\/code> from the spdep<\/code> package to compute the adjacency of regions, and then executes a ‘greedy’ algorithm. For each region it finds the smallest colour not present in any of its neighbours. If all the colours are currently being used by its neighbours then it has to start using a new one.<\/p>\n
\n
\n
\n
library<\/a>(<\/span>\"spdep\"<\/span>)<\/span>\nlibrary<\/a>(<\/span>\"RColorBrewer\"<\/span>)<\/span>\ncolourmap<\/span> <-<\/span> function<\/span>(<\/span>polys<\/span>,<\/span> snap<\/span>)<\/span> {<\/span>\n    colouring<\/span> =<\/span> list<\/a>(<\/span>)<\/span>\n    length<\/a>(<\/span>colouring<\/span>)<\/span> =<\/span> length<\/a>(<\/span>polys<\/span>)<\/span>\n    nbl<\/span> =<\/span> poly2nb<\/a>(<\/span>polys<\/span>,<\/span> snap<\/span> =<\/span> snap<\/span>)<\/span>\n    nextFreeColour<\/span> =<\/span> 1<\/span>\n    for<\/span> (<\/span>node<\/span> in<\/span> 1<\/span>:<\/span>length<\/a>(<\/span>nbl<\/span>)<\/span>)<\/span> {<\/span>\n        nbrColours<\/span> =<\/span> unlist<\/a>(<\/span>colouring<\/span>[<\/span>nbl<\/span>[[<\/span>node<\/span>]<\/span>]<\/span>]<\/span>)<\/span>\n        freeColour<\/span> =<\/span> min<\/a>(<\/span>setdiff<\/a>(<\/span>(<\/span>1<\/span>:<\/span>nextFreeColour<\/span>)<\/span>,<\/span> nbrColours<\/span>)<\/span>)<\/span>\n        colouring<\/span>[[<\/span>node<\/span>]<\/span>]<\/span> =<\/span> freeColour<\/span>\n        if<\/span> (<\/span>freeColour<\/span> == nextFreeColour<\/span>)<\/span> {<\/span>\n            nextFreeColour<\/span> =<\/span> nextFreeColour<\/span> +<\/span> 1<\/span>\n        }<\/span>\n    }<\/span>\n    unlist<\/a>(<\/span>colouring<\/span>)<\/span>\n}<\/span>\ncolouring<\/span> =<\/span> colourmap<\/a>(<\/span>neighbourhoods<\/span>,<\/span> snap<\/span> =<\/span> 100<\/span>)<\/span>\npalette<\/a>(<\/span>brewer.pal<\/a>(<\/span>max<\/a>(<\/span>colouring<\/span>)<\/span>,<\/span> \"Set3\"<\/span>)<\/span>)<\/span>\nplot<\/a>(<\/span>neighbourhoods<\/span>,<\/span> col<\/span> =<\/span> colouring<\/span>)<\/span>\ntitle<\/a>(<\/span>\"Topology colouring\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
This map shouldn’t have any two adjacent regions with the same colour. The snap<\/code> had to be adjusted to make this work, however,because the map has some problems…<\/p>\n
Digitizing Problems<\/h2>\n
Unfortunately the borders are very poorly digitized, and the edges between neighbourhoods do not match up very well. You can see this by zooming in on Carlisle and plotting with semi-transparent colours. This shows up the overlaps and the gaps.<\/p>\n
\n
\n
\n
plot<\/a>(<\/span>neighbourhoods<\/span>,<\/span> xlim<\/span> =<\/span> c<\/a>(<\/span>336000<\/span>,<\/span> 344000<\/span>)<\/span>,<\/span> ylim<\/span> =<\/span> c<\/a>(<\/span>552000<\/span>,<\/span> 560000<\/span>)<\/span>,<\/span> col<\/span> =<\/span> rainbow<\/a>(<\/span>10<\/span>,<\/span>\n    alpha<\/span> =<\/span> 0.5<\/span>)<\/span>,<\/span> border<\/span> =<\/span> NA<\/span>)<\/span>\ntitle<\/a>(<\/span>\"Overlapping and misaligned polygons\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
Point overlays<\/h2>\n
We can use the over<\/code> function to count the number of crimes in each neighbourhood. This function does various things depending on the exact class of object passed to it. It is also fussy about making sure its arguments have the same coordinate system.<\/p>\n
Because we know that our polygons aren’t well constructed, we’ll return the list of polygons that each point is in. We can then count how many by computing the length of the list element. We’ll then plot Carlisle and show those points that are not over exactly one region.<\/p>\n
\n
\n
\n
crimeN<\/span> =<\/span> over<\/a>(<\/span>SpatialPoints<\/a>(<\/span>streetCrime<\/span>)<\/span>,<\/span> SpatialPolygons<\/a>(<\/span>neighbourhoods<\/span>@<\/span>polygons<\/span>)<\/span>,<\/span> returnList<\/span> =<\/span> TRUE<\/span>)<\/span>\nnPoly<\/span> =<\/span> unlist<\/a>(<\/span>lapply<\/a>(<\/span>crimeN<\/span>,<\/span> length<\/span>)<\/span>)<\/span>\ntable<\/a>(<\/span>nPoly<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
## nPoly\n##    0    1    2 \n##    6 4287   57 \n<\/pre>\n<\/div>\n
\n
plot<\/a>(<\/span>neighbourhoods<\/span>,<\/span> xlim<\/span> =<\/span> c<\/a>(<\/span>336000<\/span>,<\/span> 344000<\/span>)<\/span>,<\/span> ylim<\/span> =<\/span> c<\/a>(<\/span>552000<\/span>,<\/span> 560000<\/span>)<\/span>,<\/span> col<\/span> =<\/span> rainbow<\/a>(<\/span>10<\/span>,<\/span>\n    alpha<\/span> =<\/span> 0.5<\/span>)<\/span>,<\/span> border<\/span> =<\/span> NA<\/span>)<\/span>\nplot<\/a>(<\/span>streetCrime<\/span>,<\/span> cex<\/span> =<\/span> 0.25<\/span>,<\/span> pch<\/span> =<\/span> 19<\/span>,<\/span> add<\/span> =<\/span> TRUE<\/span>,<\/span> col<\/span> =<\/span> \"black\"<\/span>)<\/span>\nplot<\/a>(<\/span>streetCrime<\/span>[<\/span>nPoly<\/span> ><\/span> 1<\/span>,<\/span> ]<\/span>,<\/span> add<\/span> =<\/span> TRUE<\/span>,<\/span> pch<\/span> =<\/span> 19<\/span>,<\/span> col<\/span> =<\/span> \"red\"<\/span>)<\/span>\nplot<\/a>(<\/span>streetCrime<\/span>[<\/span>nPoly<\/span> <<\/span> 1<\/span>,<\/span> ]<\/span>,<\/span> add<\/span> =<\/span> TRUE<\/span>,<\/span> pch<\/span> =<\/span> 19<\/span>,<\/span> col<\/span> =<\/span> \"blue\"<\/span>)<\/span>\ntitle<\/a>(<\/span>\"Points not in a single polygon\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
Points in polygons<\/h2>\n
Conversely we can overlay polygons on points and get a list of the same length as the number of polygons. Each element contains the indices of points in that polygon. This way we will double-count points that are in more than one polygon, and drop those outside all the polygons.<\/p>\n
\n
\n
\n
crimeL<\/span> =<\/span> over<\/a>(<\/span>SpatialPolygons<\/a>(<\/span>neighbourhoods<\/span>@<\/span>polygons<\/span>)<\/span>,<\/span> SpatialPoints<\/a>(<\/span>streetCrime<\/span>)<\/span>,<\/span> returnList<\/span> =<\/span> TRUE<\/span>)<\/span>\n# how many crimes per area<\/span>\nunlist<\/a>(<\/span>lapply<\/a>(<\/span>crimeL<\/span>,<\/span> length<\/span>)<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
##  [1] 399 129 120  93  66 100  33  31  86 145 156 244  54  62  62  47  35 126  90 115  74 205  41  65 270 122 190  40  92 127  93  81\n## [33]   7  38  43  53  32  21  24 152  33 142 195  20  38  10\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n
We can then assign this vector to the neighbourhood object and use it to produce some maps. The default colours in spplot<\/code> aren’t particularly good, so we define a new palette from gray to red and set that.<\/p>\n
\n
\n
\n
colours2<\/span> =<\/span> colorRampPalette<\/a>(<\/span>c<\/a>(<\/span>\"gray90\"<\/span>,<\/span> \"red\"<\/span>)<\/span>)<\/span>(<\/span>20<\/span>)<\/span>\nsp.theme<\/a>(<\/span>set<\/span> =<\/span> TRUE<\/span>,<\/span> regions<\/span> =<\/span> list<\/a>(<\/span>col<\/span> =<\/span> colours2<\/span>)<\/span>)<\/span>\nneighbourhoods<\/span>$<\/span>TotalCrime<\/span> =<\/span> unlist<\/a>(<\/span>lapply<\/a>(<\/span>crimeL<\/span>,<\/span> length<\/span>)<\/span>)<\/span>\nneighbourhoods<\/span>$<\/span>TotalCrimeRateK<\/span> =<\/span> 1000<\/span> *<\/span> neighbourhoods<\/span>$<\/span>TotalCrime<\/span>\/<\/span>neighbourhoods<\/span>$<\/span>Population<\/span>\nspplot<\/a>(<\/span>neighbourhoods<\/span>,<\/span> \"TotalCrimeRateK\"<\/span>,<\/span> main<\/span> =<\/span> \"Crimes\/1000ppl\/month from street data\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
neighbourhoods<\/span>@<\/span>data<\/span>[<\/span>order<\/a>(<\/span>neighbourhoods<\/span>$<\/span>TotalCrimeRateK<\/span>,<\/span> decreasing<\/span> =<\/span> TRUE<\/span>)<\/span>[<\/span>1<\/span>:<\/span>5<\/span>]<\/span>,<\/span> c<\/a>(<\/span>\"Name\"<\/span>,<\/span>\n    \"Population\"<\/span>,<\/span> \"TotalCrimeRateK\"<\/span>)<\/span>]<\/span>\n<\/pre>\n<\/div>\n
\n
##                               Name Population TotalCrimeRateK\n## 24            Carlisle City Centre       2042          132.22\n## 23        Botcherby And Durranhill       1334           48.73\n## 0                   Barrow Central      13706           29.11\n## 26                      St. Aidans       6813           27.89\n## 31 Raffles \/ Newtown \/ Willowholme       5172           15.66\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
You’ll notice that most of the county is gray – this is because the distribution is very skew. Try taking a log-transform and mapping that.<\/p>\n
Area Data<\/h2>\n
The monthly area crime data is supplied in a CSV file. We can read this in, join it to the polygons by matching the ID and Neighbourhood columns from each object. The map from this data looks broadly similar to the data from the street-level file.<\/p>\n
\n
\n
\n
monthCrime<\/span> =<\/span> read.csv<\/a>(<\/span>Datafile<\/a>(<\/span>\"2012-03-cumbria-neighbourhood.csv\"<\/span>)<\/span>)<\/span>\ndim<\/a>(<\/span>monthCrime<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
## [1] 46 15\n<\/pre>\n<\/div>\n
\n
monthCrime<\/span>[<\/span>1<\/span>:<\/span>10<\/span>,<\/span> c<\/a>(<\/span>\"Neighbourhood\"<\/span>,<\/span> \"All.crime.and.ASB\"<\/span>)<\/span>]<\/span>\n<\/pre>\n<\/div>\n
\n
##    Neighbourhood All.crime.and.ASB\n## 1         GARS07                33\n## 2         GARS06                99\n## 3         GARS05                66\n## 4         GARS04                88\n## 5         GARS03               119\n## 6         GARS02               125\n## 7         GARS01               407\n## 8         GARY01                 7\n## 9         GARS09                85\n## 10        GARS08                30\n<\/pre>\n<\/div>\n
\n
neighbourhoods<\/span>$<\/span>MonthlyCrime<\/span> =<\/span> monthCrime<\/span>[<\/span>match<\/a>(<\/span>neighbourhoods<\/span>$<\/span>ID<\/span>,<\/span> monthCrime<\/span>$<\/span>Neighbourhood<\/span>)<\/span>,<\/span>\n    \"All.crime.and.ASB\"<\/span>]<\/span>\nneighbourhoods<\/span>$<\/span>MonthlyCrimeRate<\/span> =<\/span> 1000<\/span> *<\/span> neighbourhoods<\/span>$<\/span>MonthlyCrime<\/span>\/<\/span>neighbourhoods<\/span>$<\/span>Population<\/span>\nsp.theme<\/a>(<\/span>set<\/span> =<\/span> TRUE<\/span>,<\/span> regions<\/span> =<\/span> list<\/a>(<\/span>col<\/span> =<\/span> colours2<\/span>)<\/span>)<\/span>\nspplot<\/a>(<\/span>neighbourhoods<\/span>,<\/span> \"MonthlyCrimeRate\"<\/span>,<\/span> main<\/span> =<\/span> \"Crimes\/1000ppl\/month from monthly data\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n<\/section>\n
\n
Land Use\/Cover Data<\/h1>\n
The land use data comes from the EU Corine Land Cover<\/a> data set. This uses a coordinate system designed for use across Europe. We’ll transform our neighbourhoods to this system, and plot them, using two colours for contrast. The land use is a GeoTIFF file, and it comes with a text file describing the land use categories for creating legends.<\/p>\n
The land use classification is hierarchical at three levels. We can read in the legend file and display the top level classes.<\/p>\n
\n
\n
\n
legend<\/span> =<\/span> read.csv<\/a>(<\/span>Datafile<\/a>(<\/span>\"legend.csv\"<\/span>)<\/span>,<\/span> row.names<\/span> =<\/span> \"GRID_CODE\"<\/span>,<\/span> as.is<\/span> =<\/span> TRUE<\/span>)<\/span>\nunique<\/a>(<\/span>legend<\/span>$<\/span>LABEL1<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
## [1] \"Artificial surfaces\"           \"Agricultural areas\"            \"Forest and semi natural areas\" \"Wetlands\"                     \n## [5] \"Water bodies\"                  \"NODATA\"                        \"UNCLASSIFIED\"                 \n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n
Then each of those levels divides down into subcategories. For example, the agricultural areas divides into:<\/p>\n
\n
\n
\n
unique<\/a>(<\/span>legend<\/span>$<\/span>LABEL2<\/span>[<\/span>legend<\/span>$<\/span>LABEL1<\/span> == \"Agricultural areas\"<\/span>]<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
## [1] \"Arable land\"                      \"Permanent crops\"                  \"Pastures\"                        \n## [4] \"Heterogeneous agricultural areas\"\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n
and then those classifications are further subdivided, for example:<\/p>\n
\n
\n
\n
unique<\/a>(<\/span>legend<\/span>$<\/span>LABEL3<\/span>[<\/span>legend<\/span>$<\/span>LABEL2<\/span> == \"Heterogeneous agricultural areas\"<\/span>]<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
## [1] \"Annual crops associated with permanent crops\"                                          \n## [2] \"Complex cultivation patterns\"                                                          \n## [3] \"Land principally occupied by agriculture, with significant areas of natural vegetation\"\n## [4] \"Agro-forestry areas\"                                                                   \n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n
We’ll now read in the land use file itself and plot it with the raster<\/code> package.<\/p>\n
\n
\n
\n
landUse<\/span> =<\/span> raster<\/a>(<\/span>Datafile<\/a>(<\/span>\"cumbriaLandUse.tif\"<\/span>)<\/span>)<\/span>\nplot<\/a>(<\/span>landUse<\/span>)<\/span>\nneighEU<\/span> =<\/span> spTransform<\/a>(<\/span>neighbourhoods<\/span>,<\/span> CRS<\/a>(<\/span>projection<\/a>(<\/span>landUse<\/span>)<\/span>)<\/span>)<\/span>\nplot<\/a>(<\/span>neighEU<\/span>,<\/span> add<\/span> =<\/span> TRUE<\/span>,<\/span> border<\/span> =<\/span> \"white\"<\/span>,<\/span> lwd<\/span> =<\/span> 3<\/span>)<\/span>\nplot<\/a>(<\/span>neighEU<\/span>,<\/span> add<\/span> =<\/span> TRUE<\/span>,<\/span> border<\/span> =<\/span> \"black\"<\/span>,<\/span> lwd<\/span> =<\/span> 1<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
The colours correspond to the land use categories of LABEL3<\/code> above.<\/p>\n
Extracting Values<\/h2>\n
The extract<\/code> function is the workhorse for getting information out of rasters. We will get the land use code for each of our street crime points and add the information to the object. We use as.character<\/code> in order to look up from the row names of the legend data frame.<\/p>\n
\n
\n
\n
luCrime<\/span> =<\/span> unlist<\/a>(<\/span>extract<\/a>(<\/span>landUse<\/span>,<\/span> spTransform<\/a>(<\/span>streetCrime<\/span>,<\/span> CRS<\/a>(<\/span>projection<\/a>(<\/span>landUse<\/span>)<\/span>)<\/span>)<\/span>)<\/span>)<\/span>\nluCrime<\/span>[<\/span>luCrime<\/span> == \"0\"<\/span>]<\/span> <-<\/span> NA<\/span>\nstreetCrime<\/span>$<\/span>landUse<\/span> =<\/span> luCrime<\/span>\nstreetCrime<\/span>$<\/span>landUseCode3<\/span> =<\/span> legend<\/span>[<\/span>as.character<\/a>(<\/span>luCrime<\/span>)<\/span>,<\/span> \"LABEL3\"<\/span>]<\/span>\nstreetCrime<\/span>$<\/span>landUseCode2<\/span> =<\/span> legend<\/span>[<\/span>as.character<\/a>(<\/span>luCrime<\/span>)<\/span>,<\/span> \"LABEL2\"<\/span>]<\/span>\nstreetCrime<\/span>$<\/span>landUseCode1<\/span> =<\/span> legend<\/span>[<\/span>as.character<\/a>(<\/span>luCrime<\/span>)<\/span>,<\/span> \"LABEL1\"<\/span>]<\/span>\nplot<\/a>(<\/span>streetCrime<\/span>,<\/span> col<\/span> =<\/span> legend<\/span>[<\/span>streetCrime<\/span>$<\/span>landUse<\/span>,<\/span> \"Colour\"<\/span>]<\/span>,<\/span> pch<\/span> =<\/span> 19<\/span>)<\/span>\ntitle<\/a>(<\/span>\"Street level crime by land use\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
We can also use ggplot<\/code> to produce maps by crime type, coloured by broad land use code<\/p>\n
\n
\n
\n
ggplot<\/a>(<\/span>streetCrime<\/span>@<\/span>data<\/span>)<\/span> +<\/span> geom_point<\/a>(<\/span>aes<\/a>(<\/span>x<\/span> =<\/span> Easting<\/span>,<\/span> y<\/span> =<\/span> Northing<\/span>,<\/span> col<\/span> =<\/span> landUseCode1<\/span>)<\/span>)<\/span> +<\/span>\n    coord_equal<\/a>(<\/span>)<\/span> +<\/span> facet_wrap<\/a>(<\/span>~<\/span>Crime.type<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n<\/section>\n
\n
Land Use Classification<\/h1>\n
Some neighbourhoods are clearly more rural than others. By considering the pixels that each polygon covers we can come up with a simple classification scheme.<\/p>\n
Polygon-Raster Overlay<\/h2>\n
The extract<\/code> function considers all polygons in a Spatial Polygons Data Frame and returns a list where the elements are the pixels in each polygon. We’ll then process this list to get a table of land use categories for each neighbourhood.<\/p>\n
The lutable<\/code> function first removes any pixels in the ‘missing’ category, and then computes the fraction of the pixels in each category. We then make the list into a matrix with row and column names so we can then do a stacked barplot. Note there’s not enough room for all the labels, so only a few are shown.<\/p>\n
\n
\n
\n
ext<\/span> =<\/span> extract<\/a>(<\/span>landUse<\/span>,<\/span> neighEU<\/span>)<\/span>\n\nlutable<\/span> =<\/span> function<\/span>(<\/span>x<\/span>,<\/span> ...<\/span>)<\/span> {<\/span>\n    x<\/span> =<\/span> x<\/span>[<\/span>x<\/span> !=<\/span> 0<\/span>]<\/span>\n    tx<\/span> =<\/span> table<\/a>(<\/span>factor<\/a>(<\/span>x<\/span>,<\/span> levels<\/span> =<\/span> row.names<\/a>(<\/span>legend<\/span>)<\/span>)<\/span>)<\/span>\/<\/span>length<\/a>(<\/span>x<\/span>)<\/span>\n    tx<\/span>\n}<\/span>\n\nclm<\/span> =<\/span> do.call<\/a>(<\/span>cbind<\/span>,<\/span> lapply<\/a>(<\/span>ext<\/span>,<\/span> lutable<\/span>)<\/span>)<\/span>\nrownames<\/a>(<\/span>clm<\/span>)<\/span> =<\/span> legend<\/span>$<\/span>LABEL3<\/span>\ncolnames<\/a>(<\/span>clm<\/span>)<\/span> =<\/span> neighEU<\/span>$<\/span>Name<\/span>\nbarplot<\/a>(<\/span>clm<\/span>,<\/span> col<\/span> =<\/span> as.character<\/a>(<\/span>legend<\/span>$<\/span>Colour<\/span>)<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
We can show the legend in a separate plot.<\/p>\n
\n
\n
\n
par<\/a>(<\/span>mar<\/span> =<\/span> c<\/a>(<\/span>0<\/span>,<\/span> 0<\/span>,<\/span> 0<\/span>,<\/span> 0<\/span>)<\/span>)<\/span>\nplot<\/a>(<\/span>c<\/a>(<\/span>0<\/span>,<\/span> 1<\/span>)<\/span>,<\/span> c<\/a>(<\/span>0<\/span>,<\/span> 1<\/span>)<\/span>,<\/span> type<\/span> =<\/span> \"n\"<\/span>,<\/span> axes<\/span> =<\/span> FALSE<\/span>,<\/span> xlab<\/span> =<\/span> \"\"<\/span>,<\/span> ylab<\/span> =<\/span> \"\"<\/span>)<\/span>\nlegend<\/a>(<\/span>0<\/span>,<\/span> 1<\/span>,<\/span> legend<\/span> =<\/span> substr<\/a>(<\/span>legend<\/span>$<\/span>LABEL3<\/span>,<\/span> 1<\/span>,<\/span> 50<\/span>)<\/span>,<\/span> fill<\/span> =<\/span> as.character<\/a>(<\/span>legend<\/span>$<\/span>Colour<\/span>)<\/span>,<\/span> ncol<\/span> =<\/span> 2<\/span>,<\/span>\n    title<\/span> =<\/span> \"Land Use\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
This shows our regions are either mostly red (urban) or yellow\/green (rural). For our simple classification scheme we’ll look at the number of pixels in the region that are in the broad classes defined in LABEL2<\/code>, and take the class with the most pixels. The classify<\/code> function here takes as its argument a vector of land use class codes and returns the label of the most common value of LABEL2<\/code> for those codes. We’ll then run that on the list of pixel values to get a vector classification that we’ll add to the spatial polygons.<\/p>\n
\n
\n
\n
classify<\/span> <-<\/span> function<\/span>(<\/span>x<\/span>)<\/span> {<\/span>\n    ### x is vector of raster values<\/span>\n    ## remove missing values<\/span>\n    inR<\/span> =<\/span> x<\/span>[<\/span>x<\/span> !=<\/span> 0<\/span>]<\/span>\n    tableValues<\/span> =<\/span> sort<\/a>(<\/span>table<\/a>(<\/span>legend<\/span>[<\/span>as.character<\/a>(<\/span>inR<\/span>)<\/span>,<\/span> \"LABEL2\"<\/span>]<\/span>)<\/span>,<\/span> decreasing<\/span> =<\/span> TRUE<\/span>)<\/span>\n    return<\/a>(<\/span>names<\/a>(<\/span>tableValues<\/span>)<\/span>[<\/span>1<\/span>]<\/span>)<\/span>\n\n}<\/span>\n# example<\/span>\nclassify<\/a>(<\/span>c<\/a>(<\/span>1<\/span>,<\/span> 1<\/span>,<\/span> 1<\/span>,<\/span> 22<\/span>,<\/span> 22<\/span>,<\/span> 43<\/span>,<\/span> 22<\/span>,<\/span> 22<\/span>)<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
## [1] \"Heterogeneous agricultural areas\"\n<\/pre>\n<\/div>\n
\n
neighEU<\/span>$<\/span>Class<\/span> =<\/span> factor<\/a>(<\/span>unlist<\/a>(<\/span>lapply<\/a>(<\/span>ext<\/span>,<\/span> classify<\/span>)<\/span>)<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n
Now we’ll plot the neighbourhoods with the new classes. To get colours for each of the classes we’ll average the colours for the levels in the class:<\/p>\n
\n
\n
\n
meanColour<\/span> <-<\/span> function<\/span>(<\/span>cs<\/span>)<\/span> {<\/span>\n    # convert colours to red-green-blue, and average components<\/span>\n    rgbc<\/span> =<\/span> apply<\/a>(<\/span>col2rgb<\/a>(<\/span>as.character<\/a>(<\/span>cs<\/span>)<\/span>)<\/span>,<\/span> 1<\/span>,<\/span> mean<\/span>)<\/span>\n    # convert r-g-b back to colours<\/span>\n    rgb<\/a>(<\/span>rgbc<\/span>[<\/span>1<\/span>]<\/span>,<\/span> rgbc<\/span>[<\/span>2<\/span>]<\/span>,<\/span> rgbc<\/span>[<\/span>3<\/span>]<\/span>,<\/span> max<\/span> =<\/span> 256<\/span>)<\/span>\n}<\/span>\n\n# use tapply to compute mean colours within groups defined by LABEL2<\/span>\ncolourTable<\/span> =<\/span> tapply<\/a>(<\/span>as.character<\/a>(<\/span>legend<\/span>$<\/span>Colour<\/span>)<\/span>,<\/span> as.character<\/a>(<\/span>legend<\/span>$<\/span>LABEL2<\/span>)<\/span>,<\/span> meanColour<\/span>)<\/span>\ncolours3<\/span> =<\/span> colourTable<\/span>[<\/span>levels<\/a>(<\/span>neighEU<\/span>$<\/span>Class<\/span>)<\/span>]<\/span>\n\nsp.theme<\/a>(<\/span>set<\/span> =<\/span> TRUE<\/span>,<\/span> regions<\/span> =<\/span> list<\/a>(<\/span>col<\/span> =<\/span> colours3<\/span>)<\/span>)<\/span>\nspplot<\/a>(<\/span>neighEU<\/span>,<\/span> \"Class\"<\/span>,<\/span> main<\/span> =<\/span> \"Classification system\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
This map shows us a possible classification of urban and two types of rural neighbourhoods. Notice that the regions on the west coast, although they have some urban content, have more rural land area. The one industrial area in the south contains the shipyards at Barrow.<\/p>\n<\/section>\n
\n
Crimes and Roads<\/h1>\n
In this section we’ll look at crime locations and roads. We have a shapefile of roads with classification in lat-long coordinates which we’ll convert to Ordnance Survey grid reference so we can plot thestreetCrime<\/code> points on it. We’ll also define some colours and widths so the plot looks a bit like the OpenStreetMap maps.<\/p>\n
\n
\n
\n
mainroads<\/span> =<\/span> spTransform<\/a>(<\/span>readOGR<\/a>(<\/span>Datadir<\/span>,<\/span> \"mainroads\"<\/span>)<\/span>,<\/span> CRS<\/a>(<\/span>\"+init=epsg:27700\"<\/span>)<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
## OGR data source with driver: ESRI Shapefile \n## Source: \"..\/Data\", layer: \"mainroads\"\n## with 3486 features and 3 fields\n## Feature type: wkbLineString with 2 dimensions\n<\/pre>\n<\/div>\n
\n
summary<\/a>(<\/span>mainroads<\/span>$<\/span>type<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
##       motorway  motorway_link        primary   primary_link           road      secondary secondary_link       tertiary \n##            198             75            610              9            193            629              3            945 \n##  tertiary_link          trunk     trunk_link \n##             11            757             56 \n<\/pre>\n<\/div>\n
\n
roadColours<\/span> =<\/span> c<\/a>(<\/span>\"#809CC0\"<\/span>,<\/span> \"#809CC0\"<\/span>,<\/span> \"#EC989B\"<\/span>,<\/span> \"#EC989B\"<\/span>,<\/span> \"gray\"<\/span>,<\/span> \"#FDD5A4\"<\/span>,<\/span> \"#FDD5A4\"<\/span>,<\/span>\n    \"gray\"<\/span>,<\/span> \"gray\"<\/span>,<\/span> \"#A9DAA9\"<\/span>,<\/span> \"#A9DAA9\"<\/span>)<\/span>\nroadWidths<\/span> =<\/span> c<\/a>(<\/span>4<\/span>,<\/span> 4<\/span>,<\/span> 3<\/span>,<\/span> 3<\/span>,<\/span> 1<\/span>,<\/span> 2<\/span>,<\/span> 2<\/span>,<\/span> 1<\/span>,<\/span> 1<\/span>,<\/span> 4<\/span>)<\/span>\nplot<\/a>(<\/span>mainroads<\/span>,<\/span> col<\/span> =<\/span> roadColours<\/span>[<\/span>mainroads<\/span>$<\/span>type<\/span>]<\/span>,<\/span> lwd<\/span> =<\/span> roadWidths<\/span>[<\/span>mainroads<\/span>$<\/span>type<\/span>]<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
Next we can compute the distance from each point to the nearest point on the major roads of the road network. First we extract the motorways, trunk roads, and primary roads, then work with that.<\/p>\n
\n
\n
\n
isMajor<\/span> =<\/span> mainroads<\/span>$<\/span>type<\/span> %in%<\/span> c<\/a>(<\/span>\"motorway\"<\/span>,<\/span> \"motorway_link\"<\/span>,<\/span> \"trunk\"<\/span>,<\/span> \"trunk_link\"<\/span>,<\/span> \"primary\"<\/span>,<\/span>\n    \"primary_link\"<\/span>)<\/span>\nmajorRoads<\/span> =<\/span> mainroads<\/span>[<\/span>isMajor<\/span>,<\/span> ]<\/span>\nstreetCrime<\/span>$<\/span>dRoad<\/span> =<\/span> unlist<\/a>(<\/span>lapply<\/a>(<\/span>1<\/span>:<\/span>nrow<\/a>(<\/span>streetCrime<\/span>)<\/span>,<\/span> function<\/span>(<\/span>i<\/span>)<\/span> {<\/span>\n    gDistance<\/a>(<\/span>streetCrime<\/span>[<\/span>i<\/span>,<\/span> ]<\/span>,<\/span> majorRoads<\/span>)<\/span>\n}<\/span>)<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n
We can use ggplot<\/code> to produce some graphics showing how the road distance varies for different crimes. A box plot shows the outliers well.<\/p>\n
\n
\n
\n
x45<\/span> =<\/span> opts<\/a>(<\/span>axis.text.x<\/span> =<\/span> theme_text<\/a>(<\/span>angle<\/span> =<\/span> -<\/span>45<\/span>,<\/span> hjust<\/span> =<\/span> 0<\/span>,<\/span> vjust<\/span> =<\/span> 1<\/span>)<\/span>)<\/span>\nggplot<\/a>(<\/span>streetCrime<\/span>@<\/span>data<\/span>,<\/span> aes<\/a>(<\/span>x<\/span> =<\/span> Crime.type<\/span>,<\/span> y<\/span> =<\/span> dRoad<\/span>)<\/span>)<\/span> +<\/span> x45<\/span> +<\/span> geom_boxplot<\/a>(<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
What can we get from this? Perhaps only that shoplifting happens close to major roads, but thats because shops tend to be close to major roads.<\/p>\n<\/section>\n
\n
Space-Time Regional Data<\/h1>\n
We also have a data file of the crime counts in neighbourhoods for a number of months (We would have more months of data but the neighbourhood boundaries changed, so we can only go back so far). The file has one for each month for each neighbourhood, and lists the number of crimes in each category. This form is ideal for feeding into ggplot<\/code> to get a quick time-series plot.<\/p>\n
\n
\n
\n
allCrime<\/span> =<\/span> read.csv<\/a>(<\/span>Datafile<\/a>(<\/span>\"allCrime.csv\"<\/span>)<\/span>,<\/span> as.is<\/span> =<\/span> TRUE<\/span>)<\/span>\nhead<\/a>(<\/span>allCrime<\/span>[<\/span>,<\/span> 1<\/span>:<\/span>5<\/span>]<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
##     Month                Force Neighbourhood All.crime.and.ASB Burglary\n## 1 2011-09 Cumbria Constabulary        GARS07                41        3\n## 2 2011-09 Cumbria Constabulary        GARS06               141        3\n## 3 2011-09 Cumbria Constabulary        GARS05                56        0\n## 4 2011-09 Cumbria Constabulary        GARS04                82        2\n## 5 2011-09 Cumbria Constabulary        GARS03               119        4\n## 6 2011-09 Cumbria Constabulary        GARS02                90        2\n<\/pre>\n<\/div>\n
\n
# plot time series for each nbrhood - leave the legend out<\/span>\nggplot<\/a>(<\/span>allCrime<\/span>,<\/span> aes<\/a>(<\/span>x<\/span> =<\/span> Month<\/span>,<\/span> y<\/span> =<\/span> Criminal.damage.and.Arson<\/span>,<\/span> group<\/span> =<\/span> Neighbourhood<\/span>,<\/span> colour<\/span> =<\/span> Neighbourhood<\/span>)<\/span>)<\/span> +<\/span>\n    geom_line<\/a>(<\/span>)<\/span> +<\/span> opts<\/a>(<\/span>legend.position<\/span> =<\/span> \"none\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
To plot maps of criminal damage and arson over time, we’ll reformat the data into a matrix of neighbourhood by month and attach it to the neighbourhood map data object. We have to do some work with the months to get them into something that can be used in a formula for the spplot<\/code> function. We’ll just plot the first four months.<\/p>\n
\n
\n
\n
library<\/a>(<\/span>\"reshape2\"<\/span>)<\/span>\ncriminalDamage<\/span> =<\/span> dcast<\/a>(<\/span>allCrime<\/span>,<\/span> Neighbourhood<\/span> ~<\/span> Month<\/span>,<\/span> value.var<\/span> =<\/span> \"Criminal.damage.and.Arson\"<\/span>)<\/span>\n\nmonths<\/span> =<\/span> names<\/a>(<\/span>criminalDamage<\/span>)<\/span>[<\/span>-<\/span>1<\/span>]<\/span>\nnames<\/a>(<\/span>criminalDamage<\/span>)<\/span>[<\/span>-<\/span>1<\/span>]<\/span> =<\/span> format<\/a>(<\/span>as.Date<\/a>(<\/span>sprintf<\/a>(<\/span>\"%s-01\"<\/span>,<\/span> months<\/span>)<\/span>)<\/span>,<\/span> \"%b.%Y\"<\/span>)<\/span>\n\ncriminalDamage<\/span> =<\/span> criminalDamage<\/span>[<\/span>match<\/a>(<\/span>neighbourhoods<\/span>$<\/span>ID<\/span>,<\/span> criminalDamage<\/span>$<\/span>Neighbourhood<\/span>)<\/span>,<\/span> ]<\/span>\nneighbourhoods<\/span>@<\/span>data<\/span> =<\/span> cbind<\/a>(<\/span>neighbourhoods<\/span>@<\/span>data<\/span>,<\/span> criminalDamage<\/span>)<\/span>\nnames<\/a>(<\/span>neighbourhoods<\/span>)<\/span> <-<\/span> make.names<\/a>(<\/span>names<\/a>(<\/span>neighbourhoods<\/span>)<\/span>)<\/span>\nsp.theme<\/a>(<\/span>set<\/span> =<\/span> TRUE<\/span>,<\/span> regions<\/span> =<\/span> list<\/a>(<\/span>col<\/span> =<\/span> colours2<\/span>,<\/span> lwd<\/span> =<\/span> 2<\/span>)<\/span>)<\/span>\nspplot<\/a>(<\/span>neighbourhoods<\/span>,<\/span> c<\/a>(<\/span>\"Sep.2011\"<\/span>,<\/span> \"Oct.2011\"<\/span>,<\/span> \"Nov.2011\"<\/span>,<\/span> \"Dec.2011\"<\/span>)<\/span>,<\/span> main<\/span> =<\/span> \"Criminal Damage\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n
\"\"<\/a><\/div>\n<\/div>\n
Static plots like that are fine for publications in dead-tree format, but for online work we can use technology to make animations.<\/p>\n
Animation<\/h2>\n
The animation<\/code> package provides functions for creating frame-based animations. All you need to do is write the code that produces a plot as each frame of the animation.<\/p>\n
\n
\n
\n
library<\/a>(<\/span>\"animation\"<\/span>)<\/span>\nmonths<\/span> =<\/span> names<\/a>(<\/span>criminalDamage<\/span>)<\/span>[<\/span>-<\/span>1<\/span>]<\/span>\nani.options<\/a>(<\/span>verbose<\/span> =<\/span> FALSE<\/span>,<\/span> outdir<\/span> =<\/span> getwd<\/a>(<\/span>)<\/span>)<\/span>\nsaveHTML<\/a>(<\/span>{<\/span>\n    ani.options<\/a>(<\/span>interval<\/span> =<\/span> 0.5<\/span>,<\/span> verbose<\/span> =<\/span> FALSE<\/span>)<\/span>\n    for<\/span> (<\/span>month<\/span> in<\/span> months<\/span>)<\/span> {<\/span>\n        sp.theme<\/a>(<\/span>set<\/span> =<\/span> TRUE<\/span>,<\/span> regions<\/span> =<\/span> list<\/a>(<\/span>col<\/span> =<\/span> colours2<\/span>,<\/span> lwd<\/span> =<\/span> 2<\/span>)<\/span>)<\/span>\n        print<\/a>(<\/span>spplot<\/a>(<\/span>neighbourhoods<\/span>,<\/span> month<\/span>,<\/span> at<\/span> =<\/span> seq<\/a>(<\/span>0<\/span>,<\/span> 55<\/span>,<\/span> len<\/span> =<\/span> 20<\/span>)<\/span>,<\/span> main<\/span> =<\/span> month<\/span>)<\/span>)<\/span>\n        ani.pause<\/a>(<\/span>)<\/span>\n    }<\/span>\n}<\/span>,<\/span> img.name<\/span> =<\/span> \"crime_ani\"<\/span>,<\/span> single.opts<\/span> =<\/span> \"'utf8': false\"<\/span>,<\/span> autoplay<\/span> =<\/span> FALSE<\/span>,<\/span> interval<\/span> =<\/span> 0.5<\/span>,<\/span> imgdir<\/span> =<\/span> \"crime_dir\"<\/span>,<\/span>\n    htmlfile<\/span> =<\/span> \"crimeAni.html\"<\/span>,<\/span> ani.height<\/span> =<\/span> 400<\/span>,<\/span> ani.width<\/span> =<\/span> 600<\/span>,<\/span> title<\/span> =<\/span> \"Crime animation\"<\/span>)<\/span>\n<\/pre>\n<\/div>\n
\n
## HTML file created at: \/home\/rowlings\/Work\/R\/UseR\/UseR2012\/Tutorial\/WorkshopsHTML\/crimeAni.html\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n
\n
\"\"<\/p>\n
\n
First<\/a>Previous<\/a>Play<\/a>Next<\/a>Last<\/a>\u00a0<\/a>\u00a0<\/a>\u00a0<\/a>\u00a0<\/a>\u00a0<\/a>\u00a0<\/a>\u00a0<\/a><\/span>Loop<\/a>–<\/a>+<\/a><\/span><\/form>\n<\/div>\n<\/div>\n<\/div>\n<\/section>\n
\n
Lessons Learnt<\/h1>\n
    \n
  • Creating spatial point data sets from scratch<\/li>\n
  • Making maps with base graphics, ggplot<\/code>, and spplot<\/code><\/li>\n
  • Kernel Smoothing with package:sparr<\/code><\/li>\n
  • Reading in vector data with readOGR<\/code><\/li>\n
  • Reading in raster data with raster<\/code><\/li>\n
  • Analysing raster data with extract<\/code><\/li>\n
  • Transforming coordinates with spTransform<\/code><\/li>\n
  • Making animations with the animation<\/code> package<\/li>\n<\/ul>\n<\/section>\n","protected":false},"excerpt":{"rendered":"
    For this workshop session we’ll read in some data and do some basic mapping. We need a number of packages for this. library(“sp”) library(“raster”) library(“rgdal”)… <\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[20],"tags":[],"class_list":["post-841","post","type-post","status-publish","format-standard","hentry","category-r"],"_links":{"self":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/posts\/841","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/comments?post=841"}],"version-history":[{"count":0,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/posts\/841\/revisions"}],"wp:attachment":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/media?parent=841"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/categories?post=841"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/tags?post=841"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}