13.5 Define metropolitan areas

We define metropolitan areas as pixels of 20 km² inhabited by more than 500,000 people.Pixels at this coarse resolution can rapidly be created using aggregate(), as introduced in Section 5.3.3.The command below uses the argument fact = 20 to reduce the resolution of the result twenty-fold (recall the original raster resolution was 1 km²):

pop_agg = aggregate(reclass$pop, fact = 20, fun = sum)

The next stage is to keep only cells with more than half a million people.

pop_agg = pop_agg[pop_agg > 500000, drop = FALSE] 

Plotting this reveals eight metropolitan regions (Figure 13.2).Each region consists of one or more raster cells.It would be nice if we could join all cells belonging to one region.raster’s clump() command does exactly that.Subsequently, rasterToPolygons() converts the raster object into spatial polygons, and st_as_sf() converts it into an sf-object.

polys = pop_agg %>% 
  clump() %>%
  rasterToPolygons() %>%
  st_as_sf()

polys now features a column named clumps which indicates to which metropolitan region each polygon belongs and which we will use to dissolve the polygons into coherent single regions (see also Section 5.2.6):

metros = polys %>%
  group_by(clumps) %>%
  summarize()

Given no other column as input, summarize() only dissolves the geometry.


Figure 13.2: The aggregated population raster (resolution: 20 km) with the identified metropolitan areas (golden polygons) and the corresponding names.

The resulting eight metropolitan areas suitable for bike shops (Figure 13.2; see also code/13-location-jm.R for creating the figure) are still missing a name.A reverse geocoding approach can settle this problem.Given a coordinate, reverse geocoding finds the corresponding address.Consequently, extracting the centroid coordinate of each metropolitan area can serve as an input for a reverse geocoding API.The revgeo package provides access to the open source Photon geocoder for OpenStreetMap, Google Maps and Bing.By default, it uses the Photon API.revgeo::revgeo() only accepts geographical coordinates (latitude/longitude); therefore, the first requirement is to bring the metropolitan polygons into an appropriate coordinate reference system (Chapter 6).

metros_wgs = st_transform(metros, 4326)
coords = st_centroid(metros_wgs) %>%
  st_coordinates() %>%
  round(4)

Choosing frame as revgeocode()’s output option will give back a data.frame with several columns referring to the location including the street name, house number and city.

library(revgeo)
metro_names = revgeo(longitude = coords[, 1], latitude = coords[, 2], 
                     output = "frame")

To make sure that the reader uses the exact same results, we have put them into spDataLarge:

# attach metro_names from spDataLarge
data("metro_names", package = "spDataLarge")

Overall, we are satisfied with the city column serving as metropolitan names (Table 13.2) apart from one exception, namely Wülfrath which belongs to the greater region of Düsseldorf.Hence, we replace Wülfrath with Düsseldorf (Figure 13.2).Umlauts like ü might lead to trouble further on, for example when determining the bounding box of a metropolitan area with opq() (see further below), which is why we avoid them.

metro_names = dplyr::pull(metro_names, city) %>% 
  as.character() %>% 
  ifelse(. == "Wülfrath", "Duesseldorf", .)