1. Introduction

Point-of-interest (POI) data from OpenStreetMap allow us to study how commercial and service activity is spatially organised within Warsaw. This paper examines four everyday POI categories — restaurants, cafés, pharmacies and supermarkets — and asks:

  1. Are these POI types randomly scattered across the city, or do they form spatially clustered patterns?
  2. Where are the local hotspots of each category, and how many distinct spatial clusters can be identified (k-means, DBSCAN)?
  3. How strongly is each category spatially concentrated relative to a perfectly uniform distribution (relative entropy / ETA)?
  4. Are POI categories located in similar parts of the city, or do they follow different spatial logics?
  5. How does POI density relate to population density? In particular, where are areas of high population but low POI density of a given type (potentially underserved areas), and where is the opposite true?

2. Packages, paths and constants

library(osmdata)      # OSM data download
library(sf)            # vector spatial data
library(dplyr)          # data manipulation
library(tidyr)          # reshaping
library(readr)          # csv I/O
library(ggplot2)         # visualisation
library(scales)          # readable scales / labels
library(spatstat)        # ppp objects, KDE, Clark-Evans test
library(dbscan)          # DBSCAN, kNN distances
library(viridis)         # colour palettes
library(gridExtra)       # arranging plots
library(ggdendro)        # dendrograms in ggplot format
library(purrr)           # functional helpers (map_dfr)
library(factoextra)
library(spatialWarsaw)  # QDC

Sys.setenv(LANG = "en")
options(scipen = 999)
set.seed(123)

# EPSG:2180 -- projected CRS for Poland, units in metres. Used for all
# distance-based calculations (KDE bandwidth, clustering, grid cell size).
# EPSG:4326 -- WGS84 lon/lat, used only for querying OSM.
crs_pl  <- 2180
crs_geo <- 4326

dir.create("data", showWarnings = FALSE)
dir.create("outputs", showWarnings = FALSE)

3. Study area: Warsaw

warsaw_boundary_path <- "data/warsaw_boundary.gpkg"

if (file.exists(warsaw_boundary_path)) {

  warsaw_boundary <- st_read(warsaw_boundary_path, quiet = TRUE)

} else {

  warsaw_boundary_raw <- getbb("Warszawa, Polska", format_out = "sf_polygon")

  # getbb() with format_out = "sf_polygon" can return either a single sf
  # object or a list with $multipolygon / $polygon elements, depending on
  # how many matching OSM relations exist. It handles both cases.
  if (inherits(warsaw_boundary_raw, "sf")) {
    warsaw_boundary <- warsaw_boundary_raw
  } else {
    warsaw_boundary <- warsaw_boundary_raw$multipolygon
    if (is.null(warsaw_boundary)) warsaw_boundary <- warsaw_boundary_raw[[1]]
  }

  warsaw_boundary <- warsaw_boundary %>%
    st_make_valid() %>%
    st_transform(crs_pl) %>%
    st_union() %>%
    st_sf(city = "Warszawa", geometry = .)

  st_write(warsaw_boundary, warsaw_boundary_path, delete_dsn = TRUE, quiet = TRUE)
}

ggplot() +
  geom_sf(data = warsaw_boundary, fill = "grey95", color = "grey20") +
  theme_minimal() +
  labs(title = "Study area: Warsaw administrative boundary")

4. POI data acquisition (OpenStreetMap)

poi_cache_path <- "data/warsaw_pois.gpkg"

get_poi_points <- function(key, value, category_name, bbox) {

  message("Downloading: ", category_name)

  q <- opq(bbox = bbox) %>%
    add_osm_feature(key = key, value = value)

  x <- osmdata_sf(q)

  objects <- list()

  if (!is.null(x$osm_points) && nrow(x$osm_points) > 0) {
    objects[[length(objects) + 1]] <- x$osm_points
  }

  if (!is.null(x$osm_lines) && nrow(x$osm_lines) > 0) {
    objects[[length(objects) + 1]] <- x$osm_lines %>% st_point_on_surface()
  }

  if (!is.null(x$osm_polygons) && nrow(x$osm_polygons) > 0) {
    objects[[length(objects) + 1]] <- x$osm_polygons %>% st_point_on_surface()
  }

  if (length(objects) == 0) {
    message("No objects found for: ", category_name)
    return(NULL)
  }

  # Use bind_rows to combine points, lines, and polygons regardless of extra columns
  result <- dplyr::bind_rows(objects) %>%
    distinct(osm_id, .keep_all = TRUE) %>%
    mutate(category = category_name) %>%
    # Standardize columns before returning
    select(osm_id, any_of("name"), category, geometry)

  message("Found ", nrow(result), " objects.")
  result
}

if (file.exists(poi_cache_path)) {

  pois_raw <- st_read(poi_cache_path, quiet = TRUE)

} else {

  bbox_warsaw <- st_bbox(st_transform(warsaw_boundary, crs_geo))

  poi_categories <- list(
    list(key = "amenity", value = "restaurant",  name = "restaurant"),
    list(key = "amenity", value = "cafe",        name = "cafe"),
    list(key = "amenity", value = "pharmacy",    name = "pharmacy"),
    list(key = "shop",    value = "supermarket", name = "supermarket")
  )

  poi_list <- lapply(poi_categories, function(cat) {
    get_poi_points(cat$key, cat$value, cat$name, bbox_warsaw)
  })

  # Remove any NULL entries if a category returned no features
  poi_list <- Filter(Negate(is.null), poi_list)

  # Combine all categories safely
  pois_raw <- dplyr::bind_rows(poi_list)

  st_write(pois_raw, poi_cache_path, delete_dsn = TRUE, quiet = TRUE)
}
# Keep only points that actually fall within the Warsaw boundary, and drop
# any geometry problems introduced by the point-on-surface conversion.
pois_sf <- pois_raw %>%
  st_transform(crs_pl) %>%
  st_filter(warsaw_boundary, .predicate = st_within) %>%
  mutate(category = factor(category,
                            levels = c("restaurant", "cafe",
                                       "pharmacy", "supermarket"))) %>%
  filter(!st_is_empty(geom))

