Appendix A

Data Sources

A combination of data from different sources was used in this study. Most of the data required some degree of manipulation before it could be used to build the location-allocation model.
All of the original source data was downloaded to the one folder. After manipulation, the useable data was all saved together in a separate folder to the original data, to avoid confusion.

Boundary Files

Boundary files are provided by the Central Statistics Office via their website (Central Statistics Office, 2011). They can be obtained here.

The CSO provide this data in ESRI shapefile format. Download three shapefiles: Administrative counties; Electoral Divisions (ED); and, Small Areas (SA).

Garda Stations

Locations of existing Garda stations were obtained from An Garda Síochána’s station directory web page (An Garda Síochána, 2016) available here.

This data was only available from An Garda Síochána in zipped Keyhole Markup Langaue format (KMZ). To download the KMZ file, copy the URL and paste it into your search bar.

I have subsequently found a shapefile containing Garda station locations produced by Paul Creaner of the Garda Mapping Division which can be found on his public ESRI page here.

However, converting a file from KMZ to KML and reading it into R is still a useful technique. This process is described below.

Road Network

Geofabrik provide road network data for the Republic of Ireland and Northern Ireland in shapefile format. This dataset is created using OpenStreetMap (OSM) data (Geofabrik, 2016).

This data can be accessed here.

Download ireland-and-northern-ireland-latest.shp.zip to your working directory.

Data Preparation

Preparation is divided into three sections here by source of data.

Most of the data preperation was carried out within the R environment for the purpose of reproducibility. However, some techniques required the use of a GIS. ESRI’s ArcGIS was utilised for topology error correction regarding the road network data.

Initial Configuration

# Set working directory to folder containing all original data sources
setwd("G:\\Semester 2\\NCG 606 Thesis\\R data bank\\Kildare 5.0")
dsn <- getwd()
# Load initial packages
library(sp)
library(rgdal)
library(GISTools)
library(rgeos)
library(rgdal)

Boundary Files

Step 1)

Read the shapfiles into objects in R, then select Kildare specific data (e.g. County Kildare from the Administrative counties shapefile).

# Read Administrative counties shapefile into R
County <- readOGR(dsn = dsn, layer = "Census2011_Admin_Counties_generalised20m")
# @data shows attribute data
head(County@data)
# Use COUNTYNAME to select Kildare
Kildare <- County[County@data$COUNTYNAME == "Kildare County", ]

# Read SA shapefile into R
SA <- readOGR(dsn = dsn, layer = "Census2011_Small_Areas_generalised20m")
head(SA@data)
# Create Kildare SA's
KildareSA <- SA[SA@data$COUNTYNAME == "Kildare County", ]
Step 2)

Create a SPatialPointsDataFrame (SPTDF) from the centroids of each Small Area (SA).

# byid=TRUE applies the gCentroid() function to all subgeometries
KildareSA_Centroid <- gCentroid(KildareSA, byid = TRUE)
# Not SpatialPointsDataframe yet, only Spatial Points object. Need to add
# attribute data.
KildareSA_Centroid.spdf <- SpatialPointsDataFrame(KildareSA_Centroid, data = data.frame(KildareSA))
KildareSA_Centroid <- KildareSA_Centroid.spdf

Garda Stations

Step 1)

Read in the KML file containing the locations of Garda stations in the Republic of Ireland.

# Read in Garda station KML file
gstation <- readOGR(dsn = "gstation.kml", layer = "Garda Stations")
names(gstation)
# This reveals there is a description column also. Don't need this, so get
# rid of it:
gstation <- gstation[1]
Step 2)

Before carrying out a point-in-polygon operation, set all projections for datasets loaded so far. It’s best to set up projections before using any geometry functions to select data.

# ESPG:2157 is the IRENET95 Irish Transverse Mercator
newProj <- CRS("+init=epsg:2157")
# Set projected coordinate system to ESPG:2157 for each dataset
KildareTemp <- Kildare
Kildare <- spTransform(KildareTemp, newProj)
KildareSATemp <- KildareSA
KildareSA <- spTransform(KildareSATemp, newProj)
KildareSA_CentroidTemp <- KildareSA_Centroid
KildareSA_Centroid <- spTransform(KildareSA_CentroidTemp, newProj)
gstationTemp <- gstation
gstation <- spTransform(gstationTemp, newProj)
Step 3)

Now, a point-in-polygon operation is performed to select all Garda stations within the Kildare County boundary. A new SpatialPointsDataframe is created out of this selection. This process is outlined in Brunsdon and Comber (2015).

# Use the gIntersection() function to get Kildare specific Garda stations.
kgs <- gIntersection(Kildare, gstation, byid = TRUE)
# This is not a SpatialPointsDataframe yet - neds attribute data Object
# containing string of Administrative county ID (for Kildare) and Garda
# station IDs i.e. '20 27' '20 132' '20 159'
tmp <- rownames(data.frame(kgs))
# Need to get the station IDs.  Result is: '20' '27' etc.
tmp <- strsplit(tmp, " ")
# Select Garda station IDs by their index values.  i.e. '27' '132' '159'
# etc...
gstation.id <- (sapply(tmp, "[[", 2))
gstation.id <- as.numeric(gstation.id)
df <- data.frame(gstation[gstation.id, ])
# We only need the names of Garda stations.
df <- df[1]
# Finally, create the SPDF by combining the dataframe (names) with the
# Spatial Points.
kgs <- SpatialPointsDataFrame(kgs, data = df)

