It is important that you keep all these files in the same location as each other! They all contain different bits of information about your shapefile (and they are all needed). Below is a list of the files you might see in this collection:
Section 1.10 Mapping crime data as polygons
What about our other two columns, location, and LSOAs? Well, to put these on the map, we need a geometry representation of them. We need boundary data representing the areas we want to map. We will learn in this section where you may find, download and turn them into sf objects, and how to link our dataframe as attribute data in order to be able to map them.
Subsection 1.10.1 Finding boundary data
In this section you are going to learn how you take one of the most popular data formats for spatial objects, the shapefile, and read it into R. The shapefile was introduced by ESRI, the developers and vendors of ArcGIS. And although many other formats have developed since and ESRI no longer holds the same market position it once occupied (though they’re still the player to beat), shapefiles continue to be one of the most popular formats you will encounter in your work.
We are going to obtain shapefiles for British census geographies. For this activity, we will focus on the polygon (LSOA) rather than the lines of the streets, but the logic is more or less the same.
Census boundary data are a digitised representation of the underlying geography of the census. Census Geography is often used in research and spatial analysis because it is divided into units based on population counts, created to form comparable units, rather than other administrative boundaries such as wards or police force areas. However, depending on your research question and the context for your analysis, you might be using different units.
The hierarchy of the census geographies in the UK goes from Country to Local Authority to Middle Layer Super Output Area (MSOA) to Lower Layer Super Output Area (LSOA) to Output Area (in other countries you have similar levels in the census hierarchies):
country <- rnaturalearth::ne_countries(scale = 110, geounit = "united kingdom", returnclass = "sf")
la <- rnaturalearth::ne_states(geounit = "england", returnclass = "sf") %>% filter(region == "North West")
lsoas <- st_read("data/BoundaryData/england_lsoa_2011.shp") %>% st_transform(st_crs(la))
p1 <- ggplot() +
geom_sf(data = lsoas, aes(col = "LSOAs")) +
geom_sf(data = country, aes(col = "country")) +
geom_sf(data = la, aes(col = "local authority")) +
theme_void() +
scale_colour_manual(values = c("country" = "#1b9e77",
"local authority" = "#d95f02",
"LSOAs" = "#7570b3"),
name = "Boundary: ")
p2 <- ggplot() +
geom_sf(data = la , aes(col = "local authority")) +
geom_sf(data = lsoas, aes(col = "lsoas")) + theme_void() +
scale_colour_manual(values = c("country" = "#1b9e77",
"local authority" = "#d95f02",
"LSOAs" = "#7570b3"),
name = "Boundary: ")
p3 <- ggplot() +
# geom_sf(data = la %>% filter(name == "Manchester"), aes(col = "local authority")) +
geom_sf(data = lsoas, aes(col = "LSOAs")) + theme_void() +
scale_colour_manual(values = c("country" = "#1b9e77",
"local authority" = "#d95f02",
"LSOAs" = "#7570b3"),
name = "Boundary: ")
ggpubr::ggarrange(p1, p2, p3, nrow = 1, common.legend = TRUE)
Here we will get some boundaries for Manchester. We know that our crime data has a column for LSOA, so we can proceed with this as our appropriate unit of analysis.
To get some boundary data, you will need to be able to source geographic data for your study area. In this case, since we are working with data in the UK, we can use the UK Data Service website to find this. There is a simple Boundary Data Selector tool which you could use. Other countries also have such repositories, and in some cases where they do not, other resources, such as Open Street Map, or natural earth (and associated R package
rnaturalearth) can be handy resources. In Appendix C: Sourcing geographical data for crime analysis, we elaborate more on these, with some examples.
For now, you can turn to the data folder which you have downloaded from this book’s repository. (If unsure, refer to the Preface). Specifically, we’re looking for a folder called
BoundaryData within the "data" folder. As mentioned, we will be working with a shapefile, but in the case of this format we’re actually working with many files. When we download a shapefile, we usually get many files, all bundled inside a folder which contains them all. It is important for these files to be kept together, in the same folder. If you have a look inside this BoundaryData folder, you will notice that there are 4 files with the same name "england_lsoa_2011", but with different extensions.
Warning 1.10.2.
-
.shp — shape format; the feature geometry itself, this is what you see on the map
-
.shx — shape index format; a positional index of the feature geometry to allow seeking forwards and backwards quickly
-
.dbf — attribute format; columnar attributes for each shape, in dBase IV format.
-
.prj — projection format; the coordinate system and projection information, a plain text file describing the projection using well-known text format
Sometimes there might be more files associated with your shapefile as well, but we will not cover them here. So unlike when you work with spreadsheets and data in tabular form, which typically is just all included in one file, when you work with shapefiles, you have to live with the required information living in separate files that need to be stored together. So, being tidy and organised is even more important when you carry out projects that involve spatial data.
Subsection 1.10.2 Reading shapefiles into R
To read in your data into R, you will need to know the path to where you have saved it. Ideally this will be in your data folder in your project directory.
Let’s create an object and assign it our shapefile’s name:
# Remember to use the appropriate pathfile in your case
shp_name <- "data/BoundaryData/england_lsoa_2011.shp"
Make sure that this is saved in your working directory, and you have set your working directory. Now use the
st_read() function from the sf package to read in the shapefile (the one with the .shp extension):
manchester_lsoa <- st_read(shp_name)
Notice that you don’t have to read in any of the other files (e.g.,
.proj or .dbf); it is enough that they are there silently in the folder from which you read in the .shp extension file. Now you have your spatial data file. Notice how running the function sends to the console some metadata about your data. You have a polygon, with 282 rows, and the CRS is the projected British National Grid. You can have a look at what sort of data it contains, the same way you would view a dataframe, with the View() function:
View(manchester_lsoa)
And of course, since it’s spatial data, you can map it using the
geom_sf() function, as we did with our points:
ggplot() +
geom_sf(data = manchester_lsoa)
Great, we now have an outline of the LSOAs in Manchester. Notice how the shape is different to that of the points in our crime data, since here we only obtained the data for the city of Manchester (Manchester Local Authority) rather than for the whole metropolitan area, Greater Manchester which includes all the other local authorities aside from Manchester city as well.
Subsection 1.10.3 Data wrangling with dplyr
In order to map crimes to LSOAs, we might want to take a step back and think about the unit of analysis at which our data are collected. In our original dataframe of crimes, we saw that each crime incident is one row. So the unit of analysis is each crime. Since we were looking to map each crime at the location it happened, we used the latitude and longitude supplied for each incident, providing a geometry for mapping each crime incident. However, when we are looking to map our data to LSOA level, we need to match the crime data to the geometry we wish to display.
Have a look at the "manchester_lsoa" object we mapped above. How many rows (observations) does it have? You can check this by looking in the Environment pane, or by using the
nrow() function.
nrow(manchester_lsoa)
You can see this has (number of Manchester LSOAs) rows. This means we have geometries for (number of Manchester LSOAs) LSOAs. On the other hand, our crimes dataframe has (number of crimes) rows, one for each crime (observation). So how can we match these up? The answer lies in thinking about what it is that our map using LSOAs as our unit of analysis will be able to tell us. Think of other maps of areas. What are they usually telling you? Usually we expect to see crimes per neighbourhood, something like this. So our unit of analysis needs to be LSOA, and for each one we need to know how many crimes occurred in that area.
To achieve this, we will wrangle our data using functions from the
dplyr package. This is a package for conducting all sorts of operations with dataframes. We are not going to cover the full functionality of dplyr (which you can consult in the dplyr vignette [240]), but we are going to cover three different very useful elements of dplyr: the select function, the group_by function, and the piping operator.
The
select() function provides you with a simple way of subsetting columns from a dataframe. So, say, we just want to use one variable, "lsoa_code", from the "crimes" dataframe and store it in a new object we could write the following code. This variable tells us the LSOA in which the crime took place and it is essential if we want to group our crimes by LSOA (by using this standard method of merging information from datasets):
new_object <- select(crimes, lsoa_code)
We can also use the
group_by() function for performing group operations. Essentially this function asks R to group cases within categories and then do something with those grouped cases. So, say, we want to count the number of cases within each LSOA, we could use the following code:
#First we group the cases by LSOA code and
# store this organised data into a new object
grouped_crimes <- group_by(new_object, lsoa_code)
#Then we count the number of cases within each category
# using the summarise function to print the results
summarise(grouped_crimes, count = n())
#We create a new dataframe with these results
crime_per_LSOA <- summarise(grouped_crimes, count = n())
As you can see, we can do what we wanted, create a new dataframe with the required info. But if we do this, we are creating many objects that we don’t need, one at each step. Instead there is a more efficient way of doing this, without so many intermediate steps clogging up our environment with unnecessary objects. That’s where the piping operator comes handy. The piping operator is written like
%>% and it can be read as "and then". Look at the code below:
#First we say create a new object called crime_per_lsoa,
# and then select only the LSOA.code column to exist in this object,
# and then group this object by the LSOA.code,
# and then count the number of cases within each category.
# This is what I want in the new object.
crimes_per_lsoa <- crimes %>%
group_by(lsoa_code) %>%
summarise(count=n())
Essentially we obtain the same results but with more streamlined and elegant code, and not needing additional objects in our environment. As of version 4.1.0 of R, there is also a new piping operator denoted by the
|> symbol that can also be used and that does not depend on the magrittr package (as %>% does).
And now we have a new object, "crimes_per_lsoa". If we have a look at this one, we can now see that each row represents one LSOA, and next to it we have a variable for the number of crimes from each area. We created a new dataframe from a frequency table, and as each row of the crimes data was one crime, the frequency table tells us the number of crimes which occurred in each LSOA.
Those of you playing close attention might note that there are still more observations in this dataframe ((number of LSOAs with crimes)) than in the "manchester_lsoas" one ((number of Manchester LSOAs)). Again, this is because "crimes_per_lsoa" also includes data from census areas in municipalities of the metropolitan area of Greater of Manchester other than Manchester city local authority.
Subsection 1.10.4 Join data to sf object
Our next task is to link our crimes data to our
sf spatial object to help us map this. Notice anything similar between the data from the shapefile and the frequency table data we just created? Do they share a column? Indeed, you might notice that the "lsoa_code" field in the crimes data matches the values in the "code" field in the spatial data. In theory we could join these two data tables.
So how do we do this? Well what you can do is to link one data set with another. Data linking is used to bring together information from different sources in order to create a new, richer dataset. This involves identifying and combining information from corresponding records on each of the different source datasets. The records in the resulting linked dataset contain some data from each of the source datasets. Most linking techniques combine records from different datasets if they refer to the same entity (an entity may be a person, organisation, household or even a geographic region.)
You can merge (combine) rows from one table into another just by pasting them in the first empty cells below the target table — the table grows in size to include the new rows. And if the rows in both tables match up, you can merge columns from one table with another by pasting them in the first empty cells to the right of the table — again, the table grows, this time to include the new columns.
Merging rows is pretty straightforward, but merging columns can be tricky if the rows of one table don’t always line up with the rows in the other table. By using
left_join() from the dplyr package, you can avoid some of the alignment problems. The left_join() function will return all rows from \(x\text{,}\) and all columns from \(x\) and \(y\text{.}\) Rows in \(x\) with no match in \(y\) will have NA values in the new columns. If there are multiple matches between \(x\) and \(y\text{,}\) all combinations of the matches are returned.
So we’ve already identified that both our crimes data, and the spatial data contain a column with matching values, the codes for the LSOA that each row represents. You need a unique identifier to be present for each row in all the data sets that you wish to join. This is how R knows what values belong to what row. What you are doing is matching each value from one table to the next, using this unique identified column, that exists in both tables. For example, let’s say we have two data sets from some people in Hawkins, Indiana. In one data set we collected information about their age. In another one, we collected information about their hair colour. If we collected some information that is unique to each observation, and this is the same in both sets of data, for example their names, then we can link them up, based on this information. Something like this:
library(grid)
library(gridExtra)
t1 <- data.frame(name = c("Barb", "Steve", "Mike"),
age = c(16, 16, 13))
t2 <- data.frame(name = c("Barb", "Steve", "Mike"),
hair_colour = c("Red", "Brown", "Black"))
grid.arrange(tableGrob(t1), tableGrob(t2), nrow = 1)
And by doing so, we produce a final table that contains all values, lined up correctly for each individual observation, like this:
t3 <- left_join(t1, t2)
grid.arrange(tableGrob(t3))
This is all we are doing, when merging tables: joining the two data sets using a common column. Why are we using left join though? There is a whole family of join functions as part of dplyr which join datasets. There is also a
right_join(), and an inner_join() and an outer_join() and a full_join(). But here we use left_join(), because that way we keep all the rows in \(x\) (the left-hand side dataframe), and join to it all the matched columns in \(y\) (the right-hand side dataframe).
So let’s join the crimes data to the spatial data, using
left_join(). We have to tell the function what are the dataframes we want to join, as well as the names of the columns that contain the matching values in each one. This is "code" in the "manchester_lsoa" dataframe and "lsoa_code" in the "crimes_per_lsoa" dataframe. Like so:
manchester_lsoa <- left_join(manchester_lsoa, crimes_per_lsoa,
by = c("code"="lsoa_code"))
Now if you have a look at the data again, you will see that the column of number of crimes (
count) has been added on. You may not want to have to go through this process all the time that you want to work with this data. One thing you could do is to save the "manchester_lsoa" object as a physical file in your machine. You can use the st_write() function from the sf package to do this. If we want to write into a shapefile format, we would do as shown below. Make sure you save this file, for we will come back to it in subsequent chapters.
st_write(manchester_lsoa, "data/BoundaryData/manchester_crime_lsoa.shp")
Remark 1.10.6.
Note: if you get an error which says the layer already exists, this means there is already a file with this name in that location. You can add the parameter
append = FALSE to overwrite this file, or give it a slightly different name.
Subsection 1.10.5 Putting polygon data on the map
Now that we have joined the crimes data to the geometry, you can use this to make our map. Remember our original empty map? Well now, since we have the column (variable) for number of crimes here, we can use that to shade the polygons based on how many crimes there are in each LSOA. We can do this by specifying the
fill= parameter of the geom_sf function.
ggplot() +
geom_sf(data = manchester_lsoa, aes(fill = count))
Just like we did with the ggplot bar charts we created earlier, we can adjust the opacity of our thematic map, with the
alpha = parameter, add a basemap, with the function annotation_map_tile() and adjust the colour scheme, with the scale_fill_gradient2() function.
ggplot() +
annotation_map_tile() + # add basemap
geom_sf(data = manchester_lsoa,
aes(fill = count),
alpha = 0.7) + # alpha sets the opacity
#use scale_fill_gradient2() for a different palette
# and name the variable on the legend
scale_fill_gradient2(name ="Number of crimes")
In subsequent chapters we will play around with other packages in R that you can use to produce this kind of maps. Gradually we will discuss the kind of choices you can make in order to select the adequate type of representation for the type of spatial data and question you have, and the aesthethic choices that are adequate depending on the purposes and medium in which you will publish your map. But for now, we can rejoice: here is our very first crime map of the book!