poi_levels <- c("restaurant", "cafe", "pharmacy", "supermarket")
stopifnot(all(poi_levels %in% levels(pois_sf$category)))

poi_counts <- pois_sf %>%
  st_drop_geometry() %>%
  count(category, name = "n_points") %>%
  arrange(desc(n_points))

poi_counts
##      category n_points
## 1  restaurant     2791
## 2 supermarket     2391
## 3        cafe     1416
## 4    pharmacy      628

5. Population grid

We compare POI density with a 500 m GUS population grid

pop_grid_res  <- "500m"
pop_grid_path <- "C:/Users/Adam Pochmara/Downloads/NSP2021_TOT_grid500m_SHP"

population_grid_raw <- st_read(pop_grid_path, quiet = TRUE)

names(population_grid_raw)
## [1] "fid"      "PL_code"  "code"     "TOT"      "geometry"
pop_col <- "TOT"

stopifnot(pop_col %in% names(population_grid_raw))

population_grid <- population_grid_raw %>%
  st_make_valid() %>%
  st_transform(crs_pl) %>%
  st_filter(warsaw_boundary, .predicate = st_intersects) %>%
  rename(population = all_of(pop_col)) %>%
  mutate(population = replace_na(population, 0),
         cell_area_km2 = as.numeric(st_area(geometry)) / 1e6,
         pop_density_km2 = population / cell_area_km2,
         grid_id = row_number())

ggplot() +
  geom_sf(data = population_grid, aes(fill = pop_density_km2), color = NA) +
  geom_sf(data = warsaw_boundary, fill = NA, color = "grey15", linewidth = 0.4) +
  scale_fill_viridis_c(option = "C", trans = "sqrt", name = "Residents\nper km2") +
  theme_minimal() +
  labs(title = paste0("Population density in Warsaw (", pop_grid_res, " grid)"))

6. Descriptive overview of POI distributions

Restaurants and cafés are relatively more concentrated in the city centre, pharmacies appear to be distributed similarly to the population, while supermarkets show a more dispersed spatial pattern..

poi_map_colors <- c(
  restaurant = "#C44536",
  cafe = "#3B82F6",
  pharmacy = "#2A9D8F",
  supermarket = "#8A5CF6"
)

for (cat in poi_levels) {
  p <- ggplot() +
    geom_sf(data = warsaw_boundary, fill = "grey97", color = "grey40", linewidth = 0.3) +
    geom_sf(data = pois_sf %>% filter(category == cat),
            color = poi_map_colors[[cat]], alpha = 0.5, size = 0.45) +
    theme_minimal() +
    labs(title = paste0("POI locations in Warsaw: ", cat),
         subtitle = paste0("n = ", poi_counts$n_points[match(cat, poi_counts$category)]),
         x = NULL, y = NULL)
  print(p)
}

7. Are POI patterns random? Clark-Evans test

The Clark-Evans test compares the observed mean nearest-neighbour distance with the value expected under Complete Spatial Randomness (CSR). The R statistic is close to 1 for random patterns, below 1 for clustered patterns, and above 1 for regular (dispersed) patterns.

warsaw_window <- as.owin(st_union(warsaw_boundary))

make_ppp <- function(points_sf, window) {
  coords <- st_coordinates(points_sf)
  ppp(x = coords[, 1], y = coords[, 2],
      window = window, check = TRUE, checkdup = FALSE)
}

clark_evans_results <- pois_sf %>%
  st_drop_geometry() %>%
  distinct(category) %>%
  pull(category) %>%
  map_dfr(function(cat) {
    pts <- pois_sf %>% filter(category == cat)
    ppp_obj <- make_ppp(pts, warsaw_window)
    ce <- clarkevans.test(ppp_obj, correction = "none")
    tibble(category = cat,
           n_points = pts %>% nrow(),
           R_statistic = unname(ce$statistic),
           p_value = ce$p.value)
  })

clark_evans_results
## # A tibble: 4 × 4
##   category    n_points R_statistic p_value
##   <fct>          <int>       <dbl>   <dbl>
## 1 restaurant      2791       0.384       0
## 2 cafe            1416       0.363       0
## 3 pharmacy         628       0.637       0
## 4 supermarket     2391       0.229       0

For all four categories, the results provide strong evidence against Complete Spatial Randomness (CSR), indicating that the observed POI patterns are significantly more clustered than would be expected under a random spatial distribution.

8. Spatial clustering of POIs

8.1 K-means clustering of coordinates

K-means partitions each POI category into a fixed number of spatially coherent groups and returns cluster centres. It serves as a benchmark for what we can expect in the subsequent analyses.

poi_categories <- c(
  "restaurant",
  "cafe",
  "pharmacy",
  "supermarket"
)

elbow_plots <- lapply(poi_categories, function(cat) {

  coords <- pois_sf %>%
    filter(category == cat) %>%
    st_coordinates()

  fviz_nbclust(
    coords,
    kmeans,
    method = "wss",
    k.max = 12,
    nstart = 25
  ) +
    labs(
      title = cat,
      x = "Number of clusters (k)",
      y = "Within-cluster sum of squares"
    ) +
    theme_minimal()
})

grid.arrange(
  grobs = elbow_plots,
  ncol = 2,
  top = "Optimal number of clusters – k-means"
)

kmeans_clusters <- c(
  restaurant = 5,
  cafe = 5,
  pharmacy = 11,
  supermarket = 4
)

