Spatial machine learning for air quality mapping in Poland

Reconstruction of a pollutant concentration surface from monitoring stations and population exposure assessment

Goals of the paper:

1. Reconstruct a national concentration surface from the sparse network of monitoring stations.
2. Translate the surface into a population exposure assessment.

Data files: Chief Inspectorate of Environmental Protection (GIOŚ) hourly file for the pollutant PM2.5 2024_PM25_1g.xlsx, the station metadata workbook, and the GUS 2021 population grid population_grid_2021.shp. DEGURBA is built from GHSL via flexurba.

1 Introduction

Poor air quality is one of Poland’s most serious environmental and health problems. In winter the dominant source is low-stack emission - domestic heating with solid fuels - and long-term exposure to particulate matter is associated with thousands of premature deaths per year. State monitoring (GIOŚ) relies on a sparse network of stations, only a few per voivodeship. The question is how to reconstruct a credible, continuous concentration field over the whole country from a small number of point measurements.

The aim of this project is to apply spatial machine learning methods to reconstruct the 2024 annual concentration surface for Poland from station measurements and to assess its health consequences through a population exposure analysis.

Specific objectives:

• build a PM2.5 concentration surface from station measurements, comparing deterministic interpolation (IDW), ordinary and universal kriging, and a random-forest + kriging hybrid with spatial covariates,
• cross-validate the methods on held-out stations (RMSE, MAE) and select the best reconstruction,
• map the prediction uncertainty (kriging variance),
• assess population exposure: population-weighted mean, the share of residents above WHO/EU thresholds, and the urban-rural gap (DEGURBA) by overlaying the surface with the population grid.

2 Background and related work

This study sits at the intersection of three research streams.

Spatial machine learning

Standard machine learning models assume independent observations, which spatial data violate - nearby locations are correlated, and the relationship between predictors and outcome may vary across space.

Geostatistical interpolation and ML hybrids

Reconstructing a continuous surface from point measurements is a geostatistical problem. Ordinary and universal kriging model spatial autocorrelation directly, while regression- and random forest-assisted kriging (Hengl et al., 2018) combine a covariate model with kriging of its residuals, typically improving accuracy when auxiliary variables carry signal. Random forest plus residual kriging is exactly this hybrid.

Air-quality mapping and exposure

Estimating pollution between sparse monitors is a long-standing problem, addressed by land-use regression and geostatistical interpolation that exploit covariates such as population, traffic and land cover (Hoek et al., 2008). Turning a concentration surface into a population exposure estimate requires a population layer.

3 Packages

packages <- c(
  "tidyverse",
  "readxl",
  "sf",
  "spdep",
  "gstat",
  "automap",
  "ranger",
  "Metrics",
  "viridis",
  "scales",
  "gridExtra",
  "dplyr"
)

tibble::tibble(
  package   = packages,
  installed = vapply(packages, requireNamespace, quietly = TRUE,
                     FUN.VALUE = logical(1))
)

install.packages(packages)
install.packages(c("terra", "tidyterra", "flexurba"))

library(tidyverse)
library(readxl)
library(sf)
library(spdep)
library(gstat)
library(automap)
library(ranger)
library(Metrics)
library(viridis)
library(scales)
library(gridExtra)
library(dplyr)

# EPSG:2180 — projected CRS for Poland (metres). Correct for distances, buffers, grid sizes and kriging. EPSG:4326 (lon/lat) is only used to read raw station coordinates.
crs_pl  <- 2180
crs_geo <- 4326

4 Data

4.1 Data sources

Every layer is projected to EPSG:2180. Concentrations are in µg/m³.

Layer	Source	Type	Role in the analysis
Station measurements (hourly, 2024)	GIOŚ	points, time series	reconstruction target (aggregated to the annual/winter mean)
Station metadata & coordinates	GIOŚ	table	geolocation (joined by Kod stacji)
Population grid (1 km, Census 2021)	GUS (Statistics Poland)	polygons	prediction grid + exposure weight
Degree of Urbanisation (DEGURBA)	GHSL / Eurostat via flexurba	raster to classes	urban-rural covariate and exposure split

4.2 Station measurements