Road Network

Data preparation for the Road Network is described below. Note due to its large size, the roads dataset acquired from Geofabrik takes roughly 5-10 minutes to load into R. Other operations on this dataset are also time consuming due to the nature of its initial size.

Step 1)

Create a 15km buffer around the boundary of Kildare. Read in the roads dataset and then clip it using the 15km polygon buffer around Kildare.

# Create a 15km buffer around Kildare polygon. To check the units of
# measurement:
proj4string(Kildare)
# +units=m : So give the units in gBuffer() in meters
Kildare_Buf <- gBuffer(Kildare, width = 15000)
# Read in the 'roads' dataset. This takes roughly 5-10 mins
road <- readOGR(dsn, "roads")
# Project the road correctly. This takes 5-10 mins
roadTemp <- road
road <- spTransform(roadTemp, newProj)
# Clip roads intersecting with the buffer around Kildare
kroad <- gIntersection(Kildare_Buf, road, byid = TRUE)
Step 2)

Like the kgs object selected from the gstation dataset in Appendix A.3, the kroad object selected from the roads dataset above no longer has associated attribute data. This step outlines the process of combining a dataframe to a Spatial Lines object.

head(names(kroad))
tmp <- names(kroad)
tmp <- strsplit(tmp, " ")
road.id <- (sapply(tmp, "[[", 2))
df <- data.frame(road[road.id, ])
# The row.names of data and Line IDs don't match. Can't subtract 'buffer'
# from names(kroad), so going to add 'buffer' to row.names(df)
head(rownames(df))
rownames(df) <- paste("buffer", rownames(df), sep = " ")
head(rownames(df))
# Now create SpatialLinesDataFrame
kroad.sldf <- SpatialLinesDataFrame(kroad, data = df)
Step 3)

Now that a SpatialLinesDataframe has been created for ‘Kildare specific’ roads, some SQL-styled selections are made on the SLDF to create a subset of it with only the desired road types.

# Get at all types of road in the Kildare Road Network layer
unique(kroad.sldf@data$type)
# There are 26 levels in Kildare Road Network, and 32 overall. To plot
# different road types:
x <- kroad.sldf[kroad.sldf@data$type == "motorway", ]
plot(x, col = "darkred")
title(main = "Motorway Road Segments", font.main = 2, cex.main = 1.5)

# The roads to select; motorway, motorway_link, primary, primary_link, road,
# secondary, secondary_link, services, tertiary, tertiary_link, trunk,
# trunk_link and unclassified.
x1 <- kroad.sldf[kroad.sldf@data$type == "motorway", ]
x2 <- kroad.sldf[kroad.sldf@data$type == "motorway_link", ]
x3 <- kroad.sldf[kroad.sldf@data$type == "primary", ]
x4 <- kroad.sldf[kroad.sldf@data$type == "primary_link", ]
x5 <- kroad.sldf[kroad.sldf@data$type == "road", ]
x6 <- kroad.sldf[kroad.sldf@data$type == "secondary", ]
x7 <- kroad.sldf[kroad.sldf@data$type == "secondary_link", ]
x8 <- kroad.sldf[kroad.sldf@data$type == "services", ]
x9 <- kroad.sldf[kroad.sldf@data$type == "tertiary", ]
x10 <- kroad.sldf[kroad.sldf@data$type == "tertiary_link", ]
x11 <- kroad.sldf[kroad.sldf@data$type == "trunk", ]
x12 <- kroad.sldf[kroad.sldf@data$type == "trunk_link", ]
x13 <- kroad.sldf[kroad.sldf@data$type == "unclassified", ]

# Create a new SLDF with desired road types
Kildare_Road_Network <- rbind(x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, 
    x12, x13)
Step 4)

This is an essential step before creating a location-allocation model. The network will be used to calculate the shortest-path distance matrix, therefore it needs to be connected and not separated into isolated fragments. The main connected portion of the road network is returned to the roadCheck object.

# Load the 'shp2graph' package
library(shp2graph)
# The main self-connected portion of the road network is stored in the
# object 'roadCheck'
roadCheck <- nt.connect(Kildare_Road_Network)

# Save both the Kildare Road Network Unconnected (KRNU) and the Kildare Road
# Network Connected (KRNC) as shapefiles to the working directory. They will
# be further analysed in Appendix A.6 using ArcGIS.
writeLinesShape(roadCheck, "KRNC")
writeLinesShape(Kildare_Road_Network, "KRNU")
Step 5)

Further work was carried out within the environment of ArcGIS to create a more complete road network. This process involved using ArcMap’s Feature To Line (data management) operationto split all lines at their intersections with nodes.

