Measuring distance more thoroughly

Section 2.7 Measuring distance more thoroughly

Before we end the chapter, we want to return to the spatial operation of measuring the distance between points. This may be important, for example to those researchers focusing on studying the journey to crime by offenders. In this area, a common parameter studied is the average distance to crime from their home locations. In order to estimate these parametres, we first need to have a way to generate the distances. In this section, we will use another data set (this time from Madrid, Spain) to show a simpler example to look at the issue of geographical distance.

🔗

Subsection 2.7.1 How far are police stations in Madrid?

To illustrate how to measure distance, we will download data from the city of Madrid in Spain. Specifically we will obtain a csv file with the latitude and longitude of the police stations and a geoJSON file with the administrative boundary for the city of Madrid. Both are available from the data provided with this book (see Preamble section). We will also turn the .csv into a sf object with the appropriate coordinate reference system for this data, following the steps we’ve outlined in Chapter 1 and earlier in this chapter.

🔗

#read csv data
comisarias <- read_csv("data/nationalpolice.csv")
#set crs, read into sf object, and assign crs
polCRS <- st_crs(4326)
comisarias_sf <- st_as_sf(comisarias, coords = c("X", "Y"), crs = polCRS)

#create unique id for each row
comisarias_sf$id <- as.numeric(rownames(comisarias_sf))

#Read as sf boundary data for Madrid city
madrid <- st_read("data/madrid.geojson")

We went through a lot of steps there, so it’s worth to check in and plot our data to make sure that everything looks the way we expect. To practice with leaflet map some more, let’s plot using leaflet() function:

🔗

leaflet(comisarias_sf) %>% addTiles() %>%
  addMarkers(data = comisarias_sf) %>%
  addPolygons(data = madrid)

We can clearly see here that there are areas of the municipality that are far away from any national police station, the North West part of the city which you can see is a green area but noticeably also the South East, which is mostly urban and in fact is the location of a known shanty town and open drug market ("Cañada Real", you can read about it in the award-winning book by [120]).

🔗

Subsection 2.7.2 Distance in geographical space

There are many definitions of distance in data science and spatial data science. A common definition of distance is the Euclidean distance, which simply is the length of a segment connecting two points in a two-dimensional place. Because of the distortions caused by projections on a flat surface, a straight line on a map is not necessarily the shortest distance. Thus, another common definition used in geography is the great circle distance, which corresponds to an arc linking two points on a sphere and takes into account the spherical shape of the world. The great circle distance is useful, for example, to evaluate the shortest path when intercontinental distances are concerned. For applications which require a network distance measure, Euclidean distance is generally inadequate. Instead, Manhattan distance (also called Taxicab, Rectilinear, or City block metric) is a commonly used metric. We return to dealing with spatial data on networks in Chapter 8.

🔗

We can compute both with the st_distance function of the sf package. This function can be used to measure the distance between two points, between one point and others or between all points. In the latter case we obtain a symmetric matrix of distances (\(N \times N\)), taken pairwise between the points in our dataset. In the diagonal we find the combinations between the same points giving all null values.

🔗

Say, we want to measure the distance between the main police headquarters ("Jefatura Superior de Policia", row 34) and three other stations (say row 1, row 10, and row 25 in our dataset). We could use the following code for that:

🔗

# calculate distance
dist_headquarters <- st_distance(slice(comisarias_sf, 34),
                              slice(comisarias_sf, c(1, 10, 25)))

dist_headquarters # distance in metres

## Units: [m]
##      [,1] [,2] [,3]
## [1,] 8124 5086 8555

The result is a matrix with a single row or column (depending on the order of the spatial objects) with a class of units. Often we may want to re-express these distances in a different unit. For this purpose the units package offers useful functionality, through the set_units() function.

🔗

set_units(dist_headquarters, "km")

## Units: [km]
##       [,1]  [,2]  [,3]
## [1,] 8.124 5.086 8.555

We can compute the distance between all police stations as well.

🔗

# calculate distance
m_distance <- st_distance(comisarias_sf)

# matrix dimensions
dim(m_distance)

## [1] 34 34

If you want to preview the top of the matrix, you can use:

🔗

head(m_distance)

🔗

Subsection 2.7.3 A practical example to evaluate distance

For this practical example, we will look at the Madrid data. Earlier, we read in the relevant geoJSON file and stored in the object madrid. The CRS is EPSG 4326 and therefore it is geographic projection, with distance expressed in degrees.

🔗

st_crs(madrid)

For this we will transform and reproject to EPSG 2062, which is one of the appropriate projected coordinate systems for Madrid. How do you figure out this if you don’t want to check the EPSG site? There is a convenient package crsuggest developed by Kyle [121] that does this job for you. The function suggest_crs() from this package will attempt to list of candidate projections for your sf object. And we will see that top of that list is the 2062 EPSG.

🔗

suggested_crs <- suggest_crs(madrid)
head(suggested_crs, n = 1)

## # A tibble: 1 × 6
##   crs_code crs_name          crs_type crs_gcs crs_units
##   <chr>    <chr>             <chr>      <dbl> <chr>    
## 1 2062     Madrid 1870 (Mad… project…    4903 m        
## # … with 1 more variable: crs_proj4 <chr>

Now we can transform our madrid object appropriately:

🔗

madrid_metres <- st_transform(madrid, crs = 2062)

