1 Research question

Public libraries are a basic element of urban social infrastructure. Their locations should be sufficiently dispersed to ensure access across the city, but they may also cluster in dense, central and historically developed areas. This project asks:

  1. Do public-library locations in Warsaw depart from complete spatial randomness (CSR)?
  2. Can their apparent clustering be explained by a large-scale centre–periphery trend?
  3. What does a simple inhomogeneous Poisson model reveal about the location pattern?

The units of analysis are service locations, not organisations. If several named departments share a building, they are treated as one site. This avoids coincident points, which are inappropriate for a simple point process.

2 Data and reproducibility

Locations were extracted on 31 August 2026 from OpenStreetMap through the Overpass API. Records tagged amenity=library within the administrative area of Warsaw were retained when the Polish name or operator contained a phrase characteristic of a public branch: biblioteka publiczna, wypożyczalnia, czytelnia, or dla dzieci. This operational rule excludes most university, school and specialist libraries, but it may still omit a public facility with incomplete OSM tags. The official city portal reports 202 public library facilities; it is used as a plausibility benchmark rather than as a point-level data source.

The city boundary was obtained through Nominatim/OpenStreetMap. All computations use Poland CS92 (EPSG:2180), so distances are in metres and areas in square metres. The raw extraction and boundary are bundled with this report, making the analysis reproducible without an internet connection.

raw <- read.csv("warsaw_public_libraries.csv", encoding = "UTF-8",
                check.names = FALSE)
boundary <- st_read("warsaw_boundary.geojson", quiet = TRUE) |>
  st_make_valid() |>
  st_transform(2180) |>
  summarise()

pts <- st_as_sf(raw, coords = c("longitude", "latitude"), crs = 4326,
                remove = FALSE) |>
  st_transform(2180)
pts <- pts[lengths(st_within(pts, boundary)) > 0, ]

# One point per physical site (one-metre rounding handles exact/coincident records).
xy0 <- st_coordinates(pts)
pts$site_key <- paste(round(xy0[,1]), round(xy0[,2]), sep = "_")
sites <- pts |>
  group_by(site_key) |>
  summarise(n_units = n(), names = paste(unique(name), collapse = " | "),
            do_union = FALSE) |>
  ungroup()

data.frame(
  stage = c("OSM records matching the rule", "Inside Warsaw boundary",
            "Unique service locations"),
  n = c(nrow(raw), nrow(pts), nrow(sites))
) |> knitr::kable(caption = "Cleaning and aggregation of the point data")
Cleaning and aggregation of the point data
stage n
OSM records matching the rule 206
Inside Warsaw boundary 206
Unique service locations 206

3 Exploratory map

ggplot() +
  geom_sf(data = boundary, fill = "grey97", colour = "grey45", linewidth = 0.35) +
  geom_sf(data = sites, aes(size = n_units), colour = "#087E8B", alpha = 0.78) +
  scale_size_continuous(range = c(1.5, 4.5), breaks = c(1, 2, 3),
                        name = "Units at site") +
  coord_sf(datum = NA) +
  labs(title = "Public-library service locations in Warsaw",
       subtitle = "Several units at one address are represented by one point",
       x = NULL, y = NULL,
       caption = "Data: OpenStreetMap contributors (ODbL), extraction 2026-08-31")

The map should be read together with the irregular observation window. Sparse areas near the city edge include forests, the airport and low-density neighbourhoods; a homogeneous model is therefore deliberately only a baseline.

4 Point-pattern object and descriptive statistics

W <- as.owin(st_geometry(boundary))
xy <- st_coordinates(sites)
X <- ppp(x = xy[,1], y = xy[,2], window = W, checkdup = TRUE)
unitname(X) <- c("metre", "metres")

area_km2 <- area.owin(W) / 1e6
intensity_km2 <- npoints(X) / area_km2
nnd_km <- nndist(X) / 1000

summary_table <- data.frame(
  Measure = c("Number of unique sites", "Window area (km²)",
              "Mean sites per km²", "Mean nearest-neighbour distance (km)",
              "Median nearest-neighbour distance (km)"),
  Value = c(npoints(X), round(area_km2, 1), round(intensity_km2, 3),
            round(mean(nnd_km), 3), round(median(nnd_km), 3))
)
knitr::kable(summary_table, caption = "Basic characteristics of the pattern")
Basic characteristics of the pattern
Measure Value
Number of unique sites 206.000
Window area (km²) 516.700
Mean sites per km² 0.399
Mean nearest-neighbour distance (km) 0.524
Median nearest-neighbour distance (km) 0.391 
ggplot(data.frame(distance = nnd_km), aes(distance)) +
  geom_histogram(binwidth = 0.25, boundary = 0, fill = "#087E8B", colour = "white") +
  labs(title = "Nearest-neighbour distances", x = "Distance (km)", y = "Number of sites")