Step 6)

The modified road network was saved as a new shapefile from ArcMap to the working directory. The updated road network shapefile is called RoadNetwork_FeatureToLine1. This shapefile is loaded into R and the connectivity is checked again.

# Load improved road network back into R
Kildare_Road_Network_2 <- readOGR(dsn, "RoadNetwork_FeatureToLine1")

# Check the connectivity. The main connected portion of this road network
# will be stored in the object 'roadCheck2'. A plot of the connectivity of
# the road network is produced again
roadCheck2 <- nt.connect(Kildare_Road_Network_2)

Saving Shapefiles

All of the datasets that have been prepared so far need to be saved for use with ESRI’s ArcGIS Network Analyst extension to perform location-allocaiton analysis.

writeLinesShape(roadCheck2, "Kildare Road Network Complete")
writePolyShape(Kildare, "Kildare")
writePolyShape(KildareSA, "KildareSA")
writePointsShape(KildareSA_Centroid, "KildareSA_Centroid")
writePointsShape(kgs, "Kildare_gStation")

Convert Road Network to Graph

In order to compute the shortest-path distance matrix, the road network must be converted to a graph. This process follows that of Lu (2012).

Step 1)

First, each demand/facility point needs to be mapped to the nearest node on the network.

# Extract the coordinates Existing facilities (15 Garda stations)
ExfacilityXY <- coordinates(kgs)
# Demand/Facilities (731 SA centroids)
demandXY <- coordinates(KildareSA_Centroid)
# Map the points to the nearest nodes on the network Takes roughly 2 mins
ExFpt2nt <- points2network(ntdata = roadCheck2, pointsxy = ExfacilityXY, mapping.method = 1, 
    ELComputed = T)
# Takes roughly 5 mins
Dpt2nt <- points2network(ntdata = roadCheck2, pointsxy = demandXY, mapping.method = 1, 
    ELComputed = T)
Step 2)

Conversion from SpatialLinesDataFrame to an “igraph” object.

# Existing Facilities (Garda stations). IDs of corresponding nodes for data
# points
CoorespondIDs = ExFpt2nt[[3]]
# Generate an 'igraph' object by weighting this graph with edge length
ig <- nel2igraph(ExFpt2nt[[1]], ExFpt2nt[[2]], weight = ExFpt2nt[[8]])
gstation <- as.numeric(CoorespondIDs)
# Demand (SAs)
CoorespondIDs3 = Dpt2nt[[3]]
ig3 <- nel2igraph(Dpt2nt[[1]], Dpt2nt[[2]], weight = Dpt2nt[[8]])
SA <- as.numeric(CoorespondIDs3)

Note that some SA centroid points have mapped to the same node. This was realised due to an error when trying to compute the shortest-path distance matrix between existing facilities (Garda stations) and demand (SA centroids). The error received is due to having duplicate vertices on the same node. A way of checking for this is to use the duplicated() function on the corresponding node IDs for the SA centroids which were mapped to the network, shown below:

duplicated(CoorespondIDs3)

As can be seen, a number of SA centroids have been allocated to the same node on the network. As a work-around, the following code can be implemented to increment the node ID by one and then re-checked until no duplications are left. This is making a huge assumption that: node IDs which are close numerically, are also close in geographical space. Unfortunately, this was the best solution at the time of writing. With the interest of solving a p-median problem using R, this assumption was used.

A <- as.numeric(CoorespondIDs3)
while (anyDuplicated(A) > 0) {
    i <- anyDuplicated(A)
    A[i] <- A[i] + 1
}
SA <- A

Location-Allocation Analysis with R

library(tbart)

The tbart package developed by Brunsdon (2015b) allows the user to solve Teitz and Bart’s (1968) p-median problem.

Given a set of points, Teitz and Bart’s p-median algorithm attempts to find the subset of size p such that summed distances of any point in the set to the nearest point in p is minimised. Although generally effective, this algorithm does not guarantee that a globally optimal subset is found.

Step 1)

Set up distance matrices.

# Euclidean distance matrix Between demand and existing Garda stations
EucDist1 <- euc.dists(KildareSA_Centroid, kgs)
# Between demand and candidate facilities
EucDist2 <- euc.dists(KildareSA_Centroid, KildareSA_Centroid)

# Shortest-path distance matrix between demand and existing facilities
spm1 <- shortest.paths(ig, v = c(SA), to = c(gstation), weight = get.edge.attribute(ig, 
    name = "weight"))
Step 2)

Create two p-median models, where p = 15 (i.e. where 15 facilities are to be sited).

Model_1 will use the kgs SPDF as candidate facilities. This only has 15 locations to choose from, hence this model will show the allocation of demand to existing Garda stations as it presently is. Model_2 will use the SA centroids as both candidate facilities and demand. This model will show the allocation of demand to optimally located Garda stations.