kmeans_by_category <- setNames(map(poi_levels, function(cat) {

  pts_sf <- pois_sf %>%
    filter(category == cat)

  coords <- st_coordinates(pts_sf)

  n_clusters <- kmeans_clusters[[cat]]

  km <- kmeans(
    coords,
    centers = n_clusters,
    nstart = 20
  )

  pts_sf$kmeans_cluster <- factor(km$cluster)

  centres_sf <- as.data.frame(km$centers) %>%
    st_as_sf(
      coords = c("X", "Y"),
      crs = crs_pl
    ) %>%
    mutate(category = cat)

  list(
    points = pts_sf,
    centres = centres_sf
  )

}), poi_levels)


for (cat in poi_levels) {

  n_clusters <- kmeans_clusters[[cat]]

  p <- ggplot() +
    geom_sf(
      data = warsaw_boundary,
      fill = "grey97",
      color = "grey40",
      linewidth = 0.3
    ) +
    geom_sf(
      data = kmeans_by_category[[cat]]$points,
      aes(color = kmeans_cluster),
      alpha = 0.45,
      size = 0.4
    ) +
    geom_sf(
      data = kmeans_by_category[[cat]]$centres,
      shape = 21,
      fill = "yellow",
      color = "black",
      size = 2.4
    ) +
    theme_minimal() +
    labs(
      title = paste0(
        "K-means: ", cat,
        " (k = ", n_clusters, ")"
      ),
      subtitle = "Yellow points mark cluster centres",
      x = NULL,
      y = NULL
    ) +
    guides(color = "none")

  print(p)
}

For restaurants and cafés, one central cluster and several peripheral clusters are visible. In the case of pharmacies, the pattern more closely corresponds to the city’s districts, while supermarkets do not form any clear overall pattern.

8.2 DBSCAN clustering (density-based)

DBSCAN finds dense groups without specifying the number of clusters and flags low-density observations as noise. Initially, the minimum number of points in a cluster was set at around 1% of the points for each type.

kNN_k <- c(
  restaurant = 28,
  cafe = 14,
  pharmacy = 6,
  supermarket = 24
)

knn_categories <- c(
  restaurant = 1600,
  cafe = 1400,
  pharmacy = 1300,
  supermarket = 1500
)

for (cat in poi_levels) {

  coords <- pois_sf %>%
    filter(category == cat) %>%
    st_coordinates()

  kNNdistplot(
    coords,
    k = kNN_k[[cat]]
  )

  abline(
    h = knn_categories[[cat]],
    col = "red",
    lty = 2
  )

  title(
    paste0(
      "kNN distance plot (k = ", kNN_k[[cat]], ") -- ", cat
    )
  )
}

dbscan_eps <- c(
  restaurant = 1600,
  cafe = 1400,
  pharmacy = 1300,
  supermarket = 1500
)

dbscan_minPts <- c(
  restaurant = 28,
  cafe = 14,
  pharmacy = 6,
  supermarket = 24
)

dbscan_by_category <- setNames(map(poi_levels, function(cat) {

  pts_sf <- pois_sf %>%
    filter(category == cat)

  coords <- st_coordinates(pts_sf)

  db <- dbscan(
    coords,
    eps = dbscan_eps[[cat]],
    minPts = dbscan_minPts[[cat]]
  )

  pts_sf$dbscan_cluster <- factor(db$cluster)

  pts_sf

}), poi_levels)

dbscan_points_all <- map_dfr(
  dbscan_by_category,
  identity
)

dbscan_summary <- dbscan_points_all %>%
  st_drop_geometry() %>%
  group_by(category) %>%
  summarise(
    n_points = n(),
    eps_m = dbscan_eps[[first(category)]],
    minPts = dbscan_minPts[[first(category)]],
    n_clusters = n_distinct(
      dbscan_cluster[dbscan_cluster != "0"]
    ),
    n_noise = sum(dbscan_cluster == "0"),
    share_noise = round(n_noise / n_points, 3),
    .groups = "drop"
  )

dbscan_summary
## # A tibble: 4 × 7
##   category    n_points eps_m minPts n_clusters n_noise share_noise
##   <fct>          <int> <dbl>  <dbl>      <int>   <int>       <dbl>
## 1 restaurant      2791  1600     28          4     112       0.04 
## 2 cafe            1416  1400     14          3      64       0.045
## 3 pharmacy         628  1300      6          5      41       0.065
## 4 supermarket     2391  1500     24         14      69       0.029
for (cat in poi_levels) {

  points_cat <- dbscan_points_all %>%
    filter(category == cat)

  p <- ggplot() +

    geom_sf(
      data = warsaw_boundary,
      fill = "grey97",
      color = "grey40",
      linewidth = 0.3
    ) +

    # noise
    geom_sf(
      data = points_cat %>%
        filter(dbscan_cluster == "0"),
      color = "grey70",
      alpha = 0.25,
      size = 0.35
    ) +

    # clustered points
    geom_sf(
      data = points_cat %>%
        filter(dbscan_cluster != "0"),
      aes(color = dbscan_cluster),
      alpha = 0.65,
      size = 1.0
    ) +

    theme_minimal() +

    labs(
      title = paste0("DBSCAN: ", cat),
      subtitle = paste0(
        "eps = ", dbscan_eps[[cat]],
        " m; minPts = ", dbscan_minPts[[cat]],
        " | yellow = cluster centres; grey = noise"
      ),
      x = NULL,
      y = NULL
    ) +

    guides(color = "none")

  print(p)
}

Because the city centre still dominated the clustering results, a simple tuning criterion was introduced. The combination of eps and minPts was selected by minimising the sum of two proportions: (1) the share of observations belonging to the largest cluster and (2) the share of observations classified as noise.

