14 Geomarketing

The solutions assume the following packages are attached (other packages will be attached when needed):

library(sf)
library(dplyr)
library(purrr)
library(terra)
library(osmdata)
library(spDataLarge)

E1. This exercise requires the z22 package for accessing 100 m resolution data. Install it with remotes::install_github("JsLth/z22"). Load the population data at 100 m cell resolution using z22::z22_data("population", res = "100m", year = 2022). Aggregate it to a cell resolution of 1 km using terra::aggregate() with fun = sum, and compare the result with the 1 km resolution data from census_de. Note that the 100 m data is much larger and may take some time to download.

# Coarse inhabitant raster (1 km resolution)
#*******************************************

# Load 1 km population data from spDataLarge
data("census_de", package = "spDataLarge")
pop_1km = census_de |>
  select(x, y, pop) |>
  mutate(pop = ifelse(pop < 0, NA, pop))
inh_coarse = terra::rast(pop_1km, type = "xyz", crs = "EPSG:3035")

# Fine inhabitant raster (100 m resolution)
#******************************************

# Load 100 m population data using z22 (this may take some time)
pop_100m = z22::z22_data("population", res = "100m", year = 2022, as = "df")
pop_100m = pop_100m |>
  rename(pop = cat_0) |>
  mutate(pop = ifelse(pop < 0, NA, pop))
inh_fine = terra::rast(pop_100m, type = "xyz", crs = "EPSG:3035")

# Comparing the coarse with the fine raster
#******************************************

# aggregate to the resolution of the coarse raster
inh_fine_agg = terra::aggregate(
  inh_fine, fact = terra::res(inh_coarse)[1] / terra::res(inh_fine)[1],
  fun = sum, na.rm = TRUE)
# origin has to be the same
terra::origin(inh_fine_agg) = terra::origin(inh_coarse)
# make the comparison
summary(inh_fine_agg - inh_coarse)
plot(inh_fine_agg - inh_coarse)
# Note: Since Census 2022 provides actual counts at both resolutions,
# differences should be minimal (mainly due to rounding or edge effects)
terra::global((abs(inh_fine_agg - inh_coarse) > 100), fun = "sum", na.rm = TRUE)

E2. Suppose our bike shop predominantly sold electric bikes to older people. Change the age raster accordingly, repeat the remaining analyses and compare the changes with our original result.

# Load data from spDataLarge
data("census_de", package = "spDataLarge")

input_tidy = census_de |>
  mutate(across(c(pop, women, mean_age, hh_size), ~ifelse(.x < 0, NA, .x)))
input_ras = terra::rast(input_tidy, type = "xyz", crs = "EPSG:3035")

# attach further necessary data
data("metro_names", "shops", package = "spDataLarge")

# Basically, we are assuming that especially older people will use an electric
# bike, therefore, we increase the weights for raster cells where predominantly
# older people are living.

# Reclassification matrices for continuous values
rcl_women = matrix(c(
  0, 40, 3,
  40, 47, 2,
  47, 53, 1,
  53, 60, 0,
  60, 100, 0
), ncol = 3, byrow = TRUE)

# For elderly electric bikes: give highest weights to older age groups
rcl_age = matrix(c(
  0, 40, 0,     # Young -> low weight
  40, 42, 0,    # Young-ish -> low weight
  42, 44, 1,    # Middle-aged -> some weight
  44, 47, 2,    # Older -> higher weight
  47, 120, 3    # Elderly -> highest weight
), ncol = 3, byrow = TRUE)

rcl_hh = matrix(c(
  0, 1.5, 3,
  1.5, 2.0, 2,
  2.0, 2.5, 1,
  2.5, 3.0, 0,
  3.0, 100, 0
), ncol = 3, byrow = TRUE)

rcl = list(rcl_women, rcl_age, rcl_hh)

# Separate population (used as counts for metro detection) from variables to reclassify
pop_ras = input_ras$pop
demo_vars = c("women", "mean_age", "hh_size")
reclass = input_ras[[demo_vars]]
for (i in seq_len(terra::nlyr(reclass))) {
  reclass[[i]] = terra::classify(x = reclass[[i]], rcl = rcl[[i]], right = NA)
}
names(reclass) = demo_vars

# The rest of the analysis follows exactly the code presented in the book.

# Add metro names to metros sf object
#************************************
metro_names_vec = dplyr::pull(metro_names, city) |>
  as.character() |>
  {\(x) ifelse(x == "Velbert", "Düsseldorf", x)}() |>
  {\(x) gsub("ü", "ue", x)}()

pop_agg = terra::aggregate(pop_ras, fact = 20, fun = sum, na.rm = TRUE)
pop_agg = pop_agg[pop_agg > 500000, drop = FALSE]

polys = pop_agg |>
  terra::patches(directions = 8) |>
  terra::as.polygons() |>
  sf::st_as_sf()

metros = polys |>
  dplyr::group_by(patches) |>
  dplyr::summarize()
metros$metro_names = metro_names_vec

# Create shop/poi density raster
#*******************************
shops = sf::st_transform(shops, sf::st_crs(reclass))
# create poi raster
poi = terra::rasterize(x = shops, y = reclass, field = "osm_id", fun = "length")
# 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 = terra::classify(poi, rcl = rcl_poi, right = NA)
names(poi) = "poi"
# add poi raster to demographic weights
reclass = c(reclass, poi)

# Identify suitable locations
#****************************
# calculate the total score
result = sum(reclass)

# have a look at suitable bike shop locations in Berlin
berlin = metros[metros$metro_names == "Berlin", ]
berlin_raster = terra::crop(result, berlin)
# summary(berlin_raster)
# berlin_raster
berlin_raster = berlin_raster > 9
berlin_raster[berlin_raster == 0] = NA
# make the plot
leaflet::leaflet() |>
  leaflet::addTiles() |>
  leaflet::addRasterImage(raster::raster(berlin_raster), colors = "darkgreen", opacity = 0.8) |>
  leaflet::addLegend("bottomright", colors = c("darkgreen"),
                     labels = c("potential locations"), title = "Legend")