# Existing locations
model_1 <- allocations(KildareSA_Centroid, kgs, p = 15)
# Create the spider diagram linking existing Garda stations to SA centroids
# (demand)
star.model_1 <- star.diagram(KildareSA_Centroid, kgs, alloc = model_1$allocation)
# Optimal locations
model_2 <- allocations(KildareSA_Centroid, KildareSA_Centroid, p = 15)
# Create the spider diagram linking optimal locations to SA centroids
# (demand)
star.model_2 <- star.diagram(KildareSA_Centroid, KildareSA_Centroid, alloc = model_2$allocation)

The facilities which each demand point are allocated to are returned to both model_1 and model_2 objects, along with the allocation distances. By calculating the summed allocation distance, we can measure how efficient the ‘optimal’ solution is compared to the ‘existing solution’.

# Calculate the total distance travelled between demand and facility given
# the present spatial distribution of Garda stations
sum(model_1$allocdist)
# Calculate the total distance travelled between demand and facility when
# Garda stations are optimally located
sum(model_2$allocdist)
Step 3)

Create a SpatialPointsDataframe from the optimal locations obtained from model_2.

# Returns the index number of optimally chosen facilities
optimal_loc <- unique(model_2$allocation)
optimal_loc <- KildareSA_Centroid[optimal_loc, ]
Step 4)

Plot the locations of existing Garda stations against the optimally chosen sites.

par(mar = c(5, 0, 5, 0))
par(mfrow = c(1, 2))
# Plot existing
plot(Kildare)
plot(kgs, col = "royalblue3", pch = 19, add = T)
# Labels
ExName <- kgs$Name
ExLat <- data.frame(kgs)[, 2]
ExLon <- data.frame(kgs)[, 3]
PlEx <- pointLabel(ExLat, ExLon, ExName, offset = 0, cex = 0.7)
title(main = "Existing Garda stations", font.main = 6)
map.scale(694000, 696455.8, 10000, "Kilometers", 2, 5)
# plot optimal
plot(Kildare)
plot(optimal_loc, col = "orange", cex = 1.5, lwd = 2, add = T)
# Labels
OpName <- optimal_loc$EDNAME
OpLon <- optimal_loc$X
OpLat <- optimal_loc$Y
PlOp <- pointLabel(OpLon, OpLat, OpName, offset = 0, cex = 0.7)
title(main = "Optimal Garda stations", font.main = 6)
par(mfrow = c(1, 1))
north.arrow(700021.4, 693487.2, km2ft(0.25), col = "gray20", cex = 0.7)

Step 5)

Plot the allocation of demand to both the existing and optimal spatial distribution of Garda stations.

par(mar = c(5, 0, 5, 0))
par(mfrow = c(1, 2))
# Existing
plot(Kildare)
plot(star.model_1, add = T)
plot(kgs, col = "royalblue3", pch = 19, add = T)
title(main = "Allocation of demand to existing\n Garda stations", font.main = 6)
map.scale(694000, 696455.8, 10000, "Kilometers", 2, 5)
# Optimal
plot(Kildare)
plot(star.model_2, add = T)
plot(optimal_loc, col = "orange", cex = 1.5, lwd = 2, add = T)
title(main = "Allocation of demand to optimally\n located Garda stations", font.main = 6)
par(mfrow = c(1, 1))
north.arrow(700021.4, 693487.2, km2ft(0.25), col = "gray20", cex = 0.7)

Step 6)

It would be helpful to see the third map blown up on its own to see how close the optimal v.s. existing locations are.

par(mar = c(1, 0, 5, 0))
plot(Kildare)
plot(optimal_loc, col = "orange", pch = 3, lwd = 2, add = T)
plot(kgs, col = "royalblue3", pch = 20, lwd = 1, add = T)
title(main = "Optimal and Existing locations of Garda stations,\nin Co. Kildare", 
    font.main = 6)
legend(687000, 691328, legend = c("Existing locations", "Optimal locations"), 
    bty = "n", cex = 0.9, col = c("royalblue3", "orange"), pch = c(20, 3))
north.arrow(694000, 697000, miles2ft(0.2), col = "gray20", cex = 0.9, add = T)

Appendix C

Appendix C contains twelve choropleth maps displaying the distribution of various crimes recorded at Garda Sub-District level for Co. Kildare in 2011. The maps were produced using R with data obtained from AIRO (reference).

This document was created using R Markdown which ensures reproducibility by knitting the R code into the document. Therefore the code has to run correctly before the document will publish.

The twelve crime types are listed as follows:

Category 1) Attempts or threats to murder, assaults, harassments and related offences;

Category 2) Dangerous or negligent acts;

Category 3) Kidnapping and related offences;

Category 4) Robbery, extortion and hijacking offences;

Category 5) Burglary and related offences;

Category 6) Theft and related offences;

Category 7) Fraud, deception and related offences;

Category 8) Controlled drug offences;

Category 9) Weapons and explosives offences;

Category 10) Damage to property and to the environment;

Category 11) Public order and other social code offences;

Category 12) Offences against government, justice procedures and organisation of crime.

Data Sources