# ============================================================
# DBSCAN tuning
# ============================================================

dbscan_minPts_grid <- list(
  restaurant  = c(5, 10, 15, 20, 25, 30, 60, 90, 120, 150),
  cafe        = c(5, 10, 15, 20, 25, 30, 40, 60, 80, 100),
  pharmacy    = c(5, 10, 15, 20, 25, 30, 40, 50, 75, 100),
  supermarket = c(5, 10, 15, 20, 25, 30, 50, 75, 100, 125)
)

dbscan_eps_grid <- list(
  restaurant  = c(100, 200, 300, 400,500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000),
  cafe        = c(100, 200, 300, 400,500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000),
  pharmacy    = c(100, 200, 300, 400,500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000),
  supermarket = c(100, 200, 300, 400,500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000)
)


# ------------------------------------------------------------
# Test all parameter combinations for every category
# ------------------------------------------------------------

dbscan_tuning_results <- map_dfr(poi_levels, function(cat) {

  pts_sf <- pois_sf %>%
    filter(category == cat)

  coords <- st_coordinates(pts_sf)

  parameter_grid <- expand_grid(
    eps = dbscan_eps_grid[[cat]],
    minPts = dbscan_minPts_grid[[cat]]
  )

  map_dfr(seq_len(nrow(parameter_grid)), function(i) {

    eps_value <- parameter_grid$eps[i]
    minPts_value <- parameter_grid$minPts[i]

    db <- dbscan(
      coords,
      eps = eps_value,
      minPts = minPts_value
    )

    cluster_ids <- db$cluster

    n_points <- length(cluster_ids)

    n_noise <- sum(cluster_ids == 0)
    noise_share <- n_noise / n_points

    clustered_points <- cluster_ids[cluster_ids != 0]

    if (length(clustered_points) == 0) {

      largest_cluster_share <- 0
      n_clusters <- 0

    } else {

      cluster_sizes <- table(clustered_points)

      largest_cluster_share <- max(cluster_sizes) / n_points
      n_clusters <- length(cluster_sizes)

    }

    score <- largest_cluster_share + noise_share

    tibble(
      category = cat,
      eps = eps_value,
      minPts = minPts_value,
      n_points = n_points,
      n_clusters = n_clusters,
      n_noise = n_noise,
      largest_cluster_share = largest_cluster_share,
      noise_share = noise_share,
      score = score
    )
  })
})


# ------------------------------------------------------------
# All tuning results
# ------------------------------------------------------------

dbscan_tuning_results <- dbscan_tuning_results %>%
  arrange(category, score)

dbscan_best_parameters <- dbscan_tuning_results %>%
  group_by(category) %>%
  arrange(score, .by_group = TRUE) %>%
  slice(1) %>%
  ungroup()

dbscan_best_parameters
## # A tibble: 4 × 9
##   category      eps minPts n_points n_clusters n_noise largest_cluster_share
##   <chr>       <dbl>  <dbl>    <int>      <dbl>   <int>                 <dbl>
## 1 cafe          200      5     1416         56     471                0.133 
## 2 pharmacy      800      5      628         18      81                0.339 
## 3 restaurant    400      5     2791         94     330                0.325 
## 4 supermarket   700      5     2391         67      67                0.0937
## # ℹ 2 more variables: noise_share <dbl>, score <dbl>
dbscan_eps <- setNames(
  dbscan_best_parameters$eps,
  dbscan_best_parameters$category
)

dbscan_minPts <- setNames(
  dbscan_best_parameters$minPts,
  dbscan_best_parameters$category
)

dbscan_eps
##        cafe    pharmacy  restaurant supermarket 
##         200         800         400         700
dbscan_minPts
##        cafe    pharmacy  restaurant supermarket 
##           5           5           5           5
for (cat in poi_levels) {

  points_cat <- dbscan_points_all %>%
    filter(category == cat)

  p <- ggplot() +

    geom_sf(
      data = warsaw_boundary,
      fill = "grey97",
      color = "grey40",
      linewidth = 0.3
    ) +

    # noise
    geom_sf(
      data = points_cat %>%
        filter(dbscan_cluster == "0"),
      color = "grey70",
      alpha = 0.25,
      size = 0.35
    ) +

    # clustered points
    geom_sf(
      data = points_cat %>%
        filter(dbscan_cluster != "0"),
      aes(color = dbscan_cluster),
      alpha = 0.65,
      size = 1.0
    ) +

    theme_minimal() +

    labs(
      title = paste0("DBSCAN: ", cat),
      subtitle = paste0(
        "eps = ", dbscan_eps[[cat]],
        " m; minPts = ", dbscan_minPts[[cat]],
        " | yellow = cluster centres; grey = noise"
      ),
      x = NULL,
      y = NULL
    ) +

    guides(color = "none")

  print(p)
}

On the restaurant map, Several clusters correspond spatially to the city centre, district centres, shopping malls and the airport. These features are also present on the café map, although with different emphasis, as the city centre is more dominant relative to most district centres. Pharmacy clusters more closely resemble the city’s district structure, although there is still one large central cluster. The restaurant clusters show the greatest spatial variation: some districts are visible, but there are also many points corresponding to concentrations of shops in particular locations, such as areas along major roads leading out of the city.

8.3 QDC (Quick Density Clustering)

QDC classifies points according to their local spatial density. For each point, it uses two measures: (1) the total distance to the k nearest neighbours and (2) the number of neighbours within a fixed radius eps. Both variables are standardised and then clustered using k-means into three density classes. Dense points tend to have shorter distances to their neighbours and more neighbours within the specified radius, whereas sparse points show the opposite pattern.