5 First-order intensity

A kernel estimate provides a descriptive picture of the spatially varying intensity. The bandwidth is selected by likelihood cross-validation. Edge correction is applied by density.ppp.

bw <- bw.ppl(X)
lambda_hat <- density(X, sigma = bw, edge = TRUE, at = "pixels")
plot(lambda_hat / 1e6, main = paste0("Kernel intensity (sites per km²); bandwidth = ",
                                    round(bw/1000, 2), " km"),
     ribargs = list(las = 1))
plot(X, add = TRUE, pch = 16, cex = 0.35)

To represent the centre–periphery gradient parsimoniously, distance from the centroid of Warsaw is used as a spatial covariate. It is not interpreted causally: it is a large-scale trend control.

cent <- st_coordinates(st_centroid(boundary))[1,]
dist_centre <- as.im(function(x, y) sqrt((x-cent[1])^2 + (y-cent[2])^2)/1000, W = W)

fit0 <- ppm(X ~ 1)
fit_xy <- ppm(X ~ polynom(x, y, 2))
fit_d <- ppm(X ~ d + I(d^2), covariates = list(d = dist_centre))

models <- data.frame(
  Model = c("Homogeneous Poisson", "Quadratic spatial trend",
            "Quadratic distance-to-centre trend"),
  Parameters = c(length(coef(fit0)), length(coef(fit_xy)), length(coef(fit_d))),
  AIC = c(AIC(fit0), AIC(fit_xy), AIC(fit_d))
) |>
  mutate(Delta_AIC = round(AIC - min(AIC), 2), AIC = round(AIC, 2)) |>
  arrange(AIC)
knitr::kable(models, caption = "Comparison of Poisson point-process models")
Comparison of Poisson point-process models
Model Parameters AIC Delta_AIC
Quadratic spatial trend 6 6327.40 0.00
Quadratic distance-to-centre trend 3 6362.91 35.52
Homogeneous Poisson 1 6484.87 157.48 
dgrid <- seq(0, max(dist_centre$v, na.rm = TRUE), length.out = 200)
b <- coef(fit_d)
eta <- b[1] + b[2] * dgrid + b[3] * dgrid^2
effect <- data.frame(distance = dgrid, relative_intensity = exp(eta - eta[1]))
ggplot(effect, aes(distance, relative_intensity)) +
  geom_line(linewidth = 1, colour = "#C81D25") +
  geom_hline(yintercept = 1, linetype = 3) +
  labs(title = "Fitted centre--periphery trend",
       subtitle = "Intensity relative to the fitted value at the city centroid",
       x = "Distance from Warsaw centroid (km)", y = "Relative intensity")

6 Second-order structure

Ripley’s \(K\) function compares the observed number of neighbours within distance \(r\) with the CSR expectation \(K(r)=\pi r^2\). The transformed \(L(r)-r\) is easier to read: positive values indicate more close pairs than expected under CSR. A global simulation envelope (39 simulations) is shown for the homogeneous baseline.

E <- envelope(X, fun = Kest, nsim = 39, correction = "border",
              transform = expression(sqrt(./pi) - r), global = TRUE,
              savefuns = TRUE, verbose = FALSE)
plot(E, . - r ~ r, main = "CSR envelope for Ripley's L(r) - r",
     xlab = "r (metres)", ylab = "L(r) - r (metres)")
abline(h = 0, lty = 3)

Because CSR confounds interaction with a varying intensity, the inhomogeneous \(K\)-function is the more relevant diagnostic. The fitted quadratic spatial trend is used to estimate intensity at each observed site.

lambda_i <- predict(fit_xy, locations = X, type = "trend")
Ki <- Kinhom(X, lambda = lambda_i, correction = "border")
plot(Ki, sqrt(border/pi) - r ~ r, main = "Inhomogeneous L(r) - r",
     xlab = "r (metres)", ylab = "Linhom(r) - r (metres)")
abline(h = 0, lty = 3)

7 Model diagnostics

diag_xy <- diagnose.ppm(fit_xy, type = "pearson")

plot(diag_xy, which = "smooth", main = "Smoothed Pearson residuals")

