Section 2.4 Spatial operations
Spatial operations are a vital part of geocomputation. Spatial objects can be modified in a multitude of ways based on their location and shape. For a comprehensive overview of spatial operations in R, we recommend chapter 4 of [115].
Spatial operations differ from non-spatial operations in some ways. To illustrate the point, imagine you are researching road safety. Spatial joins can be used to find road speed limits related with administrative zones, even when no zone ID is provided. But this raises the question: should the road completely fall inside a zone for its values to be joined? Or is simply crossing or being within a certain distance sufficient? When posing such questions, it becomes apparent that spatial operations differ substantially from attribute operations on dataframes: the type of spatial relationship between objects must be considered. ([115])
So you can see, we can do exciting spatial operations with our spatial data, which we cannot with the non-spatial stuff.
Subsection 2.4.1 Reprojecting coordinates
It is important to recall here some of the learning from the previous chapter on map projections and coordinate reference systems. We learned about ways of flattening out the earth, and ways of making sense of what that means for how to be able to point to specific locations in our maps. Coordinate Reference System (CRS) is this method of how to refer to locations with our data. You might use a Geographic Coordinate System, which tells you where your data are located on the surface of the Earth. The most commonly used one is the WGS 84, where we define our locations with latitude and longitude points. The other type is a Projected Coordinate System which tells the data how to draw on a flat, two-dimensional surface (such as a computer screen). In our case here, we will often encounter the British National Grid when working with British data. Here our locations are defined with Eastings and Northings.
So why are we talking about this?
It is important to note that spatial operations that use two spatial objects rely on both objects having the same coordinate reference system.
If we are looking to carry out operations that involve two different spatial objects, they need to have the same CRS! Funky weird things happen when this condition is not met, so beware! So how do we know what CRS our spatial objects are? Well the
sf package contains a handy function called st_crs() which let us check. All you need to pass into the brackets of this function is the name of the object you want to know the CRS of.
Let’s check what is the CRS of our crimes:
st_crs(crimes)
## Coordinate Reference System: NA
You can see that we get the CRS returned as
NA. Can you think why? Have we made this into a spatial object? Or is this merely a dataframe with a latitude and longitude column? The answer is really in the question here. So we need to convert this to a sf object, or a spatial object, and make sure that R knows that the latitude and the longitude columns are, in fact, coordinates.
crimes <- st_as_sf(crimes, coords = c("longitude", "latitude"),
crs = 4326, agr = "constant", na.fail = FALSE)
In the
st_as_sf() function, we specify what we are transforming (the name of our dataframe), the column names that have the coordinates in them (longitude and latitude), the CRS we are using (4326 is the code for WGS 84, which is the CRS that uses latitude and longitude coordinates (remember BNG uses Easting and Northing)), and finally agr, the attribute-geometry-relationship, specifies for each non-geometry attribute column how it relates to the geometry, and can have one of following values: "constant", "aggregate", "identity". The option "constant" is used for attributes that are constant throughout the geometry (e.g., land use), "aggregate" where the attribute is an aggregate value over the geometry (e.g., population density or population count), "identity" when the attributes uniquely identifies the geometry of particular "thing", such as a building ID or a city name. The default value, NA_agr_, implies we don’t know. Finally we can specify the parameter na.fail which decides what to do when one or more of our observations have missing (NA) values in their coordinates. The default setting is for the entire process to fail in this case. However, if we set this to FALSE it simply means we will lose those observations. Essentially, we are getting rid of those violent crimes which have not been geocoding, excluding them from our analysis.
Now let’s check the CRS of this spatial version of our licensed premises:
st_crs(crimes)
We can now see that we have this coordinate system as WGS84. We need to then make sure that any other spatial object with which we want to perform spatial operations is also in the same CRS. Let’s look at our city centre ward boundary file:
st_crs(city_centre)
We see that this is in fact in a projected coordinate system, namely the British National Grid we mentioned. To make them align, we can re-project this object into the WGS84 geographic coordinate system. To do this, we can use the
st_transform() function.
city_centre <- st_transform(city_centre, crs = 4326)
Now we can check the projection again:
st_crs(city_centre)
And we can also check whether the CRS of the two objects match:
st_crs(crimes) == st_crs(city_centre)
## [1] TRUE
It is true! Finally, to check our bar data from Open Street Map:
st_crs(osm_bar_sf)
You will see this is also in WGS84. Since all our data are in this CRS, we can now move on to carry out some spatial operations.
But before doing so, we need to raise an important point. R is still undergoing the adaptation to new standards in the way CRS data is stored and handled. Before 2020 you could specify the CRS data with a EPSG code or with a
proj4string reference, which look more like gibberish to the non-initiated but provided more flexibility (see [115] for details). But there were problems with the old PROJ4. As we noted in the first chapter, it only allows for static reference frames (when on Earth everything is moving all the time). So PROJ the main library for performing conversions between systems has changed. It now discourages the old proj4string format and has changed to the OGC WKT2. When you consult the https://epsg.io/ site, you can see whenever you select a projection it translates to different standards (including OGC WKT). This is one of those changes that long-term are really good. But because the cross-dependencies between spatial R packages were being updated at different speed, this had an impact on usability. Although we will stick to the more familiar EPSG notation, we strongly recommend reading the following technical notes by Pebesma and Bivand ([116], [117], [118]), and to be particularly alert to this issue.
Subsection 2.4.2 Subsetting points
Recall above that we wanted to focus our efforts on the City Centre ward of Manchester; however, for our bounding box to download OSM data we used Greater Manchester. If we were to plot our bars, we would see that we have many which fall outside of the City Centre ward:
plot(st_geometry(osm_bar_sf), col = 'red')
plot(st_geometry(city_centre), add = TRUE)
This is also the case for our crimes data:
plot(st_geometry(osm_bar_sf), col = 'red')
plot(st_geometry(crimes), col = 'blue', add = TRUE)
plot(st_geometry(city_centre), add = TRUE)
So if we really want to focus on City Centre, we should create spatial objects for the crimes and the bars which include only those which fall within the City Centre ward boundary.
First, we check whether they have the same CRS.
st_crs(city_centre) == st_crs(crimes)
## [1] TRUE
We do indeed, as we made sure in the previous section. Now we can move on to our spatial operation, where we select only those points within the City Centre polygon. To do this, we first make a list of intersecting points to the polygon, using the
st_intersects() function. This function takes two arguments: first, the polygon which we want to subset our points within, and second, the points which we want to subset. We then use the resulting "cc_crimes" object to subset the crimes object to include only those which intersect (return TRUE for intersects):
# intersection
cc_crimes <- st_intersects(city_centre, crimes)
# subsetting
cc_crimes <- crimes[unlist(cc_crimes),]
Have a look at this new "cc_crimes" object in your environment. How many observations does it have? Is this now fewer than the previous "crimes" object? Why do you think this is?
(Hint: you’re removing everything that is outside the City Centre polygon)
We can plot this again to have a look:
plot(st_geometry(city_centre))
plot(st_geometry(cc_crimes), col = 'blue', add = TRUE)
We have successfully performed our first spatial operation. We managed to subset our points dataset of crimes to include only those crimes which are located inside the polygon for City Centre.
We can do the same for the bars:
# intersection
cc_bars <- st_intersects(city_centre, osm_bar_sf)
# subsetting
cc_bars <- osm_bar_sf[unlist(cc_bars),]
We can see that of the previous bars, some number are within the City Centre ward. We can plot our data now:
plot(st_geometry(city_centre))
plot(st_geometry(cc_bars), col = 'red', add = TRUE)
plot(st_geometry(cc_crimes), col = 'blue', add = TRUE)
Subsection 2.4.3 Building buffers
So we now have our bars and our violent crimes in Manchester City Centre. Let’s go back to our original question. We want to know about crime in and around our areas of interest, in this case our bars. But how can we count this? We have our points that are crimes, right? Well... How do we connect them to our points that are licensed premises?
One approach is to build a buffer around our bars, and say that we will count all the crimes which fall within a specific radius of this bar. What should this radius be? Well, this is where your domain knowledge as criminologist or crime analyst comes in. How far away would you consider a crime to still be related to this pub? 400 metres? 500 metres? 900 metres? 1 km? What do you think? This is one of those it depends questions, where there is no universal right answer; instead it will depend on the environment, question, and contextual factors. Theory will have an important role to play, for example whether the processes are conceptualised as micro-level, or meso-level. Whatever buffer you choose, you should justify and make sure that you can defend when someone might ask about it, as the further you reach obviously the more crimes you will include, and these might alter your results.
So, let’s say we are interested in all crimes that occur within 400 metres of each licensed premise, within the study area. We chose 400 metres here as this is often the recommended distance for accessible bus stop guidance, so basically as far as people should walk to get to a bus stop. So in this case, we want to take our points, which represent the licensed premises, and build buffers of 400 metres around them.
You can do with the
st_buffer() function. We pass two arguments to our function: the item which we want to buffer (the points in our cc_bars object) and the size of this buffer. Let’s quickly illustrate:
prem_buffer <- st_buffer(cc_bars, 1)
You might get a warning here, saying "st_buffer does not correctly buffer longitude/latitude datadist is assumed to be in decimal degrees (arc_degrees).". This message indicates that sf assumes a distance value (our size of the buffer, specified as ’1’ above) is given in degrees. This is because we have our data in a Geographic Coordinate System (lat/long data in WGS 84).
If we want to calculate the size of our buffer in a meaningful distance on our 2D surfaces, we can transform to a Projected Coordinate System, such as British National Grid. Let’s do this now:
#The code for BNG is 27700
bars_bng <- st_transform(cc_bars, 27700)
Now we can try again, with metres, specifying our indicated 400 metres radius:
bars_buffer <- st_buffer(bars_bng, 400)
Let’s see how that looks:
plot(st_geometry(bars_buffer))
plot(st_geometry(bars_bng), add = T)
That should look nice and squiggly. We can see it looks like there is quite a lot of overlap here. Should we maybe consider smaller buffers? Let’s look at 100 metre buffers:
bar_buffer_100 <- st_buffer(bars_bng, 100) # create 100m buffer
# plot new buffers
plot(st_geometry(bar_buffer_100))
plot(st_geometry(bars_bng), add = T)
Quite a bit of overlap, but this is possibly because the licensed premises are close together in the City Centre. We will discuss how to deal with this later on. For now, let’s go with these 100 metre buffers, and where a crime falls into an area of overlap, we will count it towards both premises. Where a crime falls within multiple buffers, it will count towards all the bars associated with those buffers.
The next step will be to count the number of crimes which fall into each buffer. Before we move on though, remember the CRS for our crimes is WGS 84 here, so we will need to convert our buffer layer back to this:
buffer_WGS84 <- st_transform(bar_buffer_100, 4326)
Now let’s just have a look:
plot(st_geometry(buffer_WGS84))
plot(st_geometry(cc_crimes), col = 'blue', add = T)
Some crimes fall inside some buffers, others not so much. Well, let’s get to our next spatial operation to be able to determine how many crimes happened in the 100 metre radius of each bar in Manchester City Centre.
Subsection 2.4.4 Counting points within a polygon
When you have a polygon layer and a point layer and want to know how many or which of the points fall within the bounds of each polygon, you can use this method of analysis. In computational geometry, the point-in-polygon (PIP) problem asks whether a given point in the plane lies inside, outside, or on the boundary of a polygon. As you can see, this is quite relevant to our problem, wanting to count how many crimes (points) fall within 100 metres of our licensed premises (our buffer polygons). We can achieve this with the
st_join() function, which spatially joins the bar name to each crime, and then we could find the frequency of each name (hence the count(name)). This returns a frequency table of the number of crimes within the buffer of each bar, saved in the crimes_per_prem object.
crimes_per_prem <- cc_crimes %>%
st_join(buffer_WGS84, ., left = FALSE) %>%
count(name)
You now have a new dataframe,
crimes_per_prem, which has a column for the name of the bars, a column for the number of violent crimes that fall within the buffer, and a column for the geometry. Take a moment to look at this table. Use the View() function. Which premises have the most violent crimes? Let’s see the bar with the most crimes:
crimes_per_prem %>%
filter(n == max(n)) %>%
select(name, n)
The bar with the highest number of crimes will be shown in the output. Keep this in mind for the next section.
So in this case, we used the PIP approach, and counted the number of points which fell into each polygon. We saw earlier, with the buffers, that they often overlapped with one another. This means that a crime may have been counted multiple times. This resulting data therefore tells us: How many crimes happened within 100 metres of each bar. This is one way to approach the problem, but not the only way. In our next spatial operation, we will calculate distances in order to explore another way.
Subsection 2.4.5 Distances: Finding the nearest point
Another way to solve this problem is to assign each crime (point) to the closest possible bar (other point). That is, look at the distances between crimes’ locations and the locations of all the bars in Manchester, and then, from those, choose the bar which is the closest. Then, we can assign this bar as the location for that crime.
We can achieve this using the
st_nearest_feature() function. This function takes our two sf objects, and for each row of the first one (x = cc_crimes), simply returns the index of the nearest features from the second one (y = cc_bars). We combine with the mutate() function in order to create a new variable which contains this index for each crime. Let’s illustrate:
crime_w_bars <- cc_crimes %>%
mutate(nearest_bar = st_nearest_feature(cc_crimes, cc_bars))
If we now have a look at this new object "crime_w_bars", we can see it is our crimes data. But we have a new column, which contains the index of the closest bar in the cc_bars dataframe, right at the end. So for example, the first point there, the nearest bar is that in location 84 (vectors in R are 1-indexed, not 0-indexed like many other languages). If we wanted to look at what is on the 84th row, we may call:
cc_bars[84]
Note: you may get a different number, if since time of writing new bars have been added, or taken away in the OSM database. You can continue with a different number to follow along!
However, this returns all the 90+ variables for this row. If we want only the name, we can query for the 84th row and the 2nd column (which is
name):
cc_bars[84, 2]
You can see the name is "Dive". For that first crime, in our dataset, the nearest bar is "Dive" bar. Now, instead of going through this process manually for each point, we can use the index to subset within our
mutate() function:
crime_w_bars <- cc_crimes %>%
mutate(nearest_bar = cc_bars[st_nearest_feature(cc_crimes, cc_bars),2])
Now we have new information in this
nearest_bar column, the name of the nearest bar, and the geometry. We actually don’t need the geometry for now, as we will simply be counting the frequency of each bar, which we can join back to our cc_bars object, which has a geometry, so we can extract the $name element only, and remove the geometry. Like so:
crimes_per_prem_2 <- crime_w_bars %>% # create new crimes_per_prem_2 object
st_drop_geometry() %>% # drop (remove) the geometry
group_by(nearest_bar$name) %>% # group by to find frequency of each bar
summarise(num_crimes = n()) %>% # count number of crimes
rename(name = `nearest_bar$name`) # rename variable to 'name'
To tie this back to our spatial object "cc_bars", we can use the
left_join() function:
crimes_per_prem_2 <- left_join(cc_bars, crimes_per_prem_2,
by = c("name" = "name"))
Let’s see the bar with the most crimes with this approach:
crimes_per_prem_2 %>%
filter(num_crimes == max(num_crimes, na.rm = TRUE)) %>%
select(name, num_crimes)
## Simple feature collection with 1 feature and 2 fields ## Geometry type: POINT ## Dimension: XY ## Bounding box: xmin: -2.236 ymin: 53.48 xmax: -2.236 ymax: 53.48 ## Geodetic CRS: WGS 84 ## name num_crimes geometry ## 1 Crafty Pig 50 POINT (-2.236 53.48)
The bar with the highest number of crimes is still the same bar, but the number of crimes differs between the two approaches. This means there is a difference in the number of crimes attributed to this bar. Some crimes fell within the buffer in the first approach, but were closer to another bar in the dataset, and were instead attributed to that one using this approach.
So which is better? This is once again up to you as the researcher and analyst to decide. They do slightly different things, and so will answer slightly different questions. With the nearest feature approach, instead of talking about the number of crimes within some distance to the bar, we are instead talking about, for each crime, the closest venue. This might mean that we could be attributing crimes that happen quite far from the venue to it, just because it’s the closest within our dataset. However, we are counting each crime only once. Pros and cons need to be weighed to make decisions like these.
Subsection 2.4.6 Measuring distances
Let’s have a look at this bar called "Crafty Pig". We can select this from the
cc_bars, buffers, and crimes
cp <- cc_bars %>% filter(name == "Crafty Pig")
cp_buffer <- bar_buffer_100 %>% filter(name == "Crafty Pig") %>%
st_transform(4326) # transform CRS
cp_crimes <- crime_w_bars %>% filter(nearest_bar$name == "Crafty Pig")
We can use
mapply() and the st_union() function to draw a line between each crime and the closest bar (Crafty Pig in this case):
dist_lines <- st_sfc( # Create simple feature geometry list column
mapply( # apply function to multiple objects
function(a,b){ # specify function parametres
st_cast(st_union(a,b),"LINESTRING") # specify function
},
cp_crimes$geometry, # input a for function
cp_crimes$nearest_bar$geometry, # input b for function
SIMPLIFY=FALSE)) # don't attempt to reduce the result
We can then plot these to get an idea of what we’re looking at:
plot(st_geometry(cp_buffer))
plot(st_geometry(cp), col = "black", add = TRUE)
plot(st_geometry(cp_crimes), col = "blue", add = TRUE)
plot(st_geometry(dist_lines), add = TRUE)
Note: make sure to run all 4 lines above in one batch.
So we can see that all these crimes happened only one location, which is within the 100 metre buffer. But how far exactly are they?
You can use the
st_distance() function to answer this question. We wrap this in the mutate() function in order to create a new column called distance which will contain for each row (each crime) the distance between that and its nearest bar (in this case Crafty Pig).
cp_crimes <- cp_crimes %>%
mutate(distance = st_distance(geometry, nearest_bar$geometry))
Having a look at our newly created variable, we can see the distance of the crime locations from the Crafty Pig bar.
One thing you might find strange about the data is: why are all these crimes geocoded on top of one another? This is how the open data are released, using geo-masking by snapping crime locations to a geo-mask (a set of points). This is done to ensure anonymity in the data (see [119] for more detail on this). In non-anonymised data, you might expect to see a little less overlap in your crime locations. Then with variation in distances between crimes and their nearest bars, we could use these distances to inform a buffer width for example. Anyway, we will return to distances a little later with a better dataset. But now, let’s move on to putting these outcomes on a map, which will help us further investigate the case of the Crafty Pig!