qdc_sample_size <- 5000
qdc_nclust <- 3
qdc_k <- 10
qdc_eps <- 0.05

run_qdc_manual <- function(points_sf, sample_size = 5000, nclust = 3, k = 10, eps = 0.05) {

  points_geo <- st_transform(points_sf, crs_geo)
  n_available <- nrow(points_geo)

  if (n_available < (k + 2)) {
    stop("Too few points for QDC with k = ", k, ".")
  }

  n_use <- min(sample_size, n_available)
  selected <- sample(seq_len(n_available), n_use, replace = FALSE)
  selected_sf <- points_geo[selected, ]

  dane <- st_coordinates(selected_sf) %>%
    as.data.frame()
  names(dane) <- c("x", "y")

  # QDC spatial variable 1: sum of distances to k nearest neighbours
  knn.dist <- kNNdist(as.matrix(dane[, c("x", "y")]), k, all = TRUE)
  dane$knndist1 <- apply(knn.dist, 1, sum)

  # QDC spatial variable 2: number of neighbours inside fixed radius eps
  agg.radius <- frNN(as.matrix(dane[, c("x", "y")]), eps = eps)
  dane$frnn1 <- vapply(agg.radius$id, length, integer(1))

  # Normalisation of both spatial variables
  dane$knndist1.scaled <- as.numeric(scale(dane$knndist1))
  dane$frnn1.scaled <- as.numeric(scale(dane$frnn1))

  # K-means on the two normalised spatial variables
  km.set1 <- kmeans(dane[, c("knndist1.scaled", "frnn1.scaled")],
                    centers = nclust, nstart = 30)
  dane$km.set1 <- factor(km.set1$cluster)

  # Thresholds following the course implementation for 3 clusters
  if (nclust != 3) {
    stop("This outcome labelling follows the course implementation and requires nclust = 3.")
  }

  t1 <- max(sapply(levels(dane$km.set1), function(cl) {
    min(dane$knndist1.scaled[dane$km.set1 == cl])
  }))

  t2 <- max(sapply(levels(dane$km.set1), function(cl) {
    min(dane$frnn1.scaled[dane$km.set1 == cl])
  }))

  # High total kNN distance -> low density; high neighbour count -> high density
  dane$outcome <- ifelse(
    dane$knndist1.scaled > t1,
    "low-density",
    ifelse(dane$frnn1.scaled > t2, "high-density", "mid-density")
  )

  selected_sf$qdc_cluster <- factor(dane$outcome,
                                    levels = c("low-density", "mid-density", "high-density"))

  list(
    sampled_points = selected_sf,
    diagnostics = dane,
    thresholds = tibble(t1 = t1, t2 = t2),
    sample_size = n_use
  )
}

qdc_by_category <- setNames(
  map(poi_levels, function(cat) {
    run_qdc_manual(
      pois_sf %>% filter(category == cat),
      sample_size = qdc_sample_size,
      nclust = qdc_nclust,
      k = qdc_k,
      eps = qdc_eps
    )
  }),
  poi_levels
)

qdc_summary <- map_dfr(poi_levels, function(cat) {
  res <- qdc_by_category[[cat]]
  res$sampled_points %>%
    st_drop_geometry() %>%
    count(qdc_cluster, name = "n_points") %>%
    mutate(
      category = cat,
      sample_size = res$sample_size,
      .before = 1
    )
})

qdc_summary
##       category sample_size  qdc_cluster n_points
## 1   restaurant        2791  low-density      131
## 2   restaurant        2791  mid-density     1218
## 3   restaurant        2791 high-density     1442
## 4         cafe        1416  low-density      123
## 5         cafe        1416  mid-density      449
## 6         cafe        1416 high-density      844
## 7     pharmacy         628  low-density       82
## 8     pharmacy         628  mid-density      306
## 9     pharmacy         628 high-density      240
## 10 supermarket        2391  low-density      320
## 11 supermarket        2391  mid-density      755
## 12 supermarket        2391 high-density     1316
for (cat in poi_levels) {
  qdc_pts <- qdc_by_category[[cat]]$sampled_points

  p <- ggplot() +
    geom_sf(data = warsaw_boundary, fill = "grey97", color = "grey40", linewidth = 0.3) +
    geom_sf(data = qdc_pts, aes(color = qdc_cluster), alpha = 0.7, size = 0.65) +
    theme_minimal() +
    labs(title = paste0("QDC density classes: ", cat),
         subtitle = paste0("Sample = ", qdc_by_category[[cat]]$sample_size,
                           "; k = ", qdc_k, "; eps = ", qdc_eps,
                           "; low / mid / high density"),
         color = "QDC class", x = NULL, y = NULL)
  print(p)
}

The supermarket distribution is the most interesting, as the city centre is classified as a low-density area relative to the surrounding supermarket concentrations. Around it, there is a ring of areas with a high concentration of supermarkets, while on the outskirts there are mid-density clusters that appear to correspond to local service centres. The other POI categories have a similar distribution between each other, although pharmacies appear to be more closely aligned with the population distribution than cafés and restaurants.

9. KDE surfaces

Kernel Density Estimation turns the four discrete POI point patterns into continuous local-intensity surfaces.

kde_bandwidth <- 1000    # metres
kde_eps <- 100          # pixel size, metres

kde_to_df <- function(kde_image, value_name = "density") {
  as.data.frame(kde_image) %>% rename("{value_name}" := value)
}

kde_by_category <- setNames(map(poi_levels, function(cat) {
  pts_sf <- pois_sf %>% filter(category == cat)
  ppp_obj <- make_ppp(pts_sf, warsaw_window)
  density.ppp(ppp_obj, sigma = kde_bandwidth, kernel = "gaussian",
             edge = TRUE, eps = kde_eps)
}), poi_levels)