Some more data was required to create the crime maps for Co. Kildare. Crime data is provided by AIRO (2016) at Garda Sub-District (GSD) level for the years 2004 - 2016.

AIRO Crime Data

The crime data can be obtained from the AIRO datastore here.

The CSV file was downloaded to the working directory.

Garda Sub-Districts

The Crime data is provided at GSD level, therefore the GSD boundary file must be downloaded from the CSO (2011) here.

The GSD zip file was downloaded - Census2011_Garda_SubDistricts_Nov2013.zip (27 MB).

Data Preparation

AIRO Crime Data

Step 1

Subset the AIRO dataframe so that only Kildare-specific data is extracted for the year 2011.

# Read in the CSV file
GardaCSV <- read.csv("garda_stations_crime.csv", stringsAsFactors = F)
# Match the data with the corresponding station by station name, as IDs
# don't match
GardaCSVmatch <- match(kgs$Name, GardaCSV$Station)

# Change 'Kilcullen' to 'Killcullen' in the CSV dataframe
GardaCSV$Station[325] <- paste("Killcullen")
# Now they match
GardaCSVmatch <- match(kgs$Name, GardaCSV$Station)
GardaCSVmatch
Step 2)
# Subset GardaCSV to get KildareGardaCSV using GardaCSVmatch
KildareGardaCSV <- GardaCSV[GardaCSVmatch, ]

# Extract only the data for 2011
df <- as.data.frame(KildareGardaCSV)
# Create a dataframe with all 2011 crime data
Crime2011 <- df[, grepl("2011", names(df))]
# Create a dataframe with all Kildare stations
KildareCrime <- df[1:5]
# Bind these two df together to create a KildareCrime2011 dataframe
KildareCrime2011.df <- cbind(KildareCrime, Crime2011)

Garda Sub-Districts

Subset the GSD boundary file to extract Kildare-specific GSDs.

# Read in the Garda Sub-Districts shapefile
GSD <- readOGR(dsn, "Census2011_Garda_Subdistricts_Nov2013")
# @data shows attribute data
names(GSD@data)
# Use DIVISION to select Kildare
KildareGSD <- GSD[GSD@data$DIVISION == "Kildare", ]
Step 1)

Create a SpatialPolygonsDataframe (SPDF) by combining both modified datasets.

# Change the rownames in KildareCrime2011.df and KildareGSD to letters so
# they match
rownames(KildareCrime2011.df) <- letters[1:15]
row.names(KildareGSD) <- letters[1:15]
# Create SpatialPolygonsDataFrame
KildareCrime2011.spdf <- SpatialPolygonsDataFrame(KildareGSD, data = KildareCrime2011.df)
KildareCrime2011 <- KildareCrime2011.spdf

# Project the shapefile correctly
newProj <- CRS("+init=epsg:2157")
KildareCrime2011_Temp <- KildareCrime2011
KildareCrime2011 <- spTransform(KildareCrime2011_Temp, newProj)

Tidy up the column names before beginning to plot the choropleth maps.

# Tidy up the column names
names(KildareCrime2011@data)[6] <- paste("Attempts_Murder")
names(KildareCrime2011@data)[7] <- paste("Dangerous_Acts")
names(KildareCrime2011@data)[8] <- paste("Kidnapping")
names(KildareCrime2011@data)[9] <- paste("Robbery")
names(KildareCrime2011@data)[10] <- paste("Burglary")
names(KildareCrime2011@data)[11] <- paste("Theft")
names(KildareCrime2011@data)[12] <- paste("Fraud")
names(KildareCrime2011@data)[13] <- paste("Controlled_Drugs")
names(KildareCrime2011@data)[14] <- paste("Weapons")
names(KildareCrime2011@data)[15] <- paste("Damage")
names(KildareCrime2011@data)[16] <- paste("Public_Order")
names(KildareCrime2011@data)[17] <- paste("Offences_Gov")
Step 2)

Finally, before plotting any maps - reproject KildareGSD and create labels for each GSD to use when each individual map is plotted.

# Project Kildare GSD correctly
KildareGSD_Temp <- KildareGSD
KildareGSD <- spTransform(KildareGSD_Temp, newProj)
# Create Label names for each GSD. First, get centroid for each GSD
GSDcentroid <- gCentroid(KildareCrime2011, byid = TRUE)
LabName <- KildareCrime2011@data$Station
Lat <- GSDcentroid$x
Lon <- GSDcentroid$y

Plot Population Distribution

Population distribution at GSD level

# Choropleth map of population of Kildare at GSD level
par(mar = c(1, 0, 5, 0))
shades2 <- auto.shading(KildareGSD$Total2011, n = 6, cols = brewer.pal(6, "Purples"), 
    cutter = rangeCuts)
choropleth(KildareGSD, KildareGSD$Total2011, shading = shades2)
pl <- pointLabel(Lat, Lon, LabName, cex = 0.7)
plot(GSDcentroid, pch = 16, cex = 0.5, add = T)
choro.legend(694425.2, 704194.4, shades2, cex = 0.7, title = "Population")
map.scale(701637.9, 682556.3, 10000, "Kilometers", 2, 5)
north.arrow(705244.3, 713874.5, km2ft(0.25), col = "darkgray", cex = 0.7)
title(main = "Population of Kildare\n at Garda Sub-District Level", font.main = 6)