The GIOŚ hourly files are wide: column 1 holds the timestamp, every other column is one station, and the first rows are metadata (station code, indicator, averaging time, unit). Below code reads the raw sheet, recover the station codes from the “Kod stacji” row, and aggregate each station’s hourly series to a single 2024 value, keeping only stations with at least 75% of valid hours. For PM the natural unit is a daily/annual mean; aggregation = "winter_mean" restricts to the heating season (Jan–Mar, Nov–Dec).

# Read one GIOŚ wide hourly file and return per-station summary value
read_gios_hourly <- function(path, completeness_min = 0.75, months = NULL) {
  raw <- readxl::read_excel(path, col_names = FALSE, .name_repair = "minimal")

  # Row 2 is "Kod stacji"; hours start after row 6 (header)
  station_codes <- as.character(unlist(raw[2, -1], use.names = FALSE))
  body <- raw[7:nrow(raw), ]

  datetime <- as.POSIXct(as.character(unlist(body[, 1], use.names = FALSE)),
                         format = "%Y-%m-%d %H:%M", tz = "UTC")
  vals <- suppressWarnings(
    as.data.frame(lapply(body[, -1], function(col) as.numeric(as.character(col))))
  )
  names(vals) <- station_codes

  if (!is.null(months)) { # winter-mean
    keep <- as.integer(format(datetime, "%m")) %in% months
    vals <- vals[keep, , drop = FALSE]
  }

  out <- tibble::tibble(
    station_code = station_codes,
    completeness = colSums(!is.na(vals)) / nrow(vals),
    conc         = colMeans(vals, na.rm = TRUE)
  )
  out[!is.na(out$station_code) & out$completeness >= completeness_min, ]
}

winter_months <- if (params$aggregation == "winter_mean") c(1, 2, 3, 11, 12) else NULL

stations_value <- read_gios_hourly(
  params$hourly_path,
  completeness_min = params$completeness_min,
  months = winter_months
)

nrow(stations_value) # number of stations with a valid value

## [1] 95

summary(stations_value$conc) # concentration distribution [ug/m3]

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   6.772  12.954  14.930  15.214  17.212  23.275

4.3 Station coordinates

Coordinates come from the metadata workbook (sheet STACJE), joined to the measurements by Kod stacji.

meta_raw <- readxl::read_excel(params$meta_path, sheet = "STACJE")

nm      <- names(meta_raw)
lat_col <- nm[grepl("WGS84", nm) & grepl("N\\s*$", nm)][1]   # "WGS84 φ N"
lon_col <- nm[grepl("WGS84", nm) & grepl("E\\s*$", nm)][1]   # "WGS84 λ E"

meta <- meta_raw %>%
  transmute(
    station_code = `Kod stacji`,
    voivodeship  = Województwo,
    station_type = `Typ stacji`,
    lat = as.numeric(.data[[lat_col]]),
    lon = as.numeric(.data[[lon_col]])
  ) %>%
  filter(!is.na(lat), !is.na(lon))

stations_sf <- stations_value %>%
  inner_join(meta, by = "station_code") %>%
  st_as_sf(coords = c("lon", "lat"), crs = crs_geo) %>%
  st_transform(crs_pl)

nrow(stations_sf) # geolocated stations with a valid value

## [1] 95

4.4 Population grid (GUS 2021)

The GUS 2021 census grid gives residents per 1 km cell (column res). It is both the prediction target and the exposure weight. We keep populated cells and use their sum as the national frame.

pop_grid <- st_read(params$pop_grid_path, quiet = TRUE) %>%
  st_make_valid() %>%
  st_transform(crs_pl)

pop_grid$pop <- as.numeric(pop_grid[[params$pop_col]])
pop_grid <- pop_grid %>% filter(!is.na(pop), pop > 0)

poland <- st_union(pop_grid) # study area
nrow(pop_grid) # populated 1 km cells

## [1] 197130

sum(pop_grid$pop) # total population

## [1] 37019327

ggplot() +
  geom_sf(data = poland, fill = "grey97", color = "grey60", linewidth = 0.3) +
  geom_sf(data = stations_sf, aes(color = conc), size = 1.4, alpha = 0.85) +
  scale_color_viridis_c(option = "C") +
  theme_minimal() +
  labs(title = paste0(params$pollutant, " at monitoring stations - Poland ",
                      params$year),
       subtitle = paste0(nrow(stations_sf), " stations, ",
                         params$aggregation),
       color = "ug/m3")