Q <- quadrat.test(fit_xy, nx = 4, ny = 4)
data.frame(statistic = unname(Q$statistic), df = unname(Q$parameter),
           p_value = Q$p.value) |>
  knitr::kable(digits = 3, caption = "Pearson quadrat goodness-of-fit test")
Pearson quadrat goodness-of-fit test
statistic df p_value
22.783 9 0.013

8 Conclusions

The project separates two questions that are easily confused. The CSR analysis asks whether libraries look like points from a constant-rate Poisson process; the model comparison and inhomogeneous \(K\) analysis ask whether departures remain after a large-scale intensity trend is allowed.

Interpret the generated output as follows:

  • If the observed \(L(r)-r\) curve lies above the CSR envelope, there are more close pairs than a homogeneous Poisson process predicts.
  • A substantially lower AIC for either trend model shows that intensity is not spatially constant. The sign and shape of the relative-intensity plot describe the centre–periphery association.
  • If the inhomogeneous curve is much closer to zero than the CSR curve, much of the apparent clustering is attributable to first-order heterogeneity. Persistent positive departures suggest additional local clustering, for example co-location with civic centres or dense housing.

These results are descriptive, not an evaluation of accessibility. A stronger policy analysis would add population density, age structure, public-transport travel time, opening hours and capacity. The OSM selection rule and completeness are further limitations. Consequently, the report supports claims about the spatial pattern in the supplied data, not about causal determinants or equal access.

9 References

Baddeley, A., Rubak, E., & Turner, R. (2015). Spatial Point Patterns: Methodology and Applications with R. CRC Press.

Kopczewska, K. (2020). Applied Spatial Statistics and Econometrics: Data Analysis in R. Routledge.

Moraga, P. (2023). Spatial Statistics for Data Science: Theory and Practice with R. Chapman & Hall/CRC.

OpenStreetMap contributors (2026). OpenStreetMap data, ODbL 1.0. Data accessed via the Overpass API and Nominatim on 31 August 2026.

Warszawa 19115 (2026). Biblioteki Publiczne m.st. Warszawy. City information portal; used to benchmark the approximate number of facilities.

10 Reproducibility information

sessionInfo()
## R version 4.5.2 (2025-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=Polish_Poland.utf8  LC_CTYPE=Polish_Poland.utf8   
## [3] LC_MONETARY=Polish_Poland.utf8 LC_NUMERIC=C                  
## [5] LC_TIME=Polish_Poland.utf8    
## 
## time zone: Europe/Warsaw
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] spatstat.model_3.7-2   rpart_4.1.24           spatstat.explore_3.8-2
##  [4] nlme_3.1-168           spatstat.random_3.5-1  spatstat.geom_3.8-2   
##  [7] spatstat.univar_3.2-0  spatstat.data_3.1-9    ggplot2_4.0.3         
## [10] dplyr_1.2.1            sf_1.1-2              
## 
## loaded via a namespace (and not attached):
##  [1] s2_1.1.9              sass_0.4.10           generics_0.1.4       
##  [4] class_7.3-23          tensor_1.5.1          KernSmooth_2.23-26   
##  [7] lattice_0.22-7        digest_0.6.39         magrittr_2.0.4       
## [10] spatstat.utils_3.2-4  evaluate_1.0.5        grid_4.5.2           
## [13] RColorBrewer_1.1-3    fastmap_1.2.0         jsonlite_2.0.0       
## [16] Matrix_1.7-4          spatstat.sparse_3.2-0 e1071_1.7-17         
## [19] DBI_1.2.3             mgcv_1.9-3            scales_1.4.0         
## [22] jquerylib_0.1.4       abind_1.4-8           cli_3.6.6            
## [25] rlang_1.2.0           units_1.0-0           polyclip_1.10-7      
## [28] splines_4.5.2         withr_3.0.2           cachem_1.1.0         
## [31] yaml_2.3.12           otel_0.2.0            tools_4.5.2          
## [34] deldir_2.0-4          vctrs_0.7.3           R6_2.6.1             
## [37] proxy_0.4-29          lifecycle_1.0.5       classInt_0.4-11      
## [40] pkgconfig_2.0.3       pillar_1.11.1         bslib_0.9.0          
## [43] gtable_0.3.6          glue_1.8.0            Rcpp_1.1.1           
## [46] xfun_0.55             tibble_3.3.1          tidyselect_1.2.1     
## [49] rstudioapi_0.17.1     knitr_1.51            goftest_1.2-3        
## [52] farver_2.1.2          htmltools_0.5.9       labeling_0.4.3       
## [55] rmarkdown_2.31        wk_0.9.5              compiler_4.5.2       
## [58] S7_0.2.1