Populaiton distribution at SA level

# Plot a standard choropleth of populaiton of Kildare per SA
par(mar = c(1, 0, 5, 0))
shades5 <- auto.shading(KildareSA$TOTAL2011, n = 6, cols = brewer.pal(6, "Purples"), 
    cutter = rangeCuts)
choropleth(KildareSA, KildareSA$TOTAL2011, shading = shades5)
choro.legend(694425.2, 704194.4, shades5, cex = 0.7, title = "Population")
map.scale(701637.9, 682556.3, 10000, "Kilometers", 2, 5)
north.arrow(705244.3, 713874.5, km2ft(0.25), col = "darkgray", cex = 0.7)
title(main = "Population of Kildare,\n at Small Area Level", font.main = 6)

Population distribution at SA level (Cartogram)

library(Rcartogram)
library(getcartr)
# Plot a choropleth cartogram of Kildare SAs - Total Population from 2011
# Census
Kildare.carto <- quick.carto(KildareSA, KildareSA$TOTAL2011)
par(mfrow = c(1, 1))
par(mar = c(0, 0, 5, 0))
shades4 <- auto.shading(Kildare.carto$TOTAL2011, n = 6, cols = brewer.pal(6, 
    "Purples"), cutter = rangeCuts)
choropleth(Kildare.carto, Kildare.carto$TOTAL2011, shading = shades4)
choro.legend(93.34762, 42.83977, shades4, cex = 0.7, title = "Population")
north.arrow(112.4466, 58.34041, km2ft(5e-04), col = "darkgray", cex = 0.7)
title(main = "Population of Kildare,\n at Small Area Level\n(Cartogram)", font.main = 6)

Plot Crime Distribution

Crime Category One

Attempts or threats to murder, assaults, harrassments and related offences (2011)

par(mar = c(1, 0, 5, 0))
# Set up shading scheme
shades = auto.shading(KildareCrime2011$Attempts_Murder, n = 6, cols = brewer.pal(6, 
    "Blues"), cutter = rangeCuts)
# Plot the choropleth
choropleth(KildareCrime2011, KildareCrime2011$Attempts_Murder, shading = shades)
pl <- pointLabel(Lat, Lon, LabName, cex = 0.7)
plot(GSDcentroid, pch = 16, cex = 0.5, add = T)
choro.legend(694425.2, 704194.4, shades, cex = 0.7, title = "Recorded Crime")
map.scale(701637.9, 682556.3, 10000, "Kilometers", 2, 5)
north.arrow(705244.3, 713874.5, km2ft(0.25), col = "darkgray", cex = 0.7)
title(main = "Crime Category 1,\nat Garda Sub-District level in Co. Kildare", 
    font.main = 6)

Crime Category Two

Dangerous or negligent acts (2011)

par(mar = c(1, 0, 5, 0))
# Set up shading scheme
shades = auto.shading(KildareCrime2011$Dangerous_Acts, n = 6, cols = brewer.pal(6, 
    "Blues"), cutter = rangeCuts)
# Plot the choropleth
choropleth(KildareCrime2011, KildareCrime2011$Dangerous_Acts, shading = shades)
pl <- pointLabel(Lat, Lon, LabName, cex = 0.7)
plot(GSDcentroid, pch = 16, cex = 0.5, add = T)
choro.legend(694425.2, 704194.4, shades, cex = 0.7, title = "Recorded Crime")
map.scale(701637.9, 682556.3, 10000, "Kilometers", 2, 5)
north.arrow(705244.3, 713874.5, km2ft(0.25), col = "darkgray", cex = 0.7)
title(main = "Crime Category 2,\nat Garda Sub-District level in Co. Kildare", 
    font.main = 6)

Crime Category Three

Kidnapping and related offences (2011)

par(mar = c(0, 0, 5, 0))
# Set up shading scheme
shades = auto.shading(KildareCrime2011$Kidnapping, n = 6, cols = brewer.pal(6, 
    "Blues"), cutter = rangeCuts)
# Plot the choropleth
choropleth(KildareCrime2011, KildareCrime2011$Kidnapping, shading = shades)
pl <- pointLabel(Lat, Lon, LabName, cex = 0.7)
plot(GSDcentroid, pch = 16, cex = 0.5, add = T)
choro.legend(694425.2, 704194.4, shades, cex = 0.7, title = "Recorded Crime")
map.scale(701637.9, 682556.3, 10000, "Kilometers", 2, 5)
north.arrow(705244.3, 713874.5, km2ft(0.25), col = "darkgray", cex = 0.7)
title(main = "Crime category 3,\nat Garda Sub-District level in Co. Kildare", 
    font.main = 6)

Crime Category Four

Robbery, extortion and hijacking offences (2011)