Figure shows the input data: 95 monitoring stations with more than 75% completed data, coloured by their 2024 annual mean PM2.5, over a population grid covering about 37 million residents. Station coverage is uneven - it is denser in the south and around cities, rare across the north-east. Annual means range approximately from 8 to 22 µg/m³, with the highest values concentrated in southern Poland.

5 Spatial covariates

The reconstruction uses environmental information at each location. Below are listed three covariates that may drive pollution:

• population density (pop) - a proxy for household heating intensity,
• distance to the nearest urban centre - pollution tends to peak in and around cities,
• planar coordinates (x, y).

# The densest cells (top 1% of population)
grid_xy   <- st_coordinates(st_centroid(st_geometry(pop_grid)))
centre_th <- quantile(pop_grid$pop, 0.99, na.rm = TRUE)
centres   <- st_centroid(st_geometry(pop_grid[pop_grid$pop >= centre_th, ]))
centres   <- st_sf(geometry = centres)

# Distance to nearest centre (metres) for any set of points
dist_to_centre <- function(points_sf) {
  nf <- st_nearest_feature(points_sf, centres)
  as.numeric(st_distance(points_sf, centres[nf, ], by_element = TRUE))
}

# Station covariates: nearest populated cell giving local population, distance to centre, coordinates
st_cell <- st_nearest_feature(stations_sf, pop_grid)
stations_sf <- stations_sf %>%
  mutate(
    pop_dens    = pop_grid$pop[st_cell],
    dist_centre = dist_to_centre(stations_sf),
    x = st_coordinates(.)[, 1],
    y = st_coordinates(.)[, 2]
  )

# Grid covariates (the prediction targets)
grid_sf <- st_centroid(pop_grid) %>%
  mutate(
    pop_dens    = pop_grid$pop,
    dist_centre = dist_to_centre(.),
    x = grid_xy[, 1],
    y = grid_xy[, 2]
  )

The “distance to centre” feature depends on how centres are defined - here it is a density threshold, in next chapter it is the DEGURBA urban-centre class.

6 Urban structure with DEGURBA

DEGURBA encodes the urban–rural gradient that organises heating emissions and is the basis for the exposure breakdown. Workflow downloads the GHSL layers (built-up, population, land) and classifies the grid.

library(terra); library(tidyterra); library(flexurba)

# 1. Download global GHSL layers
download_GHSLdata(output_directory = "data/global",
                  filenames = c("built.tif", "pop.tif", "land.tif"))

# 2. Crop to Poland
poland_moll <- st_transform(poland, "ESRI:54009")
crop_GHSLdata(extent = ext(vect(poland_moll)),
              global_directory = "data/global",
              global_filenames = c("built.tif", "pop.tif", "land.tif"),
              output_directory = "data/poland",
              output_filenames = c("built_pl.tif", "pop_pl.tif", "land_pl.tif"))

# 3. Preprocess and classify (Level 1: 1 rural, 2 cluster, 3 centre)
data_pl    <- DoU_preprocess_grid(directory = "data/poland",
                                  c("built_pl.tif", "pop_pl.tif", "land_pl.tif"))
degurba_l1 <- DoU_classify_grid(data = data_pl)
terra::writeRaster(degurba_l1, "data/degurba_pl.tif", overwrite = TRUE)

# Attach DEGURBA class to each grid cell and station. If the raster exists use the real classification; otherwise fall back to population tertiles
deg_labels <- c("1" = "rural", "2" = "cluster", "3" = "centre")

if (file.exists("data/degurba_pl.tif") && requireNamespace("terra", quietly = TRUE)) {
  deg <- terra::rast("data/degurba_pl.tif")
  to_deg <- function(points_sf) {
    v <- terra::vect(st_transform(points_sf, terra::crs(deg)))
    factor(deg_labels[as.character(terra::extract(deg, v)[, 2])],
           levels = deg_labels)
  }
  grid_sf$degurba     <- to_deg(grid_sf)
  stations_sf$degurba <- to_deg(stations_sf)
} else {
  message("DEGURBA raster not found — using population-tertile")
  brks <- quantile(pop_grid$pop, c(0, .5, .9, 1), na.rm = TRUE)
  grid_sf$degurba     <- cut(grid_sf$pop_dens, brks,
                             labels = deg_labels, include.lowest = TRUE)
  stations_sf$degurba <- cut(stations_sf$pop_dens, brks,
                             labels = deg_labels, include.lowest = TRUE)
}

