Merge pull request #1175 from geocompx/refactor_location

Refactor Chapter 14 to use z22 Census 2022 data

Author: jannes-m

.gitignore

@@ -17,3 +17,6 @@ libs/
 geocompr2.rds
 *.gpkg
 figures/
+.claude/
+CLAUDE.md
+*_cache/

14-location.Rmd

@@ -11,118 +11,116 @@ library(osmdata)
 library(spDataLarge)
 ```
 
-E1. Download the csv file containing inhabitant information for a 100 m cell resolution (https://www.zensus2011.de/SharedDocs/Downloads/DE/Pressemitteilung/DemografischeGrunddaten/csv_Bevoelkerung_100m_Gitter.zip?__blob=publicationFile&v=3).
-Please note that the unzipped file has a size of 1.23 GB.
-To read it into R, you can use `readr::read_csv`.
-This takes 30 seconds on a machine with 16 GB RAM.
-`data.table::fread()` might be even faster, and returns an object of class `data.table()`.
-Use `dplyr::as_tibble()` to convert it into a tibble.
-Build an inhabitant raster, aggregate it to a cell resolution of 1 km, and compare the difference with the inhabitant raster (`inh`) we have created using class mean values.
+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.
 
 ```{r, 14-ex-e1, eval=FALSE}
 # Coarse inhabitant raster (1 km resolution)
 #*******************************************
 
-# inhabitant raster (coarse resolution); this is one of the results of the 
-# previous exercise
+# Load 1 km population data from spDataLarge
 data("census_de", package = "spDataLarge")
-input = select(census_de, x = x_mp_1km, y = y_mp_1km, pop = Einwohner,
-                      women = Frauen_A, mean_age = Alter_D, hh_size = HHGroesse_D)
-input_tidy = dplyr::mutate(input, dplyr::across(.fns = ~ifelse(. %in% c(-1, -9), NA, .)))
-input_ras = terra::rast(input_tidy, type = "xyz", crs = "EPSG:3035")
-inh_coarse = input_ras$pop
-# reclassify, i.e. convert the classes into inhabitant numbers using class means
-rcl = matrix(c(1, 1, 125, 2, 2, 375, 3, 3, 1250, 4, 4, 3000, 5, 5, 6000,
-               6, 6, 8000), ncol = 3, byrow = TRUE)
-inh_coarse = terra::classify(inh_coarse, rcl = rcl, right = NA)
+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)
 #******************************************
-url =
-  paste0("https://www.zensus2011.de/SharedDocs/Downloads/DE/Pressemitteilung/",
-         "DemografischeGrunddaten/csv_Bevoelkerung_100m_Gitter.zip", 
-         "?__blob=publicationFile&v=3")
-# download fine raster
-download.file(url = url, destfile = file.path(tempdir(), "census.zip"),
-              method = "auto", mode = "wb")
-# list the file names
-nms = unzip(file.path(tempdir(), "census.zip"), list = TRUE)
-# unzip only the csv file
-base_name = grep(".csv$", nms$Name, value = TRUE)
-unzip(file.path(tempdir(), "census.zip"), files = base_name, exdir = tempdir())
-# read in the csv file
-input = data.table::fread(file.path(tempdir(), base_name)) |>
-  dplyr::as_tibble()
-input = select(input, x = starts_with("x_mp_1"),
-                      y = starts_with("y_mp_1"), inh = Einwohner)
-# set -1 and -9 to NA
-input = dplyr::mutate(input,
-                      dplyr::across(.fns = ~ifelse(. %in% c(-1, -9), NA, .)))
-# convert table into a raster (x and y are cell midpoints)
-inh_fine = terra::rast(input, type = "xyz", crs = "EPSG:3035")
-# Note that inh_fine contains the actual number of inhabitants per raster cell
-# instead of mean class values as was the case with its coarse 1km counterpart
+
+# 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 = terra::aggregate(
-  inh_fine, fact = terra::res(inh_coarse)[1] / terra::res(inh_fine)[1], 
+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) = terra::origin(inh_coarse)
+terra::origin(inh_fine_agg) = terra::origin(inh_coarse)
 # make the comparison
-summary(inh_fine - inh_coarse)
-plot(inh_fine - inh_coarse)
-plot(abs(inh_fine - inh_coarse) > 1000)
-# the biggest deviations can be found in big cities like Berlin
-terra::global((abs(inh_fine - inh_coarse) > 1000), fun = "sum", na.rm = TRUE)
-# 18,121 cells have a deviation > 1000 inhabitants
-terra::global((abs(inh_fine - inh_coarse) > 5000), fun = "sum", na.rm = TRUE)
-# 338 cells have a deviation > 5000
+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. 
+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.
 
 ```{r, 14-ex-e2, eval=FALSE}
-# Here, we assue that you have already created `input_ras` in the first exercise.
+# 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.
-rcl_pop = matrix(c(1, 1, 127, 2, 2, 375, 3, 3, 1250, 
-                   4, 4, 3000, 5, 5, 6000, 6, 6, 8000), 
-                 ncol = 3, byrow = TRUE)
-rcl_women = matrix(c(1, 1, 3, 2, 2, 2, 3, 3, 1, 4, 5, 0), 
-                   ncol = 3, byrow = TRUE)
-# here we are giving the classes (3 to 5) containing the oldest people the
-# highest weight
-rcl_age = matrix(c(1, 1, 1, 2, 2, 1, 3, 5, 3),
-                 ncol = 3, byrow = TRUE)
-rcl_hh = rcl_women
-rcl = list(rcl_pop, rcl_women, rcl_age, rcl_hh)
-
-reclass = input_ras
-for (i in 1:terra::nlyr(reclass)) {
+
+# 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) = names(input_ras)
+names(reclass) = demo_vars
 
-# The rest of the analysis follows exactly the code presented in the book. 
+# The rest of the analysis follows exactly the code presented in the book.
 
 # Add metro names to metros sf object
 #************************************
-metro_names = dplyr::pull(metro_names, city) |>
+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(reclass$pop, fact = 20, fun = sum, na.rm = TRUE)
-pop_agg = pop_agg[pop_agg > 500000, drop = FALSE] 
+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) |>
@@ -132,7 +130,7 @@ polys = pop_agg |>
 metros = polys |>
   dplyr::group_by(patches) |>
   dplyr::summarize()
-metros$metro_names = metro_names
+metros$metro_names = metro_names_vec
 
 # Create shop/poi density raster
 #*******************************
@@ -142,23 +140,22 @@ 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), 
+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)  
+rcl_poi = cbind(rcl_poi, 0:3)
 # reclassify
-poi = terra::classify(poi, rcl = rcl_poi, right = NA) 
+poi = terra::classify(poi, rcl = rcl_poi, right = NA)
 names(poi) = "poi"
-# remove population raster and add poi raster
-reclass = reclass[[names(reclass) != "pop"]] |>
-  c(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[metro_names == "Berlin", ]
+berlin = metros[metros$metro_names == "Berlin", ]
 berlin_raster = terra::crop(result, berlin)
 # summary(berlin_raster)
 # berlin_raster
@@ -168,6 +165,6 @@ berlin_raster[berlin_raster == 0] = NA
 leaflet::leaflet() |>
   leaflet::addTiles() |>
   leaflet::addRasterImage(raster::raster(berlin_raster), colors = "darkgreen", opacity = 0.8) |>
-  leaflet::addLegend("bottomright", colors = c("darkgreen"), 
+  leaflet::addLegend("bottomright", colors = c("darkgreen"),
                      labels = c("potential locations"), title = "Legend")
 ```