14.5 Points of interest

poi (Points of Interest)

The osmdata package provides easy-to-use access to Open Street Map data

Caution: even filtering for metro areas in Germany, this query operates on 2Gb of data.

shops <- purrr::map(metro_names, function(x) {
  message("Downloading shops of: ", x, "\n")

  Sys.sleep(sample(seq(5, 10, 0.1), 1))
  query = osmdata::opq(x) |>
    osmdata::add_osm_feature(key = "shop")
  points = osmdata::osmdata_sf(query)

  iter = 2
  while (nrow(points$osm_points) == 0 && iter > 0) {
    points = osmdata_sf(query)
    iter = iter - 1
  }

  points$osm_points
})

# checking if we have downloaded shops for each metropolitan area
ind = purrr::map_dbl(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!")
}

# select only specific columns
shops = purrr::map_dfr(shops, select, osm_id, shop)

We will just use the convenient sample dataset with features in each of the metro areas.

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

shops |> 
  ggplot() +
  geom_sf(alpha = 0.01, shape = 20) +
  theme_minimal()

This spatial point object must be converted into a raster

shops <- sf::st_transform(shops, st_crs(reclass))

poi <- terra::rasterize(x = shops, y = reclass, field = "osm_id", fun = "length")

As with the other raster layers (population, women, mean age, household size) the poi raster is reclassified into four classes. Here again, some judgement is needed.

The authors choose the Fisher-Jenks natural breaks approach which minimizes within-class variance, the result of which provides an input for the reclassification matrix.

int <- classInt::classIntervals(terra::values(poi), n = 4, style = "fisher")

## Warning in classInt::classIntervals(terra::values(poi), n = 4, style =
## "fisher"): var has missing values, omitted in finding classes

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)  

poi <- terra::classify(poi, rcl = rcl_poi, right = NA) 

names(poi) = "poi"

poi

## class       : SpatRaster 
## dimensions  : 868, 642, 1  (nrow, ncol, nlyr)
## resolution  : 1000, 1000  (x, y)
## extent      : 4031000, 4673000, 2684000, 3552000  (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
## source(s)   : memory
## name        : poi 
## min value   :   0 
## max value   :   3