Spatial Analysis of Blood Test Points in Warsaw
Aleksandra Szpakowska
2025-09-04
1 Introduction
Access to diagnostic healthcare services, such as blood test collection points, is a key component of urban health equity. In cities like Warsaw—characterized by diverse population densities, urban sprawl, and private service operators—the spatial distribution of diagnostic points directly affects both service availability and public health outcomes.
This report presents a comprehensive spatial analysis of blood collection points in Warsaw, with a focus on identifying underserved areas and guiding recommendations for new locations. Using geospatial data, population grids, and a suite of analytical tools—including Voronoi diagrams, distance-based indicators, population-weighted metrics, and spatial autocorrelation statistics (Moran’s I, LISA) we aim to answer the following questions:
• Where are diagnostic collection points currently located?
• Which areas have insufficient or no coverage?
• How evenly are services distributed in relation to population?
• Where should new diagnostic points be placed to maximize accessibility?The study integrates datasets on existing diagnostic points (Alab, Diagnostyka, Synevo), population distribution on a 1×1 km grid, and spatial geometries of Warsaw. Special attention is given to equity of access—measured both in terms of geographic proximity and population service load—allowing us to pinpoint critical areas of need and generate data-driven recommendations for future expansion.
By combining visual exploration, quantitative indicators, and spatial statistical techniques, this analysis supports evidence-based planning for improving healthcare accessibility and reducing diagnostic service disparities across Warsaw.
2 Webscraping addresses
To obtain the full list of blood test collection points in Warsaw, a custom Python-based web scraping script was developed. The scraper targeted Diagnostyka, Alab and Synevo websites, specifically the pages listing collection facilities in Warsaw. Due to the dynamic structure of the pages, where addresses are loaded into a virtualized scrollable popup, the script utilized Selenium for browser automation. The scraper carefully interacted with dynamic elements, waited for content rendering, and extracted both the visible address data and its underlying geolocation or clinic metadata. Each address entry was parsed and saved for further spatial analysis. This method ensured completeness of the dataset and avoided issues arising from JavaScript-rendered content, which could not be accessed using static HTML parsers like BeautifulSoup. The resulting dataset provides the geographic foundation for spatial distribution, coverage, and accessibility analysis of medical services in Warsaw.
3 Loading and preparing the data
This section of the analysis consolidates address data for blood test collection points from three major providers: Diagnostyka, ALAB and Synevo into a unified dataset. Using the dplyr and tidyverse suite in R, addresses are harmonized and enriched with a standard suffix (“Warszawa, Polska”) to improve geocoding accuracy. The script employs the tidygeocoder package with the ArcGIS geocoding API to convert addresses into geographic coordinates (latitude and longitude). After removing incomplete geocoding results (missing coordinates), the data is transformed into a spatial object (sf class) with WGS84 CRS (EPSG:4326), enabling further spatial analysis.
The code loads two spatial layers administrative boundary of Warsaw (powiat Warszawa) from a shapefile and population grid data (tsed18337.shp).
The population layer is clipped to the Warsaw boundary using spatial intersection, and each grid cell is uniquely indexed with a GRID_ID. This prepares the base layers for next tasks such as service coverage analysis, population accessibility or Voronoi tessellation.
diag <- read_csv("punkty_diagnostyka_warszawa.csv")
alab <- read_csv("punkty_alab_warszawa.csv")
synevo <- read_csv("punkty_synevo_warszawa.csv")
all_points <- bind_rows(diag, alab, synevo) %>%
mutate(full_address = paste0(address, ", Warszawa, Polska"))
geocoded <- all_points %>%
geocode(
address = full_address,
method = "arcgis",
quiet = FALSE,
full_results = FALSE
)
geocoded_sf <- geocoded %>%
filter(!is.na(lat), !is.na(long)) %>%
st_as_sf(coords = c("long", "lat"), crs = 4326)
WAW <- st_read("data/powiaty/powiaty.shp", quiet = TRUE) %>%
filter(jpt_nazwa_ == "powiat Warszawa") %>%
st_transform(crs = 4326)
POP <- st_read("data/query/tsed18337.shp", quiet = TRUE) %>%
st_transform(crs = 4326) %>%
st_intersection(WAW)
POP <- POP %>%
mutate(GRID_ID = row_number())4 Exploring and visualising the data
To gain an initial spatial understanding of the distribution of blood test collection facilities, an interactive dot map was created using the tmap package in view mode (Leaflet-based). Each point represents a single medical collection site in Warsaw, color-coded by provider — ALAB, Diagnostyka, or Synevo. A small number of points were marked as “Missing,” likely due to incomplete or ambiguous source data.
tmap_mode("view")
tm_shape(geocoded_sf) +
tm_dots(
col = "operator",
palette = "Dark2",
size = 0.4,
title = "Operator"
)To assess the spatial distribution of blood test facilities relative to urban structure, the geocoded service points were spatially joined to a population grid covering the Warsaw area. Each grid cell represents a small territorial unit (1 km²), and the number of points falling within each cell (n_punkty) was calculated. Missing values were replaced with zeros to ensure completeness across the urban area.
joined <- st_join(geocoded_sf, POP, join = st_within, left = FALSE)
points_per_grid <- joined %>%
mutate(GRID_ID = as.character(GRID_ID)) %>%
group_by(GRID_ID) %>%
summarise(n_punkty = n(), .groups = "drop")
POP_with_counts <- POP %>%
mutate(GRID_ID = as.character(GRID_ID)) %>%
left_join(st_drop_geometry(points_per_grid), by = "GRID_ID") %>%
mutate(n_punkty = replace_na(n_punkty, 0))
tmap_mode("view")
tm_shape(POP_with_counts) +
tm_polygons("n_punkty", style = "quantile", palette = "YlOrRd", title = "number of points") +
tm_shape(geocoded_sf) +
tm_dots(col = "black", size = 0.2)The resulting choropleth map visualizes the local density of collection sites using a quantile classification:
• Darker red cells indicate higher concentrations of collection points.
• Lighter yellow areas represent lower densities or absence of facilities.
• Each individual point is also plotted in black for additional spatial clarity.without_points <- POP_with_counts %>%
filter(n_punkty == 0)
tm_shape(without_points) +
tm_polygons(col = "lightgray", border.col = "white") +
tm_shape(geocoded_sf) +
tm_dots(col = "operator", palette = "Set2", size = 0.2, title = "Operator") +
tm_layout(title = "Grid cells without collection points")This map highlights grid cells within Warsaw that lack any blood collection facilities, based on a spatial join between the geocoded locations and the population grid. Cells with zero service points (n_punkty == 0) were isolated and visualized in light gray, overlaid with the actual facility locations colored by operator.
Large continuous zones without coverage are visible, particularly in the north-eastern and south-western parts of the city. These underserved areas exist despite urbanization in some regions, suggesting a mismatch between population distribution and service accessibility. Central districts show near-complete saturation, while peripheral districts face potential service deserts.
5 Distance and density based indicators
distances <- st_nn(POP, geocoded_sf, k = 1, returnDist = TRUE)POP$nearest_dist_km <- as.numeric(distances$dist) / 1000
ggplot(POP, aes(x = nearest_dist_km)) +
geom_histogram(binwidth = 0.5, fill = "steelblue", color = "white") +
labs(
title = "Distance to nearest blood test point",
x = "Distance (km)",
y = "Number of grid cells"
) +
theme_minimal()
This histogram illustrates the distribution of distances from each population grid cell to its closest blood collection facility, based on a nearest-neighbour spatial query (st_nn() with k = 1). The distances are calculated in kilometers and grouped into 0.5 km bins.
Most of Warsaw’s population cells are located within 1-2 kilometers of the nearest facility, which indicates generally good accessibility in central and moderately urbanized areas. However, a long tail of the distribution reveals that some peripheral zones are located over 5 km from the nearest collection point, potentially posing accessibility challenges for residents without private transport or those with mobility constraints.
POP$area_km2 <- st_area(POP) / 1e6
POP$area_km2 <- as.numeric(POP$area_km2)
POP_with_counts$area_km2 <- st_area(POP_with_counts) / 1e6
POP_with_counts$area_km2 <- as.numeric(POP_with_counts$area_km2)
POP_with_counts <- POP_with_counts %>%
mutate(
density_points_km2 = n_punkty / area_km2,
density_points_per_capita = ifelse(tot > 0, n_punkty / tot, NA)
)
tm_shape(POP_with_counts) +
tm_polygons("density_points_km2", palette = "Blues", style = "jenks", title = "points per km²")This choropleth map presents the spatial density of blood collection points per square kilometer, calculated by dividing the number of points (n_punkty) within each grid cell by the cell’s surface area.
The highest densities (3+ points per km², darkest blue) are concentrated in central and western Warsaw, where service providers tend to cluster. Peripheral districts generally show low or zero density, confirming earlier findings on service gaps. The density metric highlights redundancy in central zones (overlapping service coverage), versus scarcity in outlying areas, providing evidence for targeted network optimization.
tm_shape(POP_with_counts) +
tm_polygons("density_points_per_capita", palette = "Greens", style = "jenks", title = "points per capita")This map illustrates the number of blood test collection points per resident in each grid cell, offering a population-adjusted view of service coverage. The metric was computed as the ratio of points to total population (density_points_per_capita).
Darker green areas represent higher service availability per resident, often found in low-population zones that still host facilities (e.g. commercial/office districts). Lighter areas (light green to white) signal low service-per-person ratios, highlighting where demand may exceed local supply - even if geographic density seems acceptable. Notably, some densely populated grid cells have poor per capita access, revealing hidden service imbalances not visible in maps based solely on spatial density.
6 Voronoi Analysis
Voronoi diagrams are a powerful spatial tool for estimating catchment areas of point-based services. In this context, the Voronoi method divides Warsaw into non-overlapping zones, where each area is assigned to the nearest collection point, assuming straight-line (Euclidean) access.
geocoded_m <- st_transform(geocoded_sf, 2180)
WAW_m <- st_transform(WAW, 2180)
POP_m <- st_transform(POP, 2180)
geocoded_m_unique <- geocoded_m %>%
distinct(geometry, .keep_all = TRUE) %>%
tibble::rowid_to_column("point_id")
bbox_waw <- st_bbox(st_transform(WAW, 2180)) %>% st_as_sfc()
voronoi <- st_voronoi(do.call(c, st_geometry(geocoded_m_unique)), envelope = bbox_waw) %>%
st_collection_extract("POLYGON") %>%
st_sf() %>%
st_set_crs(2180) %>%
st_intersection(st_transform(WAW, 2180))
voronoi$point_id <- geocoded_m_unique$point_id
voronoi$operator <- geocoded_m_unique$operator
voronoi_wgs <- st_transform(voronoi, 4326)
tm_shape(voronoi_wgs) +
tm_polygons(col = "operator", palette = "Set2", alpha = 0.4, title = "Catchment areas") +
tm_shape(geocoded_sf) +
tm_dots(size = 0.2, col = "black") +
tm_layout(title = "Blood test points catchment areas (Voronoi)")Each colored polygon represents the Voronoi cell of a collection point. The coloring reflects the service provider (ALAB, Diagnostyka, Synevo). Black dots mark the actual facility locations.
Voronoi zones vary significantly in size, revealing that some points serve only a hyperlocal area, while others (especially in outer districts) cover vast underserved zones. Diagnostic coverage is relatively well-balanced in central Warsaw, while peripheral areas show dominance by a single provider or large exclusive zones. This spatial segmentation can be used to detect overlapping catchments, monopolized areas, or coverage voids between facilities.
voronoi <- voronoi %>%
mutate(area_km2 = as.numeric(st_area(geometry)) / 1e6)
intersected <- st_intersection(POP_m, voronoi)
intersected <- intersected %>%
mutate(
inter_area = as.numeric(st_area(.)),
pop_weighted = tot * inter_area / as.numeric(st_area(geometry))
)
voronoi_stats <- intersected %>%
group_by(point_id, operator) %>%
summarise(
pop_sum = sum(pop_weighted, na.rm = TRUE),
area_km2 = unique(voronoi$area_km2[voronoi$point_id == cur_group()$point_id]),
.groups = "drop"
)
voronoi_stats <- voronoi_stats %>%
mutate(
pop_density = pop_sum / area_km2
)
operator_summary <- voronoi_stats %>%
group_by(operator) %>%
summarise(
total_area = sum(area_km2),
total_pop = sum(pop_sum),
n_points = n(),
avg_pop_per_point = total_pop / n_points,
avg_density = total_pop / total_area,
.groups = "drop"
)
#print(operator_summary)
operator_summary %>% arrange(desc(avg_pop_per_point))## Simple feature collection with 3 features and 6 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 626505.9 ymin: 472229.5 xmax: 655260.3 ymax: 502172.4
## Projected CRS: ETRF2000-PL / CS92
## # A tibble: 3 × 7
## operator total_area total_pop n_points avg_pop_per_point avg_density
## <chr> <dbl> <dbl> <int> <dbl> <dbl>
## 1 Alab 143. 2411657 46 52427. 16836.
## 2 Synevo 157. 602158 13 46320. 3829.
## 3 Diagnostyka 216. 1894042 44 43046. 8759.
## # ℹ 1 more variable: geometry <GEOMETRY [m]>Alab covers the highest number of people per point and operates in the densest regions, suggesting it is concentrated in urban cores. Synevo, despite fewer points, shows lowest population density - indicating its service zones are spread over larger, likely peripheral areas. Diagnostyka appears more balanced, with moderate density and population per point, reflecting a middle-ground spatial strategy.
tm_shape(voronoi_stats) +
tm_polygons("pop_sum", palette = "YlOrRd", style = "quantile", title = "Number of people assinged to each collection point") +
tm_shape(geocoded_sf) +
tm_dots(col = "black", size = 0.1)This map visualizes the total number of people assigned to each blood collection point, using Voronoi polygons to define theoretical catchment areas. The map uses a quantile-based classification to show varying population loads per point.
Dark red zones indicate facilities responsible for serving over 60,000 people, suggesting potentially overburdened locations in the eastern and southern parts of Warsaw. Lighter areas correspond to service points with relatively low population loads — often in central or over-served zones where multiple facilities compete within a small area. The uneven distribution reveals spatial imbalances in service coverage: some points operate in dense, underserved areas, while others cover smaller populations or share their catchments.
top_points <- voronoi_stats %>%
arrange(desc(pop_sum)) %>%
slice_head(n = 5)
#top_points
candidates_for_expansion <- voronoi_stats %>%
filter(pop_sum > 50000, area_km2 > 10) %>%
mutate(recommendation = "Add new diagnostic point")
#print(candidates_for_expansion)
tmap_mode("view")
tm_shape(voronoi_stats) +
tm_polygons("pop_sum", palette = "Greys", style = "quantile", title = "population") +
tm_shape(candidates_for_expansion) +
tm_borders(col = "red", lwd = 2) +
tm_text("recommendation", size = 0.8) +
tm_layout(title = "Recommenadations for new diagnostic points")This map identifies optimal zones for expanding the network of blood test collection facilities in Warsaw, based on previous Voronoi analysis and population load indicators. The areas outlined in red are Voronoi cells with the highest assigned populations (over 62,000 people per point), signaling potential service overload.
These high-priority zones are predominantly located on the urban periphery — especially in north-east, east, south, and far west Warsaw - areas that were repeatedly flagged throughout the analysis for high distance, low density, and population pressure. Each labeled cell represents a geographically underserved zone, where adding a new facility would rebalance accessibility, reduce catchment sizes, and lower individual travel burdens. This decision layer provides an evidence-based recommendation for network expansion, aligned with spatial equity and public health accessibility goals.
tmap_mode("view")
tm1 <- tm_shape(voronoi_stats) +
tm_polygons("pop_density", palette = "YlGnBu", style = "quantile", title = "Population density (people/km²)") +
tm_shape(geocoded_sf) +
tm_dots(col = "operator", palette = "Dark2", size = 0.15) +
tm_layout(title = "Population density in catchments")
tm2 <- tm_shape(voronoi_stats) +
tm_polygons("pop_sum", palette = "YlOrRd", style = "quantile", title = "Number of residents within polygons") +
tm_shape(geocoded_sf) +
tm_dots(col = "operator", palette = "Dark2", size = 0.15) +
tm_layout(title = "Population assigned to each facility")
tmap_arrange(tm1, tm2, ncol = 2)This side-by-side interactive map visualizes two complementary but distinct aspects of healthcare service accessibility within Voronoi-based catchment areas:
- Left panel — Population density in catchments
This map shows the raw population density (people/km²) inside each Voronoi zone. It helps highlight areas of urban intensity, such as central Warsaw, where land use is dense and populations are concentrated. These areas are usually expected to have higher demand for frequent and close access to medical services.
- Right panel — Population assigned to each facility
This map shows the total number of residents within each Voronoi polygon - that is, the effective population “served” by each collection point. It provides a measure of service load on each point.
candidates_for_expansion <- voronoi_stats %>%
filter(pop_sum > 50000, area_km2 > 10) %>%
mutate(
recommendation = paste0("Add new diagnostic point for operator: ", operator)
)
tmap_mode("view")
tm_shape(voronoi_stats) +
tm_polygons("pop_sum", palette = "Greys", style = "quantile", title = "Population") +
tm_shape(candidates_for_expansion) +
tm_borders(col = "red", lwd = 2) +
tm_text("operator", size = 0.7, col = "red", auto.placement = TRUE) +
tm_layout(title = "Recommended places for new diagnostic points")This map refines the location-allocation strategy by identifying priority catchment areas for expansion, with direct reference to which operator should be responsible for the new diagnostic facility.
The filtering criteria were: Voronoi zones with a population above 50,000 and area larger than 10 km², indicating both high demand and broad spatial responsibility.
Synevo is assigned two zones for expansion, both covering large eastern and southeastern regions with high population counts and low density of competing providers. Diagnostyka is recommended to expand in the south and northwest, where it currently holds large zones with sparse facilities. These recommendations are generated automatically using a reproducible filter on service burden metrics, supporting data-driven planning for infrastructure equity and efficiency.
POP_m <- st_transform(POP, 2180)
points_m <- st_transform(geocoded_sf, 2180)
POP_m$dist_to_point <- st_distance(POP_m, points_m) %>%
apply(1, min) %>%
as.numeric() / 1000
new_candidate_areas <- POP_m %>%
filter(tot > 2000, dist_to_point > 2.5) %>%
arrange(desc(tot))
tmap_mode("view")
tm_shape(POP_m) +
tm_polygons("dist_to_point", palette = "Blues", style = "quantile", title = "Distance to points (km)") +
tm_shape(new_candidate_areas) +
tm_borders(col = "red", lwd = 2) +
tm_text("pop", size = 0.7) +
tm_shape(geocoded_sf) +
tm_dots(col = "black", size = 0.1) +
tm_layout(title = "Recommended locations for a new operator")This map proposes specific grid cells where a new healthcare provider could strategically enter the Warsaw market, based on a simple, data-driven rule: population greater than 2,000 and distance greater than 2.5 km from the nearest existing collection point.
Recommended areas are predominantly located in southern, eastern, and far western parts of Warsaw, aligning with earlier Voronoi and density analyses. The selected zones are characterized by both remoteness and demographic weight, representing valuable white space for a new operator to enter without competing directly with existing networks. These recommendations offer clear spatial guidance for expansion planning, public–private partnerships, or competitive market entry.
POP_m <- POP_m %>%
mutate(
demand_score = tot * dist_to_point
)
top10 <- POP_m %>%
arrange(desc(demand_score)) %>%
slice_head(n = 10)
# top 10
tmap_mode("view")
tm_shape(POP_m) +
tm_polygons("demand_score", palette = "Greens", style = "quantile", title = "demand score") +
tm_shape(top10) +
tm_borders(col = "red", lwd = 2) +
tm_text("pop", size = 0.7, col = "red") +
tm_layout(title = "TOP 10 locations for a new operator")