13.6 Points of interest

The osmdata package provides easy-to-use access to OSM data (see also Section 7.2).Instead of downloading shops for the whole of Germany, we restrict the query to the defined metropolitan areas, reducing computational load and providing shop locations only in areas of interest.The subsequent code chunk does this using a number of functions including:

• map() (the tidyverse equivalent of lapply()), which iterates through all eight metropolitan names which subsequently define the bounding box in the OSM query function opq() (see Section 7.2).
• add_osm_feature() to specify OSM elements with a key value of shop (see wiki.openstreetmap.org for a list of common key:value pairs).
• osmdata_sf(), which converts the OSM data into spatial objects (of class sf).
• while(), which tries repeatedly (three times in this case) to download the data if it fails the first time.74Before running this code: please consider it will download almost 2GB of data.To save time and resources, we have put the output named shops into spDataLarge.To make it available in your environment either attach the spDataLarge package or run data("shops", package = "spDataLarge").

shops = map(metro_names, function(x) {
  message("Downloading shops of: ", x, "\n")
  # give the server a bit time
  Sys.sleep(sample(seq(5, 10, 0.1), 1))
  query = opq(x) %>%
    add_osm_feature(key = "shop")
  points = osmdata_sf(query)
  # request the same data again if nothing has been downloaded
  iter = 2
  while (nrow(points$osm_points) == 0 & iter > 0) {
    points = osmdata_sf(query)
    iter = iter - 1
  }
  points = st_set_crs(points$osm_points, 4326)
})

It is highly unlikely that there are no shops in any of our defined metropolitan areas.The following if condition simply checks if there is at least one shop for each region.If not, we recommend to try to download the shops again for this/these specific region/s.

# checking if we have downloaded shops for each metropolitan area
ind = map(shops, nrow) == 0
if (any(ind)) {
  message("There are/is still (a) metropolitan area/s without any features:\n",
          paste(metro_names[ind], collapse = ", "), "\nPlease fix it!")
}

To make sure that each list element (an sf data frame) comes with the same columns, we only keep the osm_id and the shop columns with the help of another map loop.This is not a given since OSM contributors are not equally meticulous when collecting data.Finally, we rbind all shops into one large sf object.

# select only specific columns
shops = map(shops, dplyr::select, osm_id, shop)
# putting all list elements into a single data frame
shops = do.call(rbind, shops)

It would have been easier to simply use map_dfr().Unfortunately, so far it does not work in harmony with sf objects.Please note that the shops object is also available in the spDataLarge package:

data("shops", package = "spDataLarge")

The only thing left to do is to convert the spatial point object into a raster (see Section 5.4.3).The sf object, shops, is converted into a raster having the same parameters (dimensions, resolution, CRS) as the reclass object.Importantly, the count() function is used here to calculate the number of shops in each cell.

If the shop column were used instead of the osm_id column, we would have retrieved fewer shops per grid cell.This is because the shop column contains NA values, which the count() function omits when rasterizing vector objects.

The result of the subsequent code chunk is therefore an estimate of shop density (shops/km²).st_transform() is used before rasterize() to ensure the CRS of both inputs match.

shops = st_transform(shops, proj4string(reclass))
# create poi raster
poi = rasterize(x = shops, y = reclass, field = "osm_id", fun = "count")

As with the other raster layers (population, women, mean age, household size) the poi raster is reclassified into four classes (see Section 13.4).Defining class intervals is an arbitrary undertaking to a certain degree.One can use equal breaks, quantile breaks, fixed values or others.Here, we choose the Fisher-Jenks natural breaks approach which minimizes within-class variance, the result of which provides an input for the reclassification matrix.

# construct reclassification matrix
int = classInt::classIntervals(values(poi), n = 4, style = "fisher")
int = round(int$brks)
rcl_poi = matrix(c(int[1], rep(int[-c(1, length(int))], each = 2), 
                   int[length(int)] + 1), ncol = 2, byrow = TRUE)
rcl_poi = cbind(rcl_poi, 0:3)  
# reclassify
poi = reclassify(poi, rcl = rcl_poi, right = NA) 
names(poi) = "poi"