9.1 Raw intensity surfaces

kde_raw_df <- map_dfr(poi_levels, function(cat) {
  kde_to_df(kde_by_category[[cat]], "density_m2") %>%
    mutate(density_km2 = density_m2 * 1e6, category = cat)
})

for (cat in poi_levels) {
  p <- ggplot() +
    geom_raster(data = kde_raw_df %>% filter(category == cat),
                aes(x = x, y = y, fill = density_km2)) +
    geom_sf(data = warsaw_boundary, fill = NA, color = "grey15", linewidth = 0.3,
            inherit.aes = FALSE) +
    coord_sf(crs = st_crs(crs_pl), expand = FALSE) +
    scale_fill_viridis_c(option = "C", trans = "sqrt", name = "Points\nper km2") +
    theme_minimal() +
    labs(title = paste0("KDE intensity surface: ", cat),
         subtitle = paste0("Bandwidth = ", kde_bandwidth, " m; Gaussian kernel; eps = ", kde_eps, " m"),
         x = NULL, y = NULL)
  print(p)
}

9.2 Standardised surfaces (density per point)

Because category sizes differ, raw intensity partly reflects how many POIs exist in each category. Dividing each surface by its own point count lets us compare the spatial shape rather than the overall volume.

poi_n_lookup <- setNames(poi_counts$n_points, as.character(poi_counts$category))

kde_std_df <- map_dfr(poi_levels, function(cat) {
  img <- kde_by_category[[cat]]
  std_img <- eval.im(img / poi_n_lookup[[cat]])
  kde_to_df(std_img, "density_per_point") %>%
    mutate(category = cat)
})