table(grid_sf$degurba)

## 
##   rural cluster  centre 
##  183906   10319    2779

7 Surface reconstruction

The concentration field on the population-grid cells is reconstructed in three ways:

1. Ordinary and universal kriging - model the residual spatial autocorrelation, universal kriging adding a covariate trend (population density, distance to the nearest centre),
2. RF + kriging - a random forest learns the covariate signal, then the residual spatial structure is kriged and added back.

The surface and uncertainty maps shown below are produced with the RF-kriging hybrid. All methods are then benchmarked against held-out stations by cross-validation, which selects the best reconstruction.

7.1 Train/test split and cross-validation

set.seed(2026)
stations_sf$u <- runif(nrow(stations_sf))
train_sf <- stations_sf[stations_sf$u <  0.8, ]
test_sf  <- stations_sf[stations_sf$u >= 0.8, ]
c(train = nrow(train_sf), test = nrow(test_sf))

## train  test 
##    77    18

A single random split is optimistic for spatial data, because nearby stations end up in both train and test. With 95 stations it is an acceptable baseline, but the ranking below is confirmed with ten-fold cross-validation.

7.2 Variogram

plot(autofitVariogram(conc ~ 1, train_sf)) # empirical + fitted variogram

A large nugget and a weak rise indicate limited residual spatial structure - which justifies the methods below.

7.3 Universal kriging

uk <- autoKrige(conc ~ pop_dens + dist_centre, train_sf, grid_sf, verbose = FALSE)

## [using universal kriging]

uk_out <- st_as_sf(uk$krige_output)

pop_grid$pred_uk    <- uk_out$var1.pred
pop_grid$pred_uk_sd <- sqrt(uk_out$var1.var)

g_pred_uk <- ggplot(pop_grid) +
  geom_sf(aes(fill = pred_uk), color = NA) +
  scale_fill_viridis_c(option = "C") + theme_minimal() +
  labs(title = "Reconstructed PM2.5 surface (universal kriging)", fill = "ug/m3")

g_unc_uk <- ggplot(pop_grid) +
  geom_sf(aes(fill = pred_uk_sd), color = NA) +
  scale_fill_viridis_c(option = "D") + theme_minimal() +
  labs(title = "Prediction uncertainty (kriging SD)", fill = "ug/m3")

grid.arrange(g_pred_uk, g_unc_uk, ncol = 2)

Reconstructed surface (left): the universal kriging reconstruction reproduces the expected Polish PM2.5 pattern: annual means are highest across central–southern Poland (Silesia, the Kraków and Łódź regions, roughly 18–20°E and 50–52°N) and lowest across the north and north-east, ranging from about 12 to 18 µg/m³. Because universal kriging models a covariate trend (population density, distance to the nearest urban centre) on top of the spatial structure, the surface captures the south to north gradient that drives PM2.5 - consistent with the dominant role of residential heating and urban structure.

Prediction uncertainty (right): the map shows the kriging standard deviation (about 2.3–3.0 µg/m³). In contrast to the concentration surface, uncertainty follows the monitoring network rather than the pollution level: it is lowest at and around stations and rises over gaps between them, so the largest values appear where stations are rare. The map indicates where the interpolated surface is least reliable and where additional monitoring would most reduce uncertainty.

Together the panels support the modelling choice. For PM2.5 at the national scale the signal is carried mainly by environmental covariates (urban–rural context, population), which universal kriging captures explicitly as a trend; the residual spatial structure is weak (consistent with the near-flat variogram), so kriging mainly interpolates smoothly between stations. Cross-validation (below) confirms universal kriging as the most accurate of the compared methods.

7.4 RF-assisted kriging

# 1) RF on covariates
rf <- ranger(conc ~ pop_dens + dist_centre + x + y,
             data = st_drop_geometry(train_sf),
             num.trees = 500, set.seed(2026), keep.inbag = TRUE)

# 2) RF trend + RF standard error at stations and on the grid
train_sf$rf_fit <- predict(rf, st_drop_geometry(train_sf))$predictions
grid_pred       <- predict(rf, st_drop_geometry(grid_sf), type = "se")
grid_sf$rf_fit  <- grid_pred$predictions
grid_sf$rf_se   <- grid_pred$se