par(mar = c(0, 0, 5, 0))
# Set up shading scheme
shades = auto.shading(KildareCrime2011$Robbery, n = 6, cols = brewer.pal(6, 
    "Blues"), cutter = rangeCuts)
# Plot the choropleth
choropleth(KildareCrime2011, KildareCrime2011$Robbery, shading = shades)
pl <- pointLabel(Lat, Lon, LabName, cex = 0.7)
plot(GSDcentroid, pch = 16, cex = 0.5, add = T)
choro.legend(694425.2, 704194.4, shades, cex = 0.7, title = "Recorded Crime")
map.scale(701637.9, 682556.3, 10000, "Kilometers", 2, 5)
north.arrow(705244.3, 713874.5, km2ft(0.25), col = "darkgray", cex = 0.7)
title(main = "Crime category 4,\nat Garda Sub-District level in Co. Kildare", 
    font.main = 6)

Crime Category Five

Burglary and related offences (2011)

par(mar = c(0, 0, 5, 0))
# Set up shading scheme
shades = auto.shading(KildareCrime2011$Burglary, n = 6, cols = brewer.pal(6, 
    "Blues"), cutter = rangeCuts)
# Plot the choropleth
choropleth(KildareCrime2011, KildareCrime2011$Burglary, shading = shades)
pl <- pointLabel(Lat, Lon, LabName, cex = 0.7)
plot(GSDcentroid, pch = 16, cex = 0.5, add = T)
choro.legend(694425.2, 704194.4, shades, cex = 0.7, title = "Recorded Crime")
map.scale(701637.9, 682556.3, 10000, "Kilometers", 2, 5)
north.arrow(705244.3, 713874.5, km2ft(0.25), col = "darkgray", cex = 0.7)
title(main = "Crime category 5,\nat Garda Sub-District level in Co. Kildare", 
    font.main = 6)

Crime Category Six

Theft and related offences (2011)

par(mar = c(0, 0, 5, 0))
# Set up shading scheme
shades = auto.shading(KildareCrime2011$Theft, n = 6, cols = brewer.pal(6, "Blues"), 
    cutter = rangeCuts)
# Plot the choropleth
choropleth(KildareCrime2011, KildareCrime2011$Theft, shading = shades)
pl <- pointLabel(Lat, Lon, LabName, cex = 0.7)
plot(GSDcentroid, pch = 16, cex = 0.5, add = T)
choro.legend(694425.2, 704194.4, shades, cex = 0.7, title = "Recorded Crime")
map.scale(701637.9, 682556.3, 10000, "Kilometers", 2, 5)
north.arrow(705244.3, 713874.5, km2ft(0.25), col = "darkgray", cex = 0.7)
title(main = "Crime category 6,\nat Garda Sub-District level in Co. Kildare", 
    font.main = 6)

Crime Category Seven

Fraud, deception and related offences (2011)

par(mar = c(0, 0, 5, 0))
# Set up shading scheme
shades = auto.shading(KildareCrime2011$Fraud, n = 6, cols = brewer.pal(6, "Blues"), 
    cutter = rangeCuts)
# Plot the choropleth
choropleth(KildareCrime2011, KildareCrime2011$Fraud, shading = shades)
pl <- pointLabel(Lat, Lon, LabName, cex = 0.7)
plot(GSDcentroid, pch = 16, cex = 0.5, add = T)
choro.legend(694425.2, 704194.4, shades, cex = 0.7, title = "Recorded Crime")
map.scale(701637.9, 682556.3, 10000, "Kilometers", 2, 5)
north.arrow(705244.3, 713874.5, km2ft(0.25), col = "darkgray", cex = 0.7)
title(main = "Crime category 7,\nat Garda Sub-District level in Co. Kildare", 
    font.main = 6)

Crime Category Eight

Controlled drugs offences (2011)

par(mar = c(0, 0, 5, 0))
# Set up shading scheme
shades = auto.shading(KildareCrime2011$Controlled_Drugs, n = 6, cols = brewer.pal(6, 
    "Blues"), cutter = rangeCuts)
# Plot the choropleth
choropleth(KildareCrime2011, KildareCrime2011$Controlled_Drugs, shading = shades)
pl <- pointLabel(Lat, Lon, LabName, cex = 0.7)
plot(GSDcentroid, pch = 16, cex = 0.5, add = T)
choro.legend(694425.2, 704194.4, shades, cex = 0.7, title = "Recorded Crime")
map.scale(701637.9, 682556.3, 10000, "Kilometers", 2, 5)
north.arrow(705244.3, 713874.5, km2ft(0.25), col = "darkgray", cex = 0.7)
title(main = "Crime category 8,\nat Garda Sub-District level in Co. Kildare", 
    font.main = 6)

Crime Category Nine

Weapons and explosive offences (2011)

par(mar = c(0, 0, 5, 0))
# Set up shading scheme
shades = auto.shading(KildareCrime2011$Weapons, n = 6, cols = brewer.pal(6, 
    "Blues"), cutter = rangeCuts)