for (cat in poi_levels) {
  p <- ggplot() +
    geom_raster(data = kde_std_df %>% filter(category == cat),
                aes(x = x, y = y, fill = density_per_point)) +
    geom_sf(data = warsaw_boundary, fill = NA, color = "grey15", linewidth = 0.3,
            inherit.aes = FALSE) +
    coord_sf(crs = st_crs(crs_pl), expand = FALSE) +
    scale_fill_viridis_c(option = "C") +
    theme_minimal() +
    labs(title = paste0("Standardised KDE surface: ", cat),
         subtitle = "Density divided by the category point count",
         fill = "Relative
intensity", x = NULL, y = NULL)
  print(p)
}

Restaurants and cafés have a centralised distribution. The distribution of pharmacies is broadly consistent with the distribution of the population. Supermarkets, in contrast, have many locations distributed throughout the city.

9.3 Bandwidth sensitivity (illustration)

A bandwidth of 1,000 m was selected for the subsequent analyses because it provided a useful balance between local detail and the broader spatial pattern in the bandwidth sensitivity analysis.

restaurant_ppp <- make_ppp(pois_sf %>% filter(category == "restaurant"), warsaw_window)

bandwidths <- c(200, 500, 1000, 1500, 2000)

bandwidth_df <- map_dfr(bandwidths, function(bw) {
  img <- density.ppp(restaurant_ppp, sigma = bw, kernel = "gaussian",
                     edge = TRUE, eps = kde_eps)
  kde_to_df(img, "density_m2") %>%
    mutate(density_km2 = density_m2 * 1e6,
           bandwidth = paste0(bw, " m"))
})

for (bw in unique(bandwidth_df$bandwidth)) {
  p <- ggplot() +
    geom_raster(data = bandwidth_df %>% filter(bandwidth == bw),
                aes(x = x, y = y, fill = density_km2)) +
    geom_sf(data = warsaw_boundary, fill = NA, color = "grey15", linewidth = 0.3,
            inherit.aes = FALSE) +
    coord_sf(crs = st_crs(crs_pl), expand = FALSE) +
    scale_fill_viridis_c(option = "C", trans = "sqrt") +
    theme_minimal() +
    labs(title = paste0("Restaurant KDE: bandwidth = ", bw),
         subtitle = "Bandwidth sensitivity",
         fill = "Points
per km2", x = NULL, y = NULL)
  print(p)
}

10. Degree of agglomeration: relative entropy (ETA)

The ETA idea: build a Voronoi tessellation from the points, treat tile areas as probabilities, and compute Shannon entropy. Perfectly uniform points give tiles of similar size and high (relative) entropy close to 1; strongly clustered points give very unequal tile sizes and low relative entropy.

compute_relative_entropy <- function(points_sf, boundary_sf) {

  points_union   <- st_union(st_geometry(points_sf))
  boundary_union <- st_union(st_geometry(boundary_sf))

  tess      <- st_voronoi(points_union, boundary_union)
  tess_clip <- st_intersection(st_cast(tess), boundary_union)

  areas  <- st_area(tess_clip)
  shares <- as.numeric(areas / sum(areas))
  shares <- shares[shares > 0]

  shannon_entropy <- -sum(shares * log(shares))
  n_tiles         <- length(shares)
  max_entropy     <- log(n_tiles)

  tibble(n_points = n_tiles,
         shannon_entropy = shannon_entropy,
         max_entropy = max_entropy,
         relative_entropy = shannon_entropy / max_entropy)
}
eta_sample_size <- 3000

eta_results <- pois_sf %>%
  st_drop_geometry() %>%
  distinct(category) %>%
  pull(category) %>%
  map_dfr(function(cat) {
    pts <- pois_sf %>% filter(category == cat)
    if (nrow(pts) > eta_sample_size) {
      pts <- pts %>% slice_sample(n = eta_sample_size)
    }
    compute_relative_entropy(pts, warsaw_boundary) %>%
      mutate(category = cat, .before = 1)
  })

eta_results %>% arrange(relative_entropy)
## # A tibble: 4 × 5
##   category    n_points shannon_entropy max_entropy relative_entropy
##   <fct>          <int>           <dbl>       <dbl>            <dbl>
## 1 cafe            1416            5.49        7.26            0.756
## 2 supermarket     2391            6.29        7.78            0.809
## 3 restaurant      2782            6.49        7.93            0.819
## 4 pharmacy         628            5.58        6.44            0.867
ggplot(eta_results, aes(x = reorder(category, relative_entropy), y = relative_entropy)) +
  geom_col(fill = "#2B6CB0") +
  geom_text(aes(label = round(relative_entropy, 2)), vjust = -0.4) +
  coord_cartesian(ylim = c(0, 1)) +
  theme_minimal() +
  labs(title = "Relative entropy (ETA) by POI category",
       subtitle = "Lower value = stronger spatial agglomeration",
       x = NULL, y = "Relative entropy")

Cafés have the lowest relative entropy and therefore the strongest spatial agglomeration, while pharmacies have the highest relative entropy and the most even spatial distribution. This is consistent with what was observed earlier and with the intuition derived from the previous results.

11. POI density vs. population density

11.1 Joining each POI type to the 500 m population grid

poi_grid_long <- st_join(
  pois_sf,
  population_grid %>% select(grid_id),
  join = st_within
) %>%
  st_drop_geometry() %>%
  filter(!is.na(grid_id)) %>%
  count(grid_id, category, name = "n_poi")

grid_categories <- population_grid %>%
  select(grid_id, population, cell_area_km2, pop_density_km2) %>%
  st_drop_geometry() %>%
  tidyr::crossing(category = poi_levels) %>%
  left_join(poi_grid_long, by = c("grid_id", "category")) %>%
  mutate(n_poi = replace_na(n_poi, 0),
         poi_density_km2 = n_poi / cell_area_km2)

grid_categories_sf <- population_grid %>%
  select(grid_id, geometry) %>%
  left_join(grid_categories, by = "grid_id")

11.2 Separate population-vs-POI maps

for (cat in poi_levels) {
  grid_cat <- grid_categories_sf %>% filter(category == cat)

  p_population <- ggplot() +
    geom_sf(data = grid_cat, aes(fill = pop_density_km2), color = NA) +
    geom_sf(data = warsaw_boundary, fill = NA, color = "grey20", linewidth = 0.3) +
    scale_fill_viridis_c(option = "C", trans = "sqrt", name = "Residents\nper km2") +
    theme_minimal() +
    labs(title = "Population density", x = NULL, y = NULL)

  p_poi <- ggplot() +
    geom_sf(data = grid_cat, aes(fill = poi_density_km2), color = NA) +
    geom_sf(data = warsaw_boundary, fill = NA, color = "grey20", linewidth = 0.3) +
    scale_fill_viridis_c(option = "C", trans = "sqrt", name = paste0(cat, "\nper km2")) +
    theme_minimal() +
    labs(title = paste0(cat, " density"), x = NULL, y = NULL)

  gridExtra::grid.arrange(
    p_population, p_poi, ncol = 2,
    top = paste0("500 m grid: population vs. ", cat)
  )
}

11.3 Mismatch maps: high population, low POI density

For every category, we standardise population density and the category’s POI density separately and calculate population z - POI z. Positive values identify cells with relatively high population density but relatively low density of the given POI type.

grid_mismatch <- grid_categories %>%
  group_by(category) %>%
  mutate(
    pop_z = as.numeric(scale(pop_density_km2)),
    poi_z = as.numeric(scale(poi_density_km2)),
    mismatch = pop_z - poi_z
  ) %>%
  ungroup()

for (cat in poi_levels) {
  p <- ggplot(grid_categories_sf %>%
                filter(category == cat) %>%
                select(grid_id, category, pop_density_km2, poi_density_km2) %>%
                left_join(grid_mismatch %>%
                            filter(category == cat) %>%
                            select(grid_id, mismatch), by = "grid_id")) +
    geom_sf(aes(fill = mismatch), color = NA) +
    geom_sf(data = warsaw_boundary, fill = NA, color = "grey20", linewidth = 0.3) +
    scale_fill_gradient2(low = "#2B6CB0", mid = "grey95", high = "#B83227",
                         midpoint = 0, name = "Pop. z -\nPOI z") +
    theme_minimal() +
    labs(title = paste0("Population vs. ", cat, " supply mismatch"),
         subtitle = "Positive = relatively high population and relatively low POI density",
         x = NULL, y = NULL)
  print(p)
}

# Top mismatch cells for every category
top_mismatch <- grid_mismatch %>%
  filter(population > 0) %>%
  group_by(category) %>%
  arrange(desc(mismatch), .by_group = TRUE) %>%
  slice_head(n = 10) %>%
  ungroup() %>%
  select(category, grid_id, population, n_poi, pop_density_km2, poi_density_km2, mismatch)

top_mismatch
## # A tibble: 40 × 7
##    category grid_id population n_poi pop_density_km2 poi_density_km2 mismatch
##    <chr>      <int>      <dbl> <int>           <dbl>           <dbl>    <dbl>
##  1 cafe        1012       7291     0          29191.            0        5.13
##  2 cafe         320       5746     0          23006.            0        3.95
##  3 cafe        1079       5692     0          22789.            0        3.91
##  4 cafe         201       5626     0          22526.            0        3.86
##  5 cafe         300       5550     0          22222.            0        3.80
##  6 cafe         653       5917     1          23691.            4.00     3.77
##  7 cafe         980       5503     0          22033.            0        3.76
##  8 cafe         285       5497     0          22009.            0        3.76
##  9 cafe         895       5423     0          21712.            0        3.70
## 10 cafe         377       5765     1          23082.            4.00     3.66
## # ℹ 30 more rows

11.4 Weighted KDE: population vs. each POI type

# Build population centroids from the 500 m grid. Population counts are
# used as weights for the continuous population surface.
pop_centroids <- st_centroid(population_grid) %>%
  select(population, geometry)

# Population centroids inside Warsaw
warsaw_sf <- sf::st_as_sf(warsaw_window) %>%
  sf::st_set_crs(sf::st_crs(pop_centroids))

pop_centroids_inside <- sf::st_filter(pop_centroids, warsaw_sf)

pop_ppp <- make_ppp(pop_centroids_inside, warsaw_window)
marks(pop_ppp) <- pop_centroids_inside$population

pop_kde <- density.ppp(
  pop_ppp,
  sigma = kde_bandwidth,
  weights = marks(pop_ppp),
  kernel = "gaussian",
  edge = TRUE,
  eps = kde_eps
)

pop_kde_share <- eval.im(pop_kde / sum(as.matrix(pop_kde), na.rm = TRUE))

for (cat in poi_levels) {
  poi_ppp <- make_ppp(pois_sf %>% filter(category == cat), warsaw_window)
  poi_kde <- density.ppp(poi_ppp, sigma = kde_bandwidth, kernel = "gaussian",
                         edge = TRUE, eps = kde_eps)
  poi_kde_share <- eval.im(poi_kde / sum(as.matrix(poi_kde), na.rm = TRUE))

  kde_diff_image <- eval.im(poi_kde_share - pop_kde_share)
  kde_diff_df <- kde_to_df(kde_diff_image, "difference")

  p <- ggplot(kde_diff_df, aes(x = x, y = y, fill = difference)) +
    geom_raster() +
    geom_sf(data = warsaw_boundary, fill = NA, color = "black", inherit.aes = FALSE) +
    coord_sf(crs = st_crs(crs_pl), expand = FALSE) +
    scale_fill_gradient2(low = "#2B6CB0", mid = "grey95", high = "#B83227",
                         midpoint = 0, name = "POI - pop.\nshare") +
    theme_minimal() +
    labs(title = paste0("KDE difference: ", cat, " vs. population"),
         subtitle = "Positive = POI spatial share exceeds population spatial share",
         x = NULL, y = NULL)
  print(p)
}

For supermarkets, there are areas in peripheral districts with more supermarkets than would be expected based on population. The opposite is true in the city centre. In the case of pharmacies, apart from two anomalies, there are more pharmacies relative to the population, while elsewhere the two distributions broadly coincide. For restaurants and cafés, an overrepresentation is visible in the city centre.

12. Do POI categories share the same spatial pattern?

Each standardised KDE surface (Section 9.2) is read pixel by pixel into a vector. Because all four surfaces share the same window, bandwidth, kernel and pixel size, pixel i refers to the same location in every category, so the vectors are directly comparable. Hierarchical clustering of these vectors shows which POI types tend to peak in the same parts of the city.

kde_matrix_to_row_vector <- function(kde_image) {
  v <- as.vector(t(as.matrix(kde_image)))
  v[!is.finite(v)] <- 0
  v
}

normalise_spatial_series <- function(x) {
  x_shifted <- x - min(x, na.rm = TRUE)
  x_shifted / sum(x_shifted, na.rm = TRUE)
}

spatial_series_list <- map(kde_by_category, function(img) {
  std_img <- eval.im(img / sum(as.matrix(img), na.rm = TRUE))
  normalise_spatial_series(kde_matrix_to_row_vector(std_img))
})

spatial_series_df <- as.data.frame(spatial_series_list)

spatial_series_distance <- dist(t(spatial_series_df), method = "euclidean")
spatial_series_hclust   <- hclust(spatial_series_distance, method = "complete")

ggdendro::ggdendrogram(spatial_series_hclust, rotate = TRUE) +
  labs(title = "Hierarchical clustering of POI KDE surfaces",
       subtitle = "POI types close together have similarly-shaped spatial distributions")

Restaurants and pharmacies form a similar spatial pattern, possibly because their main areas of concentration have a comparable spatial extent. The café distribution appears to be more concentrated around a smaller central area. Supermarkets form a separate group, reflecting their distinctly different spatial pattern.

13. Discussion and conclusions

Restaurants and cafés show the strongest centralisation, with the largest concentrations located in central Warsaw. This area may represent an important entertainment and leisure centre of the city, although further analysis of other relevant POI categories, such as bars, would be needed to support this interpretation.

Supermarkets show a more dispersed pattern, with relatively strong concentrations in peripheral areas and a weaker presence in the city centre. One possible explanation is that supermarkets may serve both local residents and commuters, particularly when located along major travel routes. This interpretation should be treated as a hypothesis rather than a direct finding of the analysis.

Pharmacies show the most even spatial distribution among the four categories and appear to be most closely aligned with population distribution. Nevertheless, their concentration in central Warsaw remains higher than would be expected from population density alone. Across the different methods, the results provide a broadly consistent picture: restaurants and cafés are strongly centralised, pharmacies are more evenly distributed, and supermarkets follow a distinct and more dispersed spatial pattern.

Some local anomalies may partly result from shopping centres, where multiple POIs of the same type can be concentrated within a relatively small area. This highlights the importance of considering the underlying structure of urban environments when interpreting purely spatial clustering results.

Limitations

  • OSM coverage is crowd-sourced and not guaranteed to be complete or consistently tagged; some POIs (especially small, independent businesses) may be missing or mis-tagged.
  • All results describe a single snapshot in time; POI locations in OSM change continuously as the underlying dataset is edited.
  • The supermarket category probably includes many food stores, many of which are not supermarkets in the strict sense.

Sources:

materials provided during classes

OpenStreetMap (OSM) data

GUS population data shp

text & code editing was assisted by ChatGPT