# 3) Krige the residuals and add them back
train_sf$resid <- train_sf$conc - train_sf$rf_fit
res_k   <- autoKrige(resid ~ 1, train_sf, grid_sf, verbose = FALSE)

## [using ordinary kriging]

res_out <- st_as_sf(res_k$krige_output)

grid_sf$pred    <- grid_sf$rf_fit + res_out$var1.pred # reconstructed surface
grid_sf$pred_sd <- sqrt(grid_sf$rf_se^2 + res_out$var1.var) # total SD (RF + kriging)

pop_grid$pred    <- grid_sf$pred
pop_grid$pred_sd <- grid_sf$pred_sd

g_pred <- ggplot(pop_grid) +
  geom_sf(aes(fill = pred), color = NA) +
  scale_fill_viridis_c(option = "C") + theme_minimal() +
  labs(title = "Reconstructed PM2.5 surface (RF-kriging)", fill = "ug/m3")

g_unc <- ggplot(pop_grid) +
  geom_sf(aes(fill = pred_sd), color = NA) +
  scale_fill_viridis_c(option = "D") + theme_minimal() +
  labs(title = "Prediction uncertainty (total SD: RF + kriging)", fill = "ug/m3")

grid.arrange(g_pred, g_unc, ncol = 2)

Reconstructed surface (left): the RF-kriging reconstruction reproduces the expected Polish PM2.5 pattern: annual means are highest across central–southern part of Poland (Silesia, Kraków and Łódź regions, roughly 18–20°E and 50–52°N) and around major agglomerations, and lowest across the north and north-east, ranging from about 12 to 18 µg/m³. Because the random forest exploits the covariates (population density, distance to the nearest urban centre), the surface captures city-scale hotspots that ordinary kriging alone would miss - consistent with the dominant role of residential heating and urban structure.

Prediction uncertainty (right): the map shows the total predictive standard deviation, combining the random forest error and the residual-kriging variance (about 1.4–2.1 µg/m³). Uncertainty is highest in the high-concentration, high-gradient central–southern part of country and near the southern border, and lowest over the more homogeneous north and north-east. It therefore follows the variability of the pollution field - where the random forest disagrees most between trees - rather than simply the density of the monitoring network.

Together the panels support the modelling choice. For PM2.5 at the national scale the signal is carried mainly by environmental covariates (urban–rural context, population), so the reconstruction is sharpest, but also least certain, exactly where concentrations are highest and most heterogeneous. This is the regime in which an RF-assisted kriging hybrid is most useful, and it highlights where additional monitoring would most reduce uncertainty.

7.5 Accuracy comparison (cross-validation)

# Single 80/20 split: score every method on the held-out stations
idw_t <- gstat::idw(conc ~ 1, train_sf, test_sf, idp = 2, debug.level = 0)$var1.pred
ok_t  <- autoKrige(conc ~ 1, train_sf, test_sf, verbose = FALSE)$krige_output$var1.pred

## [using ordinary kriging]

uk_t  <- autoKrige(conc ~ pop_dens + dist_centre, train_sf, test_sf,
                   verbose = FALSE)$krige_output$var1.pred

## [using universal kriging]

test_sf$rf_fit <- predict(rf, st_drop_geometry(test_sf))$predictions
rfk_t <- test_sf$rf_fit +
  autoKrige(resid ~ 1, train_sf, test_sf, verbose = FALSE)$krige_output$var1.pred

## [using ordinary kriging]

scoreboard <- data.frame(
  Method = c("IDW", "Ordinary kriging", "Universal kriging", "RF + kriging"),
  RMSE = c(rmse(test_sf$conc, idw_t), rmse(test_sf$conc, ok_t),
           rmse(test_sf$conc, uk_t),  rmse(test_sf$conc, rfk_t)),
  MAE  = c(mae(test_sf$conc, idw_t),  mae(test_sf$conc, ok_t),
           mae(test_sf$conc, uk_t),   mae(test_sf$conc, rfk_t))
)
scoreboard

# Ten-fold cross-validation
set.seed(2026)
k <- 10
folds <- sample(rep(1:k, length.out = nrow(stations_sf)))

