## @knitr setup1 options(width=132) opts_chunk$set(cache=TRUE,cache.path="crimecache/",autodep=TRUE,fig.path="figure/crime-",fig.width=9,fig.height=9) set.seed(310366) library=function(x){suppressPackageStartupMessages(base:::library(x,character.only=TRUE))} ## @knitr setup build_dep() # figure out dependencies automatically library("xtable") ## @knitr packages library("sp") library("raster") library("rgdal") library("rgeos") ## @knitr setdd Datadir = "../Data" Datafile = function(f){file.path(Datadir,f)} ## @knitr readstreet streetCrime = read.csv(Datafile("2012-03-cumbria-street.csv")) head(streetCrime) coordinates(streetCrime)=~Easting+Northing proj4string(streetCrime)=CRS("+init=epsg:27700") plot(streetCrime);title("Street-level crime locations") ## @knitr osmness library("OpenStreetMap") latLongBox = bbox(spTransform(streetCrime,CRS("+init=epsg:4326"))) map = openmap(c(latLongBox[2,2],latLongBox[1,1]),c(latLongBox[2,1],latLongBox[1,2]),type="osm") plot(map) points(spTransform(streetCrime,osm())) ## @knitr ggmap library("ggplot2") streetCrime@data=cbind(streetCrime@data,coordinates(streetCrime)) ggplot(streetCrime@data)+geom_point(aes(x=Easting,y=Northing,colour=Crime.type))+coord_equal()+opts(title="Street Crime By Type") ## @knitr smoothone library("spatstat") library("sparr") library("raster") coords = coordinates(streetCrime) cHull = as.matrix(coords[rev(chull(coords)),]) data = ppp(x=coords[,1],y=coords[,2],window=owin(poly=cHull), marks=streetCrime$Crime.type) makeIntensity <- function(data,type){ d = bivariate.density(data=data, ID=type, intensity=TRUE, comment=FALSE, pilotH=diff(range(data$x))/10) raster(list(x=d$X,y=d$Y,z=d$Zm)) } allCrime = makeIntensity(data,NULL) asbos = makeIntensity(data,"Anti-social behaviour") drugs = makeIntensity(data,"Drugs") ukgaz = readOGR(Datadir,"ukgaz",verbose=FALSE) topTowns = ukgaz[order(ukgaz$pop,decreasing=TRUE)[1:8],] colours = colorRampPalette(c("blue","gray","red"))(20) ## @knitr asbomap plot(asbos/allCrime,col=colours) text(topTowns,topTowns$Name) title("ASBOs as fraction of all crime") ## @knitr drugsmap plot(drugs/allCrime,col=colours) text(topTowns,topTowns$Name) title("Drug offences as fraction of all crime") ## @knitr unnamed-chunk-1 neighbourhoods = spTransform(readOGR(Datadir,"neighbourhoods"),CRS(proj4string(streetCrime))) plot(neighbourhoods) plot(streetCrime,add=TRUE) title("Neighbourhood polygons") ## @knitr topocolour library("spdep") library("RColorBrewer") colourmap <- function(polys,snap){ colouring = list() length(colouring)=length(polys) nbl = poly2nb(polys,snap=snap) nextFreeColour = 1 for(node in 1:length(nbl)){ nbrColours = unlist(colouring[nbl[[node]] ]) freeColour = min(setdiff((1:nextFreeColour),nbrColours)) colouring[[node]]=freeColour if(freeColour == nextFreeColour){ nextFreeColour = nextFreeColour + 1 } } unlist(colouring) } colouring = colourmap(neighbourhoods,snap=100) palette(brewer.pal(max(colouring),"Set3")) plot(neighbourhoods,col=colouring) title("Topology colouring") ## @knitr unnamed-chunk-2 plot(neighbourhoods,xlim=c(336000,344000),ylim=c(552000,560000),col=rainbow(10,alpha=0.5),border=NA) title("Overlapping and misaligned polygons") ## @knitr unnamed-chunk-3 crimeN = over(SpatialPoints(streetCrime),SpatialPolygons(neighbourhoods@polygons),returnList=TRUE) nPoly = unlist(lapply(crimeN,length)) table(nPoly) plot(neighbourhoods,xlim=c(336000,344000),ylim=c(552000,560000),col=rainbow(10,alpha=0.5),border=NA) plot(streetCrime,cex=0.25,pch=19,add=TRUE,col="black") plot(streetCrime[nPoly>1,],add=TRUE,pch=19,col="red") plot(streetCrime[nPoly<1,],add=TRUE,pch=19,col="blue") title("Points not in a single polygon") ## @knitr unnamed-chunk-4 crimeL = over(SpatialPolygons(neighbourhoods@polygons),SpatialPoints(streetCrime),returnList=TRUE) # how many crimes per area unlist(lapply(crimeL,length)) ## @knitr unnamed-chunk-5 colours2 = colorRampPalette(c("gray90","red"))(20) sp.theme(set=TRUE,regions=list(col=colours2)) neighbourhoods$TotalCrime = unlist(lapply(crimeL,length)) neighbourhoods$TotalCrimeRateK = 1000*neighbourhoods$TotalCrime/neighbourhoods$Population spplot(neighbourhoods,"TotalCrimeRateK",main="Crimes/1000ppl/month from street data") neighbourhoods@data[order(neighbourhoods$TotalCrimeRateK,decreasing=TRUE)[1:5],c("Name","Population","TotalCrimeRateK")] ## @knitr unnamed-chunk-6 monthCrime = read.csv(Datafile("2012-03-cumbria-neighbourhood.csv")) dim(monthCrime) monthCrime[1:10,c("Neighbourhood","All.crime.and.ASB")] neighbourhoods$MonthlyCrime= monthCrime[match(neighbourhoods$ID,monthCrime$Neighbourhood),"All.crime.and.ASB"] neighbourhoods$MonthlyCrimeRate = 1000 * neighbourhoods$MonthlyCrime / neighbourhoods$Population sp.theme(set=TRUE,regions=list(col=colours2)) spplot(neighbourhoods,"MonthlyCrimeRate",main="Crimes/1000ppl/month from monthly data") ## @knitr readlegend legend = read.csv(Datafile("legend.csv"),row.names="GRID_CODE",as.is=TRUE) unique(legend$LABEL1) ## @knitr label2 unique(legend$LABEL2[legend$LABEL1=="Agricultural areas"]) ## @knitr label3 unique(legend$LABEL3[legend$LABEL2=="Heterogeneous agricultural areas"]) ## @knitr readplotlu landUse = raster(Datafile("cumbriaLandUse.tif")) plot(landUse) neighEU = spTransform(neighbourhoods,CRS(projection(landUse))) plot(neighEU,add=TRUE,border="white",lwd=3) plot(neighEU,add=TRUE,border="black",lwd=1) ## @knitr plotcrimelu luCrime = unlist(extract(landUse,spTransform(streetCrime,CRS(projection(landUse))))) luCrime[luCrime=="0"] <- NA streetCrime$landUse = luCrime streetCrime$landUseCode3 = legend[as.character(luCrime),"LABEL3"] streetCrime$landUseCode2 = legend[as.character(luCrime),"LABEL2"] streetCrime$landUseCode1 = legend[as.character(luCrime),"LABEL1"] plot(streetCrime,col=legend[streetCrime$landUse,"Colour"],pch=19) title("Street level crime by land use") ## @knitr ggplot2waymap ggplot(streetCrime@data)+geom_point(aes(x=Easting,y=Northing,col=landUseCode1))+ coord_equal() + facet_wrap(~Crime.type) ## @knitr extractlanduse ext = extract(landUse,neighEU) lutable = function(x,...){ x = x[x!=0] tx = table(factor(x,levels=row.names(legend)))/length(x) tx } clm = do.call(cbind,lapply(ext,lutable)) rownames(clm)=legend$LABEL3 colnames(clm)=neighEU$Name barplot(clm,col=as.character(legend$Colour)) ## @knitr legend par(mar=c(0,0,0,0)) plot(c(0,1),c(0,1),type="n",axes=FALSE,xlab="",ylab="") legend(0,1,legend=substr(legend$LABEL3,1,50),fill=as.character(legend$Colour),ncol=2,title="Land Use") ## @knitr unnamed-chunk-7 classify <- function(x){ ### x is vector of raster values ## remove missing values inR = x[x!=0] tableValues = sort(table(legend[as.character(inR),"LABEL2"]),decreasing=TRUE) return(names(tableValues)[1]) } # example classify(c(1,1,1,22,22,43,22,22)) neighEU$Class = factor(unlist(lapply(ext,classify))) ## @knitr plotnewclasses meanColour <- function(cs){ # convert colours to red-green-blue, and average components rgbc = apply(col2rgb(as.character(cs)),1,mean) # convert r-g-b back to colours rgb(rgbc[1],rgbc[2],rgbc[3],max=256) } # use tapply to compute mean colours within groups defined by LABEL2 colourTable = tapply(as.character(legend$Colour),as.character(legend$LABEL2),meanColour) colours3 = colourTable[levels(neighEU$Class)] sp.theme(set=TRUE,regions=list(col=colours3)) spplot(neighEU,"Class",main="Classification system") ## @knitr roadmap mainroads=spTransform(readOGR(Datadir,"mainroads"),CRS("+init=epsg:27700")) summary(mainroads$type) roadColours = c("#809CC0","#809CC0", "#EC989B","#EC989B", "gray", "#FDD5A4","#FDD5A4", "gray","gray", "#A9DAA9","#A9DAA9") roadWidths = c(4,4,3,3,1,2,2,1,1,4) plot(mainroads,col=roadColours[mainroads$type],lwd=roadWidths[mainroads$type]) ## @knitr roaddistance isMajor = mainroads$type %in% c("motorway","motorway_link","trunk","trunk_link","primary","primary_link") majorRoads = mainroads[isMajor,] streetCrime$dRoad = unlist(lapply(1:nrow(streetCrime),function(i){gDistance(streetCrime[i,],majorRoads)})) ## @knitr ggplotroaddist x45=opts(axis.text.x=theme_text(angle=-45,hjust=0,vjust=1)) ggplot(streetCrime@data,aes(x=Crime.type,y=dRoad))+x45+geom_boxplot() ## @knitr tscrime allCrime = read.csv(Datafile("allCrime.csv"),as.is=TRUE) head(allCrime[,1:5]) # plot time series for each nbrhood - leave the legend out ggplot(allCrime,aes(x=Month,y=Criminal.damage.and.Arson,group=Neighbourhood,colour=Neighbourhood))+geom_line()+opts(legend.position="none") ## @knitr unnamed-chunk-8 library("reshape2") criminalDamage = dcast(allCrime,Neighbourhood~Month,value.var="Criminal.damage.and.Arson") months=names(criminalDamage)[-1] names(criminalDamage)[-1]=format(as.Date(sprintf("%s-01",months)),"%b.%Y") criminalDamage = criminalDamage[match(neighbourhoods$ID,criminalDamage$Neighbourhood),] neighbourhoods@data = cbind(neighbourhoods@data,criminalDamage) names(neighbourhoods) <- make.names(names(neighbourhoods)) sp.theme(set=TRUE,regions=list(col=colours2,lwd=2)) spplot(neighbourhoods,c("Sep.2011","Oct.2011","Nov.2011","Dec.2011"),main="Criminal Damage") ## @knitr crimeanimation library("animation") months = names(criminalDamage)[-1] ani.options(verbose=FALSE,outdir=getwd()) saveHTML({ ani.options(interval = 0.5,verbose=FALSE) for (month in months ) { sp.theme(set=TRUE,regions=list(col=colours2,lwd=2)) print(spplot(neighbourhoods,month,at=seq(0,55,len=20),main=month)) ani.pause() } }, img.name = "crime_ani", single.opts = "'utf8': false", autoplay = FALSE, interval = 0.5, imgdir = "crime_dir", htmlfile = "crimeAni.html", ani.height = 400, ani.width = 600, title = "Crime animation")