# Plot the choropleth
choropleth(KildareCrime2011, KildareCrime2011$Weapons, shading = shades)
pl <- pointLabel(Lat, Lon, LabName, cex = 0.7)
plot(GSDcentroid, pch = 16, cex = 0.5, add = T)
choro.legend(694425.2, 704194.4, shades, cex = 0.7, title = "Recorded Crime")
map.scale(701637.9, 682556.3, 10000, "Kilometers", 2, 5)
north.arrow(705244.3, 713874.5, km2ft(0.25), col = "darkgray", cex = 0.7)
title(main = "Crime category 9,\nat Garda Sub-District level in Co. Kildare", 
    font.main = 6)

Crime Category Ten

Damage to the property and to the environment (2011)

par(mar = c(0, 0, 5, 0))
# Set up shading scheme
shades = auto.shading(KildareCrime2011$Damage, n = 6, cols = brewer.pal(6, "Blues"), 
    cutter = rangeCuts)
# Plot the choropleth
choropleth(KildareCrime2011, KildareCrime2011$Damage, shading = shades)
pl <- pointLabel(Lat, Lon, LabName, cex = 0.7)
plot(GSDcentroid, pch = 16, cex = 0.5, add = T)
choro.legend(694425.2, 704194.4, shades, cex = 0.7, title = "Recorded Crime")
map.scale(701637.9, 682556.3, 10000, "Kilometers", 2, 5)
north.arrow(705244.3, 713874.5, km2ft(0.25), col = "darkgray", cex = 0.7)
title(main = "Crime category 10,\nat Garda Sub-District level in Co. Kildare", 
    font.main = 6)

Crime Category Eleven

Public order and other social code offences (2011)

par(mar = c(0, 0, 5, 0))
# Set up shading scheme
shades = auto.shading(KildareCrime2011$Public_Order, n = 6, cols = brewer.pal(6, 
    "Blues"), cutter = rangeCuts)
# Plot the choropleth
choropleth(KildareCrime2011, KildareCrime2011$Public_Order, shading = shades)
pl <- pointLabel(Lat, Lon, LabName, cex = 0.7)
plot(GSDcentroid, pch = 16, cex = 0.5, add = T)
choro.legend(694425.2, 704194.4, shades, cex = 0.7, title = "Recorded Crime")
map.scale(701637.9, 682556.3, 10000, "Kilometers", 2, 5)
north.arrow(705244.3, 713874.5, km2ft(0.25), col = "darkgray", cex = 0.7)
title(main = "Crime category 11,\nat Garda Sub-District level in Co. Kildare", 
    font.main = 6)

Crime Category Twelve

Offences against government, justice procedures and organisation of crime (2011)

par(mar = c(0, 0, 5, 0))
# Set up shading scheme
shades = auto.shading(KildareCrime2011$Offences_Gov, n = 6, cols = brewer.pal(6, 
    "Blues"), cutter = rangeCuts)
# Plot the choropleth
choropleth(KildareCrime2011, KildareCrime2011$Offences_Gov, shading = shades)
pl <- pointLabel(Lat, Lon, LabName, cex = 0.7)
plot(GSDcentroid, pch = 16, cex = 0.5, add = T)
choro.legend(694425.2, 704194.4, shades, cex = 0.7, title = "Recorded Crime")
map.scale(701637.9, 682556.3, 10000, "Kilometers", 2, 5)
north.arrow(705244.3, 713874.5, km2ft(0.25), col = "darkgray", cex = 0.7)
title(main = "Crime category 12,\nat Garda Sub-District level in Co. Kildare", 
    font.main = 6)

References

An Garda Síochána (2016) Station Directory. [online]
Available at: http://www.garda.ie/Stations/Map.aspx?gardaKMZFilename=Garda_Stations.kmz&Mode=gardaViewAll&Args=undefined (accessed 01 May 2016).

Brunsdon, C. (2015a) Package ‘getCartr’ [online] Avaliable at: https://github.com/chrisbrunsdon/getcartr (accessed 1 August 2016).

Brunsdon, C. (2015b) Package ‘tbart’. [online]
Available at: https://cran.r-project.org/web/packages/tbart/tbart.pdf (accessed 1 Jube 2016).

Brunsdon, C. and Singleton, A. (2015) Reproducible Research: Concepts, Techniques and Issues. In: Brunsdon, C. and Singleton, A., eds. Geocomputation: A Practical Primer. London: SAGE Publications Ltd, pp. 254 - 263.

Central Statistics Office (2011) 2011 Census Boundaries. [online]
Available at: http://census.cso.ie/censusasp/saps/boundaries/ED_SA%20Disclaimer1.htm (accessed 1 May 2016).

Geofabrik (2016) Ireland and Northern Ireland. [online]
Available at: http://download.geofabrik.de/europe/ireland-and-northern-ireland.html (accessed 1 May 2016).

Lu, B. (2014) Package ‘shp2graph’. [online]
Avaiable at: https://cran.r-project.org/web/packages/shp2graph/shp2graph.pdf (accessed 01 July 2016).

M.Sc. Geocomputation