cv_one <- function(train, test) {
  idw <- gstat::idw(conc ~ 1, train, test, idp = 2, debug.level = 0)$var1.pred
  ok  <- autoKrige(conc ~ 1, train, test, verbose = FALSE)$krige_output$var1.pred
  uk  <- autoKrige(conc ~ pop_dens + dist_centre, train, test,
                   verbose = FALSE)$krige_output$var1.pred
  rf  <- ranger(conc ~ pop_dens + dist_centre + x + y,
                data = st_drop_geometry(train), num.trees = 500)
  rfk <- predict(rf, st_drop_geometry(test))$predictions +
    autoKrige(resid ~ 1,
              train %>% mutate(resid = conc - predict(rf, st_drop_geometry(train))$predictions),
              test, verbose = FALSE)$krige_output$var1.pred
  tibble(IDW = idw, OK = ok, UK = uk, RFK = rfk, obs = test$conc)
}

cv <- purrr::map_dfr(1:k, function(i)
  cv_one(stations_sf[folds != i, ], stations_sf[folds == i, ]))

## [using ordinary kriging]
## [using universal kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using universal kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using universal kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using universal kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using universal kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using universal kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using universal kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using universal kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using universal kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using universal kriging]
## [using ordinary kriging]

cv %>% summarise(across(IDW:RFK, ~ Metrics::rmse(obs, .x)))  # CV-RMSE

Ten-fold cross-validation identified universal kriging as the best-performing method (CV-RMSE = 2.37 µg/m³), ahead of inverse distance weighting, ordinary kriging and RF-assisted kriging, which were similar (CV-RMSE ≈ 2.53–2.55 µg/m³). The advantage of universal kriging - the only method that models the covariate trend explicitly while retaining a spatial component - indicates that the reconstruction of PM2.5 across Poland is driven mainly by an environmental-covariate trend (population density and distance to the nearest urban centre) rather than by residual spatial autocorrelation. This is consistent with the near-flat residual variogram and the covariate-driven uncertainty surface. The non-linear random forest did not improve on the simple linear trend - with only 95 monitoring stations the relationship is smooth enough that a linear trend is sufficient and more stable, whereas the random forest is prone to overfitting.

Inverse distance weighting served as the deterministic baseline. In the single 80/20 split it gave the lowest error, but in ten-fold cross-validation it was average (CV-RMSE = 2.54 µg/m³), comparable to ordinary kriging and RF-kriging and behind universal kriging. This pattern was expected. IDW uses only distance to the observations and ignores the covariate trend and the spatial correlation structure, so it performs adequately where the field is smooth but cannot exploit the population/urban-structure signal that drives PM2.5.

8 Population exposure assessment

The reconstructed surface lives on the same cells as the population grid, so exposure follows directly. There are three questions to answer:

1. What is the average person breathing?
2. How many people are above the limits?
3. How does this differ between cities and the countryside?

exposure <- st_drop_geometry(pop_grid) %>%
  transmute(pop, pred_uk, degurba = grid_sf$degurba)

# 1. Population-weighted mean concentration (average exposure)
pop_weighted_mean <- with(exposure, sum(pred_uk * pop) / sum(pop))

# 2. Share of residents above WHO and EU thresholds
share_above <- function(th) with(exposure, sum(pop[pred_uk > th]) / sum(pop))
share_who   <- share_above(params$who_threshold)
share_eu    <- share_above(params$eu_threshold)
share_eu_2030 <- share_above(params$eu_threshold_2030)

tibble::tibble(
  metric = c("Population-weighted mean [ug/m3]",
             paste0("% population > WHO (", params$who_threshold, ")"),
             paste0("% population > EU current (", params$eu_threshold, ")"),
             paste0("% population > EU 2030 (",  params$eu_threshold_2030, ")")),
  value  = c(round(pop_weighted_mean, 1),
             percent(share_who,      accuracy = 0.1),
             percent(share_eu,       accuracy = 0.1),
             percent(share_eu_2030,  accuracy = 0.1))
)

# 3. Urban–rural inequality - exposure by DEGURBA class
exposure_by_class <- exposure %>%
  group_by(degurba) %>%
  summarise(
    population    = sum(pop),
    weighted_mean = sum(pred_uk * pop) / sum(pop),
    pct_above_who      = sum(pop[pred_uk > params$who_threshold])      / sum(pop),
    pct_above_eu       = sum(pop[pred_uk > params$eu_threshold])       / sum(pop),
    pct_above_eu_2030  = sum(pop[pred_uk > params$eu_threshold_2030])  / sum(pop),
    .groups = "drop"
  ) %>%
  mutate(across(starts_with("pct"), ~ percent(.x, accuracy = 0.1)),
         weighted_mean = round(weighted_mean, 1))