Before we saw that, some areas of Madrid are nowhere near a police station. Let’s say we want to get precise about this and we want to know how far different parts of the city of Madrid are from a police station, and we want to be able to show this in a map. Solving this means we have to define "parts of the city". What we will do is to divide the city of Madrid into different cells of 250 metres within a grid using the st_make_grid function of the sf package.

🔗

madrid_grid <- st_make_grid(madrid_metres,  cellsize = 250)

#only extract the points in the limits of Madrid
madrid_grid <- st_intersection(madrid_grid, madrid_metres)

We can plot the results, to see that everything has gone according to plan:

🔗

plot(madrid_grid)

A black, filled shape of Madrid. — Figure 2.7.1. Madrid divided into 250 metres grid cells
🔗

With so many "small" cells, the paper printed version of this map may just look completely black. But it is simply composed of many small cells. So how do we look at distance from police stations here? We can measure the distance between each grid cell to all 39 police stations. To estimate the distance to the nearest police station, we will find the minimum distance value for each grid, i.e., the distance to the nearest station.

🔗

comisarias_sf_metres <- st_transform(comisarias_sf, crs = 2062)

distances <- st_distance(comisarias_sf_metres,
                         st_centroid(madrid_grid)) %>%
  as_tibble()

If you view the new object "distances", you will see there is a row for each police station and a column representing each of the 10082 cells in our grid. For using these distances in a leaflet map, we will project back into 4326. And then we will compute the shortest distance for each cell.

🔗

# Compute distances
police_distances <- data.frame(
  # We want grids in a WGS 84 CRS:
  us = st_transform(madrid_grid, crs = 4326),
  # Extract minimum distance for each grid
  distance_km = map_dbl(distances, min)/1000,
  # Extract the value's index for joining with the location info
  location_id = map_dbl(distances, function(x) match(min(x), x))) %>%
  # Join with the police station table
  left_join(comisarias_sf, by = c("location_id" = "id"))

We now have a dataframe of distances, for each grid, to the nearest police station. We can have a look at the distribution of these distances through plotting them on a histogram

🔗

# Plot and examine distances
hist(police_distances$distance_km, main = 'Distance to nearest police station')

A bar chart, titled ’Distance to nearest police station’. The vertical axis is labeled Frequency, and ranges from zero to fifteen hundred, with notches at steps of five hundred. The horizontal axis represents distance in kilometres, with one bar for each kilometre range, and labels every five steps from zero to fifteen. The plot peaks at the second bar, one to two kilometres, at around seventeen hundred. The first and third bars are next, at around thirteen hundred. From there, the bars taper down as distance increases, with possible plateaus from five to nine kilometres, and again at nine to fourteen. The last bar, corresponding to eighteen to nineteen kilometres is very short. — Figure 2.7.2. Histogram of distance to nearest police station
🔗

Now we are ready to use this data to plot a map. We can get creative, and use the makeIcon() function from leaflet to assign an image as our icon, rather than these regular pointer icons we’ve been seeing. To do this, first we will adjust some aesthetics.

🔗

We get the URL for the icon (we break into part 1 (_pt1) and part 2 (_pt2) so it can be seen in the textbook), and pass this to the makeIcon() function.

🔗

# Create more appropriate icon, taking it from Wikipedia commons
icon_url_pt1 <- "https://upload.wikimedia.org/wikipedia/commons/"
icon_url_pt2 <- "a/ad/189-woman-police-officer-1.svg"

# create icon, adjusting size
police_icon <- makeIcon(paste0(icon_url_pt1, icon_url_pt2),
                        iconWidth = 12, iconHeight = 20)

Now we also want to create a colour scale. We can, for example, group our grids into quantiles. To do so, we can use the quantile() function:

🔗

# Bin ranges for a nicer colour scale
bins <- quantile(police_distances$distance_km)
# Create a binned colour palette
pal <- colorBin(c("#0868AC", "#43A2CA", "#7BCCC4", "#BAE4BC", "#F0F9E8"),
                domain = police_distances$distance_km,
                bins = bins, reverse = TRUE)

Now let’s create the map.

🔗

full_map <- leaflet() %>%
  addTiles() %>%
  addMarkers(data = comisarias_sf, icon = ~police_icon,
             group = "Police stations") %>%
  addPolygons(data = police_distances[[1]],
              fillColor = pal(police_distances$distance_km),
              fillOpacity = 0.8, weight =0, opacity = 1, color = "transparent",
              group = "Distances",
              highlight = highlightOptions(weight = 2.5, color = "#666",
                                           bringToFront = TRUE, opacity= 1),
              popupOptions = popupOptions(autoPan = FALSE, closeOnClick = TRUE,
                                          textOnly = T)) %>%
   addLegend(pal = pal, values = (police_distances$distance_km),
             opacity = 0.8, title = "Distance (Km)", position= "bottomright")

# print the map
full_map

And there you go. Just remember something. It is easy to misinterpret data and maps. You always need to care a great deal about measurement, quality of your data, and other potential issues affecting interpretation. When it comes to distance, and the movements of people and law enforcement personnel, for example, physical distance is not trivial, but time to arrival is also important and this is determined by factors other than Euclidean distance (e.g., availability and speed of transport, physical barriers, etc.). Our representation is always as good as the data we have. In Spain there are two other police forces (Guardia Civil, that patrols rural areas, and municipal civil, with jurisdiction for local administrative enforcement) that we are not representing here (that is, our data is incomplete). And we are not plotting the police stations in the nearby municipalities that are part of Madrid metropolitan area, around the edges.

🔗

Prev Top Next