exposure_by_class

pop_grid$band <- cut(
  pop_grid$pred_uk,
  breaks = c(-Inf, params$who_threshold, params$eu_threshold_2030,
             15, 20, Inf),
  labels = c("≤5 (WHO)", "5–10", "10–15", "15–20", ">20"),
  right = TRUE
)

ggplot(pop_grid) +
  geom_sf(aes(fill = band), color = NA) +
  scale_fill_viridis_d(option = "C", direction = -1, name = "PM2.5 [ug/m3]") +
  theme_minimal() +
  labs(title = "PM2.5 relative to WHO (5) and EU-2030 (10) thresholds",
       subtitle = "Whole population above WHO and EU-2030")

pop_grid$excess_2030 <- pmax(pop_grid$pred_uk - params$eu_threshold_2030, 0)

ggplot(pop_grid) +
  geom_sf(aes(fill = excess_2030), color = NA) +
  scale_fill_viridis_c(option = "B", name = "ug/m3 over 10") +
  theme_minimal() +
  labs(title = "PM2.5 exceedance above the EU-2030 limit (10 ug/m3)")

The exposure assessment overlays the reconstructed PM2.5 surface on the population grid, broken down by Degree of Urbanisation (DEGURBA). Average exposure is high and remarkably uniform across settlement types: the population-weighted annual mean is 14.7 µg/m³ in rural areas, 15.2 in urban clusters and 15.4 in urban centres. In every class the entire population lives above both the WHO 2021 guideline (5 µg/m³) and the future EU 2030 annual limit (10 µg/m³), while none is above the current EU limit (25 µg/m³). Because these binary indicators saturate at 100%, the spatial distribution is more informative: concentration bands and the exceedance margin above the EU 2030 limit both show that levels are highest in central–southern Poland (Silesia, Kraków and Łódź regions and southern Mazovia), where annual means reach 17–18 µg/m³ and exceed the 2030 limit by 7–8 µg/m³, and lowest across the north (10–14 µg/m³, exceeding the 2030 limit by only ~2 µg/m³).

9 Conclusions

This project reconstructed the 2024 PM2.5 concentration surface for Poland from 95 monitoring stations, using spatial machine learning, and turned that surface (reconstructed with universal kriging, the cross-validation-best method) into a population exposure assessment. Ten-fold cross-validation ranked universal kriging first (CV-RMSE = 2.37 µg/m³), ahead of IDW, ordinary kriging and RF-kriging (2.53–2.55), showing that the field is driven by an environmental-covariate trend (population, distance to urban centres) rather than by residual spatial autocorrelation. The population-weighted mean is 15 µg/m³ - no residents exceed the current EU limit (25 µg/m³), yet effectively 100% exceed both the WHO guideline (5 µg/m³) and the future EU 2030 limit (10 µg/m³), with the highest concentrations and exceedance margins concentrated in the central–southern Poland. Exposure is high but uniform between cities and countryside (14.7–15.3 µg/m³ across DEGURBA classes), reflecting widespread residential heating. The study delivers a transparent, reproducible reconstruction with an uncertainty layer and a population-weighted, urban–rural exposure breakdown.

10 References

Hengl, T., Nussbaum, M., Wright, M. N., Heuvelink, G. B. M., & Gräler, B. (2018). Random forest as a generic framework for predictive modeling of spatial and spatio-temporal variables. PeerJ, 6, e5518.

Hoek, G., Beelen, R., de Hoogh, K., Vienneau, D., Gulliver, J., Fischer, P., & Briggs, D. (2008). A review of land-use regression models to assess spatial variation of outdoor air pollution. Atmospheric Environment, 42(33), 7561–7578.

World Health Organization (2021). WHO Global Air Quality Guidelines. Geneva: WHO.

Eurostat (2021). Applying the Degree of Urbanisation: A Methodological Manual to Define Cities, Towns and Rural Areas for International Comparisons. Publications Office of the European Union.

CIEP air-quality measurements (https://powietrze.gios.gov.pl/pjp/archives).

GUS population grid, Census 2021 (https://geo.stat.gov.pl/).

GHSL/DEGURBA (https://human-settlement.emergency.copernicus.eu/).