arXiv CSopen access

Air Quality Downscaling with Station-Guided Pseudo-Supervision

Air Quality Downscaling with Station-Guided Pseudo-Supervision Guorun Wang1 , Simone Foti1 , Andreas D. Demou2 , Leonidas Kotoulas1 , Theodoros Christoudias2 , Alexandros Koliousis3 , Mihalis Nicolaou4 , and Stefanos Zafeiriou1 Imperial College London, UK The Cyprus Institute, Cyprus 3 Northeastern University London, UK 4 University of Cyprus, Cyprus

arXiv:2607.05292v1 [cs.LG] 6 Jul 2026

1

2

Abstract. Super-resolving coarse atmospheric fields to local PM2.5 variations is uniquely challenged by a mismatch in spatial support: while pixels represent regional averages, ground-truth observations are discrete, unaligned samples of a continuous spatial signal. To bridge this gap, we present a station-guided framework for high-resolution PM2.5 downscaling over Europe. Taking coarse CAMS atmospheric composition fields alongside heterogeneous side information (i.e., human activity, land cover, elevation, satellite aerosol observations, and wind fields) our framework jointly super-resolves (×40, ≈ 1 km) and bias-corrects CAMS rasters, without relying on temporal sequence modelling. To address the challenge of densely supervising our multi-scale transformer network with sparse in-situ data, we introduce a time-agnostic propagation strategy that utilises spatial Gaussian blending of interpolated OpenAQ observations. Extensive qualitative and station-level evaluations across Europe demonstrate that our model recovers fine-grained spatial structures and effectively mitigates localised CAMS biases. Keywords: PM2.5 downscaling · Air-quality super-resolution · Stationguided learning · Gaussian kernel interpolation · Earth-system AI

1

Introduction

Climate change is reshaping atmospheric conditions in ways that directly influence air quality through complex interactions between meteorology, atmospheric chemistry, and human activities. As climatic conditions evolve, changes in temperature, circulation patterns, and the frequency of extreme events can affect the formation, transport, and dispersion of air pollutants, creating new challenges for air-quality management and forecasting. Air quality remains a major environmental and public health concern, affecting billions of people worldwide and contributing significantly to respiratory and cardiovascular disease burdens [17]. Accurate and timely predictions are therefore essential for supporting public health interventions, informing environmental policy, and enhancing our understanding of Earth system processes. In this context, machine-learning-based

2

G. Wang et al.

forecasting systems offer a promising pathway for improving the prediction of pollution episodes and informing public health decision-making. The rapid development of machine learning (ML) systems for scientific prediction and decision-making has ushered in a new era for environmental forecasting [37, 45, 60]. In meteorology, this progress has been particularly evident in the emergence of advanced AI-based weather forecasting systems, which have achieved remarkable performance in predicting dynamical atmospheric variables [4, 6, 7, 13, 20, 24, 26, 28, 29, 34, 35, 52, 53]. However, extending these successes from weather prediction to atmospheric-composition forecasting remains challenging. Unlike conventional weather forecasting, atmospheric-composition prediction must account not only for meteorological transport but also for emissions, chemical transformations, deposition, and human activities, making it substantially more complex and computationally demanding [5]. This challenge is especially important for air pollution, which directly affects human health and ecosystems [11, 38]. To tackle this, Aurora [5], a foundation model for the Earth system, shows a pathway for moving beyond weather forecasting towards more general Earthsystem prediction, including atmospheric composition. ML has demonstrated strong capabilities in learning the complex nonlinear relationships between meteorology, emissions, and pollutant concentrations, an approach that has shown promising improvements in both predictive accuracy and computational efficiency, marking a paradigm shift from purely physics-based models toward MLbased frameworks [21, 23, 27, 42, 43, 50, 54, 59]. Despite these promising developments, current methods such as Aurora typically operate at spatial resolutions of 0.25^\circ \text {â•fi}0.4^\circ (\approx 25\text {â•fi}40\,\text {km}), which is still too coarse to capture pollution hotpspots and regional variations. Hence, an important direction is spatial downscaling (super-resolution), which aims to produce point-level predictions from coarse global reanalyses. On the data side, air quality monitoring and forecasting rely on two main data sources with complementary characteristics. On the one hand, forecast and analysis data such as those from the Copernicus Atmosphere Monitoring Service (CAMS) provide spatially complete and consistent estimates of atmospheric composition, particularly where observational data coverage is low or for atmospheric pollutants for which no direct observations are available [9, 10, 25]. Reanalysis data combine physical models with past observations through a process called data assimilation. On the other hand, ground-based monitoring stations such as EBAS [55] and OpenAQ [1] provide highly accurate point-level measurements of pollutants like fine particulate matter (PM2.5 ), regarded as the ground truth for model evaluation. However, each data source has its own limitations. For reanalysis and forecast products, the spatial resolution is relatively coarse (approximately 0.4◦ , or about 40 km). Moreover, studies have reported that CAMS exhibits systematic biases relative to ground-based measurements, often showing overestimation of particulate matter concentrations [14, 18]. In contrast, the station measurements are spatially sparse and unevenly distributed, with higher densities in particular areas and limited coverage in rural or remote

Air Quality Downscaling with Station-Guided Pseudo-Supervision

3

areas. This spatial heterogeneity makes it difficult to construct high-resolution air-quality fields solely from station data. The key challenge, therefore, lies in bridging these two data sources and combining the spatial completeness of CAMS with the precision of ground stations to produce high-resolution, bias-corrected air-quality estimates. Auxiliary variables that encode physical and environmental characteristics can provide additional information in this process. Incorporating high-resolution side information enables the model to learn physically-grounded relationships, and to propagate point observations from stations into the coarse CAMS grids, thereby improving downscaling performance. In this study, we focus on PM2.5 downscaling over Europe using OpenAQ station measurements as the ground-truth reference. Starting from the 0.4◦ CAMS forecast fields as the prior, we integrate multiple auxiliary information sources: human activity, land use, elevation, aerosol satellite observations, and wind to guide the learning process. The objective is to construct a model that effectively propagates the accuracy of the OpenAQ dataset to the broad spatial coverage of CAMS, achieving physically consistent and spatially detailed PM2.5 estimates, thereby enabling a resolution enhancement from 0.4◦ to 0.01◦ . Our main contributions are summarized as follows: – Task formulation. We recast fine-grained PM2.5 estimation as a timeagnostic, joint downscaling and bias-correction problem rather than temporal forecasting. Given coarse CAMS fields and time-matched side information, our model produces ×40 (≈1 km) bias-corrected PM2.5 maps at any single timestamp, without temporal sequence modeling, allowing it to generalize to unseen times and locations. – Gaussian-kernel pseudo-label propagation. To densely supervise the network from spatially sparse, grid-unaligned station data, we introduce a propagation scheme that spreads OpenAQ observations into a dense target and blends it with the CAMS prior where station support is weak. This turns a few thousand point measurements into pixel-wise supervision representing the ground-truth accuracy directly into the coarse field. – Continental-scale framework and evidence. We propose STARQ (STation AiR Quality), our continental-scale framework over all Europe, fusing eight heterogeneous geospatial and atmospheric sources within a multi-scale SegFormer backbone. Under a strict spatiotemporal split (unseen timestamps and unseen stations), it consistently outperforms both CAMS and S-MESH∗ , our implementation of S-MESH [49], the XGBoost-based stateof-the art method for station-guided downscaling.

2

Related Work

Downscaling for PM2.5 . High spatiotemporal-resolution PM2.5 concentration data are essential for air-pollution assessment and research. Over the past decade, PM2.5 downscaling methods have evolved from traditional statistical approaches

4

G. Wang et al.

C

D

12

10

10 6

C

D 10

8

15

Pred. closeup @ Stations

Downscaled Pred.

6

A

B A C

9

B 11 7

C

11

D

13

PM2.5 (μg m-3)

B A

CAMS closeup + GT Stations

CAMS Baseline

12

B

A

4

2

D 0

Fig. 1: PM2.5 over Europe from CAMS forecast (top) and downscaled by STARQ (bottom). Close-ups (right) are reported alongside ground truth station values and our corresponding predictions.

to modern ML and deep learning frameworks. Early PM2.5 downscaling methods typically use coarse-scale model outputs as inputs and transform them into highresolution estimates through statistical relationships. For example, [56] downscale CAMS PM2.5 using CAMS as a prior together with urban monitoring data from Budapest. However, the inherent linearity of such approaches may introduce systematic biases. In contrast, geostatistical methods such as kriging [36] explicitly model spatial autocorrelation, although they require dense monitoring networks to perform reliably. Some methods also address missing auxiliary data. MODIS aerosol optical depth (AOD), a commonly used source of side information, often contains missing values. To address this issue, [32] employ a Bayesian-based statistical downscaler to model the spatio-temporal linear relationships between AOD and PM2.5 and to fill missing AOD values. In parallel, [64] downscale the 0.1° ACAG PM2.5 product to approximately 300 m using a cascade random forest model driven by elevation, land-cover data, and AOD. [41] impute 1 km MAIAC AOD using a quantile regression forest and then use a gradient boosting machine algorithm to estimate ground-level PM2.5 . Overall, these methods improve resolution and accuracy to some extent, but they typically rely on extensive feature engineering and strong modelling assumptions. Recent studies further explore PM2.5 downscaling and forecasting across different regions and temporal scales. [50] analyze hourly air-pollutant data across 11 global cities during the COVID-19 lockdowns, using station-level spatial resolution and hourly temporal resolution. [54] develop a hybrid GNN–LSTM model for 72-hour PM2.5 forecasting over the Beijing–Tianjin–Hebei region, using hourly data from 2016–2020 at station-level spatial resolution. In the same

Air Quality Downscaling with Station-Guided Pseudo-Supervision

5

region, [42] use data from 2020–2022 at 3-hour temporal and 9 km spatial resolution for PM2.5 forecasting. [27] focus on PM2.5 forecasting in Seoul, South Korea, operating at 6-hour temporal and 27 km spatial resolution. [23] conduct hourly PM2.5 forecasting across 18 stations in the Taipei metropolitan area, using data from 2014–2020 at 1-hour temporal resolution. [43] estimate global daily PM2.5 concentrations at 0.5° × 0.625° resolution from 1980–2023, achieving high accuracy and revealing persistent pollution extremes in South and East Asia. However, these studies either focus on relatively small regions or produce results at comparatively coarse spatial resolution. Many also depend on temporal information or forecasting-specific settings. These limitations motivate our time-independent downscaling task over a large European domain. Recent European-scale downscaling methods further incorporate station observations. [49] employs station-supervised XGBoost [8] to downscale CAMS regional PM2.5 forecasts into 1 km daily PM2.5 maps, while [19] trains a framework to generate 500 m annual air-quality maps for multiple pollutants across Europe. However, neither method leverages hourly OpenAQ station observations for timestamp-level supervision, which would enable PM2.5 downscaling at arbitrary hours. In contrast, our work addresses this gap by incorporating hourly station observations to support timestamp-level PM2.5 downscaling. Convolutional and Transformer Architectures. Neural networks can automatically extract complex spatial feature relationships and demonstrate strong nonlinear fitting capabilities in downscaling tasks. For PM2.5 , [63] successfully applied spatial convolution to the downscaling of meteorological variables, demonstrating that spatial convolution is suitable for refining continuous fields such as air pollution. Convolutional architectures are powerful tools for image processing and have been widely applied to super-resolution and climate data downscaling [2, 3, 15, 58, 61]. Encoder-decoder networks, such as U-Net and Transformer-based models, can effectively fuse multi-scale features [30, 44, 57, 62]. We recognize the strong performance of convolutional and transformer architectures in downscaling tasks, and our approach builds on their respective advantages. However, existing studies often fail to achieve sufficiently fine spatial resolution, such as around 1 km, are limited in their ability to cover large regions, such as intercontinental or global domains, or rely on temporal dependencies that require time-series contextual information. In contrast, our method achieves a very high spatial resolution of 0.01° while maintaining broad spatial coverage across the whole of Europe, and it does not depend on temporal inputs. This means that, given the required spatial information at any point in time, our model can perform downscaling accordingly. Gaussian Kernel Interpolation for PM2.5 . The Gaussian downscaling method, originally introduced by [48, 51] to characterize the spatial dispersion of pollutants, was adopted by [47] to refine coarse 50 km WRF/Chem PM2.5 outputs to a 5 km resolution. Thus, the importance of Gaussian-related formulas for climate studies is evident. Moreover, the reliability of interpolation methods has already been demonstrated in reconstructing temperature fields [31]. In this study, we

6

G. Wang et al.

Dataset

Spatial Resolution

Temporal Resolution

CAMS forecast GHSL Built-up Surface (GHS-BUILT-S) GHSL Built-up Volume (GHS-BUILT-V) GHSL Population Grid (GHS-POP) Grouping Land Use CLC EU-DEM MODIS MCD19A2 ERA5 wind (single levels) OpenAQ

0.4◦ 3 arcsec 3 arcsec 3 arcsec 100 m 1 arcsec 1 km 0.25◦ point measurements

hourly 2015, 2020, 2025 2015, 2020, 2025 2015, 2020, 2025 2012, 2018 static daily hourly provider-dependent

Table 1: Spatial and temporal resolutions of the datasets.

used the Gaussian kernel to perform interpolation on the propagation of station measurements.

3

Data

Data Sources. We use eight geospatial and atmospheric raster datasets as inputs. Specifically, CAMS global atmospheric composition forecasts [10] serve as the low-resolution PM2.5 rasters to downscale and bias-correct. Global Human Settlement Layer (GHSL) provides anthropogenic activity information such as the built-up surface (GHS-BUILT-S) [39], built-up volume (GHS-BUILTV) [40], and population (GHS-POP) [46]. In addition, CORINE Land Cover (CLC) [12]characterises the land-use conditions, EU-DEM [16] captures the topography, MODIS MCD19A2 Aerosol Optical Depth (AOD) [33] measures the satellite-observed aerosol loading, and ERA5 10-m [22] provides the wind components for inferring atmospheric transport and dispersion. The spatial and temporal resolutions of these datasets are summarised in Table 1. All input variables are spatially aligned and resampled onto a common 0.01◦ European grid, while time-varying variables are matched to the corresponding OpenAQ observation timestamps. Detailed descriptions of the datasets and their preprocessing procedures are provided in Sections A and B of the Supplementary Material. While all the aforementioned input datasets are regular raster data, the ground-truth PM2.5 measurements used for supervision are sparse point observations obtained from OpenAQ, an open-source global air-quality data aggregation and sharing platform. The station measurements in OpenAQ come from various providers, including governments, research institutions, and open community networks. For the purpose of this study we select stations within Europe, regardless of their provider. In addition, we curate a significantly smaller and computationally tractable dataset from 2018, 2020, and 2023, covering the Italian peninsula and part of the Balkans, (spanning 36◦ N–47◦ N and 6◦ E–19◦ E). We hereafter refer to this subset as the Italian dataset. Since our method produces a super-resolved PM2.5 rasters with a 0.01◦ spatial resolution, we project station measurements to this resolution by averaging temporally concurrent observations from all stations within the same pixel. Fol-

Air Quality Downscaling with Station-Guided Pseudo-Supervision

7

(a) Spatial distribution of OpenAQ monitoring stations (b) Distribution of OpenAQ PM2.5 obover the European domain. servations.

Fig. 2: Spatial coverage and value distribution of the OpenAQ PM2.5 observations used in this study.

lowing this projection, the dataset contains 3,952 unique grid-aligned stations over Europe. As can be observed in Figs. 2a and 2b, not only is the stations’ spatial distribution highly non uniform, but the PM2.5 observations are also strongly right-skewed, with a median of 8.00 µg/m3 and a mean of 13.65 µg/m3 . Most observations lie in low-to-moderate concentration ranges, while a small number of high-concentration and extreme values form a long tail. These properties make the learning problem challenging: the model must generalize across uneven station coverage while remaining robust to rare high-pollution events.

4

Methodology

4.1

Model Architecture

Since the input is a multi-channel geospatial raster, we first apply a shallow convolutional network to project the input channels into a fixed feature dimension. Then, the projected feature map is passed to a four-stage hierarchical Transformer based on the SegFormer encoder. The multi-scale features computed by the four encoder stages are projected to a common decoder dimension, upsampled to the same spatial resolution, concatenated, and fused. The fused representation is then upsampled to the original patch resolution and passed through a convolutional regression head to predict a single-channel PM2.5 concentration map. 4.2

Gaussian-Kernel Pseudo-Label Propagation

Sparse station observations provide accurate PM2.5 measurements only at a limited number of monitoring sites, whereas our goal is to produce dense pixel-wise predictions. To bridge this gap, we introduce a station-guided pseudo-label propagation strategy. The method uses a Gaussian kernel to propagate point-level

8

G. Wang et al.

station observations to nearby pixels and blends the resulting interpolated field with a coarse baseline estimate in regions with station support. Let Ω denote the spatial domain and let S = {si }N i=1 be the set of monitoring stations with valid PM2.5 observations yi . For each pixel location p ∈ Ω, we denote the model prediction by ŷ(p) and the baseline estimate, e.g., CAMS, by b(p). Gaussian station influence. Each station contributes to nearby pixels through a Gaussian influence function: w_i(\mathbf {p}) = \exp \left ( -\frac {\|\mathbf {p}-\mathbf {s}_i\|_2^2}{2\sigma ^2} \right ), \label {eq:gaussian_weight}

(1)

where σ controls the spatial decay of station influence. The station-interpolated field is then computed as: \tilde {y}(\mathbf {p}) = \frac { \sum _{i=1}^{N} w_i(\mathbf {p}) y_i }{ \sum _{i=1}^{N} w_i(\mathbf {p}) + \epsilon }, \label {eq:station_interpolation}

(2)

where ϵ is a small constant for numerical stability. Pseudo-label generation. We further define a cumulative station-confidence map: W(\mathbf {p}) = \max \left (0, \min \left (1, \sum _{i=1}^{N} w_i(\mathbf {p}) \right ) \right ), \label {eq:station_confidence}

(3)

which measures how strongly a pixel is supported by nearby stations. The final pseudo-label is obtained by blending the station-interpolated field with the baseline estimate: y_{\mathrm {pseudo}}(\mathbf {p}) = W(\mathbf {p})\tilde {y}(\mathbf {p}) + \left (1-W(\mathbf {p})\right )b(\mathbf {p}). \label {eq:pseudo_label}

(4)

At station locations, the pseudo-label is replaced by the observed ground truth, and these locations do not contribute to the pseudo-label supervision loss: y_{\mathrm {pseudo}}(\mathbf {s}_i) = y_i. \label {eq:station_override}

(5)

Training objective. The training loss combines station-level supervision from observed measurements with dense pseudo-label supervision over the remaining spatial domain: \mathcal {L}_{\mathrm {station}} &= \frac {1}{N} \sum _{i=1}^{N} \left ( \hat {y}(\mathbf {s}_i) - y_i \right )^2, \label {eq:station_loss} \\ \mathcal {L}_{\mathrm {pseudo}} &= \frac {1}{|\Omega \setminus \mathcal {S}|} \sum _{\mathbf {p}\in \Omega \setminus \mathcal {S}} \left ( \hat {y}(\mathbf {p}) - y_{\mathrm {pseudo}}(\mathbf {p}) \right )^2, \label {eq:pseudo_loss} \\ \mathcal {L}_{\mathrm {total}} &= \lambda _s \mathcal {L}_{\mathrm {station}} + \lambda _p \mathcal {L}_{\mathrm {pseudo}}, \label {eq:total_loss}

(8)

Air Quality Downscaling with Station-Guided Pseudo-Supervision

9

(a) Effect of kernel width σ.

(b) Pseudo-label and loss-weight construction for the best hyperparameters.

(c) One-dimensional cross-section. Fig. 3: Visualization of Gaussian pseudo-label propagation. Increasing σ expands station influence, while the pseudo-label blends station interpolation with the CAMS baseline.

where λs and λp control the relative contributions of the station-level supervision and the dense pseudo-label supervision, respectively. We visualize the behaviour of the Gaussian pseudo-label propagation in Fig. 3. The one-dimensional cross-section further illustrates how the pseudolabel transitions from station observations back to the CAMS baseline.

10

G. Wang et al.

5

Experiments

5.1

Implementation Details

Architecture and hyperparameters. The encoder adopts a custom fourstage SegFormer configuration with hidden dimensions of [256, 512, 1280, 2048] and transformer depths of [3, 8, 27, 3]. Because the input channels correspond to meteorological, chemical, and geographical variables rather than natural RGB images, the model is trained from scratch by minimising Eq. (8). The best performing model is selected based on the epoch with the lowest validation MAE with respect to station data. Hyperparameter tuning was performed on the Italian dataset using Ray with Optuna, and using the validation MAE as optimization objective. The best configuration, subsequently used also on the full European dataset, consisted of a background sampling ratio of 0.176, a learning rate of 1.31 × 10−5 , a weight decay of 1.95 × 10−5 , a Gaussian standard deviation of σ = 12.32 pixels, as well as station- and pseudo-label-supervision weights of λs = 0.173 and λp = 0.970, respectively. Data splits and patch extraction. To rigorously evaluate generalization, we perform a spatiotemporal split. Timestamps are split into train, validation, and test sets with a ratio of 70:20:10, while station locations are split with a ratio of 80:10:10. This ensures the model is evaluated on unseen side information and station locations. We extract 64 × 64 input patches (covering approximately 64 × 64 km) cropped around monitoring stations with random spatial offsets. We also sample background patches without station supervision to encourage the model to preserve the large-scale spatial priors learned from the CAMS baseline. While during training patches are not necessarily centred with respect to a station, during validation and testing, the target station is always placed at the centre of the patch, and metrics are computed only at the site location. Bucketed timestamp sampling. Fully random patch sampling across all timestamps creates a severe I/O bottleneck, while sampling from a single timestamp destroys the temporal diversity required for stable gradients. We balance data locality and batch diversity using a bucketed timestamp sampling strategy (K = 4). Timestamps are randomly shuffled and grouped into buckets of four. All patches associated with these timestamps are pooled, randomly shuffled, and sequentially yielded as mini-batches. This restricts data access to a small, cached subset of files while retaining temporal variation within each batch. Hardware. All experiments were run on 4 NVIDIA H200 NVL GPUs. For the main experiments, we used a per-GPU batch size of BATCH_SIZE=384. For hyperparameter tuning, we launched 20 concurrent Optuna trials on a single NVIDIA H200 NVL GPU. 5.2

EU Downscaling Results

The evaluation results for the EU dataset shown in Table 2 are measured in terms of mean absolute error (MAE), root mean squared error (RMSE), and

Air Quality Downscaling with Station-Guided Pseudo-Supervision

11

Table 2: Performance comparison on the EU test set between the CAMS baseline, S-MESH∗ , and our Gaussian-kernel pseudo-label propagation method. MAE ↓ RMSE ↓ R2 ↑ CAMS Forecast 6.562 13.723 0.029 S-MESH∗ 6.485 13.080 0.118 STARQ (ours) 5.873 12.126 0.242 Table 3: Performance comparison across broad land-use categories. MAE and RMSE are reported in µg/m3 . Lower MAE and RMSE are better, while higher R2 is better. The best result for each metric within each land-use category is highlighted in bold. S-MESH∗

STARQ (ours) Land use

N MAE ↓ RMSE ↓ R ↑ MAE ↓ RMSE ↓ 2

Urban 199,020 6.297 14.274 0.210 6.893 Industrial 70,342 5.250 9.275 0.320 5.929 Agriculture 62,444 5.431 8.797 0.254 5.891 Natural 28,506 5.565 9.008 0.259 6.348 Water 16,621 4.986 7.872 0.219 5.538 No Data/Unknown 1,105 15.222 27.055 0.366 19.711

CAMS Forecast R2 ↑ MAE ↓ RMSE ↓

15.213 0.103 6.945 10.288 0.163 6.005 9.391 0.150 5.740 10.476 −0.002 6.348 8.602 0.067 6.611 34.391 −0.024 24.535

R2 ↑

15.862 0.025 10.946 0.053 9.787 0.076 10.395 0.013 10.506 −0.392 40.260 −0.404

coefficient of determination (R2 ) between predicted values at station locations and corresponding ground-truth observations. Our model outperforms both the baseline CAMS Forecast and S-MESH∗ , our implementation of S-MESH trained on the same data used to train STARQ. We exhibit lower MAE and RMSE, as well as higher R2 . Following the methodology adopted in [2], we divided the CLC dataset into 15 classes for model training and evaluation. To more clearly demonstrate the advantages of our model, we further grouped these classes into five broader categories and evaluated the improvements achieved across heterogeneous land-use conditions. As shown in Table 3, STARQ consistently achieves the best performance across all broad land-use categories in terms of MAE, RMSE, and R2 . The improvement is observed not only in the dominant urban category, which contains more than half of the test observations, but also across industrial, agricultural, natural, and water-covered regions. The largest improvement over CAMS Forecast is observed in the water category, where STARQ reduces MAE from 6.611 to 4.986 and improves R2 from −0.392 to 0.219. S-MESH∗ generally improves upon CAMS Forecast, but its gains are less consistent and substantially smaller than those achieved by STARQ. These results indicate that explicitly modelling spatial context provides more robust improvements across heterogeneous land-use conditions than point-wise nonlinear regression. For full-domain visualization, we apply overlapping sliding-window inference with window size Ws and overlap O. Predictions from overlapping windows are merged using normalized weighted averaging, where pixels near patch boundaries receive lower weights than those near the patch centre. This boundary-aware blending suppresses patch artifacts and produces a spatially smooth prediction

12

G. Wang et al.

Fig. 4: PM2.5 distribution over Europe.

Fig. 5: Regional PM2.5 predictions over Europe.

over the entire European domain. Further details of the padding, weighting, and aggregation procedures are provided in Section C of the Supplementary Material. Figures 4 and 5 show the qualitative results of our high-resolution PM2.5 estimation. Compared with the original CAMS Forecast field at a spatial resolution of 0.4◦ , our model produces predictions on a 0.01◦ grid, revealing much finer local spatial structures. In addition to increasing the spatial granularity, the model also uses station observations to locally correct the CAMS baseline, leading to more realistic concentration patterns around observed regions. Figure 1 further demonstrates that the proposed method can generate visually consistent high-resolution PM2.5 maps over large-scale geographical regions. To assess whether our model can track temporal variations at fixed sites, we present time-series case studies at two representative validation stations, shown in Figure 6. This comparison includes the CAMS baseline, the proposed model prediction, and the station ground-truth observations at the same locations. Figure 6(a) presents the results for a station with pronounced pollution episodes. In this case, CAMS systematically underestimates high-PM2.5 values, especially during winter peaks. In contrast, the proposed model better follows the temporal variability of the ground-truth observations and reduces the MAE from 17.43 to 8.89. Figure 6(b) focuses on a station with relatively low concentration levels. Here, CAMS tends to overestimate PM2.5 over the time period. The proposed model corrects this positive bias and produces predictions closer to the observations, reducing the MAE from 7.05 to 4.89.

Air Quality Downscaling with Station-Guided Pseudo-Supervision

13

These results suggest that CAMS errors are location-dependent: the coarseresolution baseline can underestimate polluted episodes at some locations while overestimating concentrations at others. By learning spatial and temporal side information from other monitoring stations and side information, the proposed model mitigates these site-specific biases and improves prediction accuracy at unseen timestamps and locations.

(a) High-concentration case.

(b) Low-concentration case.

Fig. 6: Temporal comparison between ground-truth PM2.5 observations, CAMS, and STARQ at two representative validation stations.

5.3

Channel importance

We assess local channel importance on a representative validation region using three complementary methods: (i) leave-one-channel-out ablation, monitoring the increase in MAE after removing a channel, (ii) gradient-based sensitivity, measuring model output sensitivity to input channel perturbations, and (iii) permutation-based importance, evaluating how much the model relies on spatially and temporally aligned information in each channel. Each method is described in Section D of the Supplementary Material. Figure 7 summarizes the three local channel-importance analyses. The ablation and permutation results both measure changes in MAE and therefore directly reflect whether a channel improves predictive performance in this local region. Both methods identify CAMS PM2.5 and GHSL built-up features as the dominant inputs. This is consistent with the role of CAMS as a coarse atmospheric prior and GHSL built-up variables as proxies for urban morphology and emission-related spatial structure. Other auxiliary variables, such as CLC land cover, ERA5 wind components, and DEM elevation, provide smaller but generally positive contributions. In contrast, MODIS AOD and GHSL population show slightly negative importance in this local test, suggesting that these channels may be noisy, redundant, or less aligned with the local station-level PM2.5 correction. The gradient-based sensitivity assigns non-negligible sensitivity to most channels, indicating that the model responds to multiple inputs in this region. However, this metric should be interpreted differently from the MAE-based scores: a high gradient sensitivity reflects local responsiveness to a channel, but does not necessarily mean that the channel improves predictive accuracy. A channel can

14

G. Wang et al.

(a) Ablation.

(b) Gradient sensitivity.

(c) Permutation.

Fig. 7: Local channel-importance analysis on a representative validation region.

strongly influence the model’s output while still contributing noise or redundant information, which would not be captured by gradient-based measures alone. Therefore, the three methods are best viewed as complementary: ablation and permutation importance reveal which channels are beneficial for performance, while gradient sensitivity highlights which channels the model relies on most, regardless of their effect on error.

6

Conclusion

We introduce STARQ, a station-guided downscaling framework for bias-corrected PM2.5 estimation over Europe. The proposed method combines the spatial completeness of CAMS with the point-level accuracy of OpenAQ observations, while incorporating auxiliary variables that describe human activity, land cover, topography, aerosol loading, and meteorological transport. By using Gaussian-kernel interpolation and pseudo-label propagation, the framework converts sparse station measurements into dense supervision and enables time-agnostic spatial correction at arbitrary timestamps. Experiments on the European test set demonstrate that the proposed model consistently improves upon the CAMS baseline in station-level evaluation, achieving lower MAE and RMSE and a substantially higher R2 . Qualitative results further show that the model generates spatially detailed PM2.5 fields at 0.01◦ resolution and produces smoother, more realistic local structures than the original coarse CAMS field. The time-series case studies also suggest that the model can mitigate location-dependent CAMS biases, including both underestimation during high-pollution episodes and overestimation in low-concentration regions. Additional ablation and channel-importance analyses indicate that CAMS PM2.5 and built-up features provide particularly important information, while other environmental and meteorological variables contribute complementary spatial context. Several limitations remain. The quality of supervision depends on the coverage and reliability of OpenAQ station measurements, which are spatially uneven and may contain outliers. In addition, the current pseudo-label design relies on fixed kernel-based propagation, which may not fully capture complex pollutant transport under varying meteorological conditions. Future work will explore adaptive and uncertainty-aware kernels, stronger quality-control strategies

Air Quality Downscaling with Station-Guided Pseudo-Supervision

15

for station observations, and broader evaluation across additional pollutants, regions, and extreme pollution events.

Acknowledgements GW, SF, and SZ were supported by the EPSRC Turing AI Fellowship (Grant Ref: EP/Z534699/1): Generative Machine Learning Models for Data of Arbitrary Underlying Geometry (MAGAL). AD and MN were supported through the TensorICE project (EXCELLENCE/0524/0407), which is implemented under the social cohesion programme "THALIA 2021-2027", co-funded by the European Union through the Research and Innovation Foundation of Cyprus.

16

G. Wang et al.

References 1. OpenAQ : Retrieved from https://api.openaq.org . https://api.openaq.org (2023) 2. Ashiotis, G., Tsigkanos, E., Christoudias, T., Nicolaou, M.A.: Ai for air quality: Leveraging data fusion for deep downscaling of atmospheric pollutants. In: MACLEAN@ PKDD/ECML (2022) 3. Baño-Medina, J., Manzanas, R., Gutiérrez, J.M.: Configuration and intercomparison of deep learning neural models for statistical downscaling. Geoscientific Model Development 13(4), 2109–2124 (2020) 4. Bi, K., Xie, L., Zhang, H., Chen, X., Gu, X., Tian, Q.: Accurate medium-range global weather forecasting with 3d neural networks. Nature 619(7970), 533–538 (2023) 5. Bodnar, C., Bruinsma, W.P., Lucic, A., Stanley, M., Allen, A., Brandstetter, J., Garvan, P., Riechert, M., Weyn, J.A., Dong, H., et al.: A foundation model for the earth system. Nature pp. 1–8 (2025) 6. Chen, K., Han, T., Ling, F., Gong, J., Bai, L., Wang, X., Luo, J.J., Fei, B., Zhang, W., Chen, X., et al.: The operational medium-range deterministic weather forecasting can be extended beyond a 10-day lead time. Communications Earth & Environment 6(1), 518 (2025) 7. Chen, L., Zhong, X., Zhang, F., Cheng, Y., Xu, Y., Qi, Y., Li, H.: Fuxi: a cascade machine learning forecasting system for 15-day global weather forecast. npj climate and atmospheric science 6(1), 190 (2023) 8. Chen, T., Guestrin, C.: Xgboost: A scalable tree boosting system. In: Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining. pp. 785–794 (2016) 9. Copernicus Atmosphere Monitoring Service (CAMS): Cams european air quality reanalyses. Copernicus Atmosphere Monitoring Service (CAMS) Atmosphere Data Store (2021). https://doi.org/10.24381/7cc0465a, https://ads.atmosphere. copernicus.eu/cdsapp#!/dataset/cams-europe-air-quality-reanalyses, accessed on 25-Oct-2025 10. Copernicus Atmosphere Monitoring Service (CAMS): Cams global atmospheric composition forecasts. Copernicus Atmosphere Monitoring Service (CAMS) Atmosphere Data Store (2021). https://doi.org/10.24381/04a0b097, https://ads. atmosphere . copernicus . eu / cdsapp # ! / dataset / cams - global - atmospheric composition-forecasts, accessed on 25-Oct-2025 11. Copernicus Atmosphere Monitoring Service (CAMS): Cams assessment report on european air quality in 2024. https://atmosphere.copernicus.eu/node/1330 (2025), doi:10.24380/pdtn-dc12 12. Copernicus Land Monitoring Service: Corine land cover 2018 (raster 100 m), europe, 6-yearly – version 2020_20u1, may 2020 (2020). https://doi.org/10.2909/ 960998c1- 1870- 4e82- 8051- 6485205ebbac, https://land.copernicus.eu/en/ products/corine-land-cover/clc2018, european Union’s Copernicus Land Monitoring Service information. Accessed on 01.11.2025. 13. Deo, M., Naidu, C.S.: Real time wave forecasting using neural networks. Ocean engineering 26(3), 191–203 (1998) 14. Eskes, H., Tsikerdekis, A., Ades, M., Alexe, M., Benedictow, A.C., Bennouna, Y., Blake, L., Bouarar, I., Chabrillat, S., Engelen, R., et al.: Evaluation of the copernicus atmosphere monitoring service cy48r1 upgrade of june 2023. Atmospheric Chemistry and Physics 24(16), 9475–9514 (2024)

Air Quality Downscaling with Station-Guided Pseudo-Supervision

17

15. Geiss, A., Silva, S.J., Hardin, J.C.: Downscaling atmospheric chemistry simulations with physically consistent deep learning. Geoscientific Model Development 15(17), 6677–6694 (2022) 16. GISCO, E.: European digital elevation model (eu-dem), version 1.1. Dataset downloaded from Eurostat GISCO website (2016), https://ec.europa.eu/eurostat/ web/gisco/geodata/digital-elevation-model/eu-dem, accessed 2025-11-02 17. Goldsborough III, E., Gopal, M., McEvoy, J.W., Blumenthal, R.S., Jacobsen, A.P.: Pollution and cardiovascular health: a contemporary review of morbidity and implications for planetary health. American Heart Journal Plus: Cardiology Research and Practice 25, 100231 (2023) 18. Gualtieri, G., Brilli, L., Carotenuto, F., Cavaliere, A., Gioli, B., Giordano, T., Putzolu, S., Vagnoli, C., Zaldei, A.: Assessing capability of copernicus atmosphere monitoring service to forecast pm2. 5 and pm10 hourly concentrations in a european air quality hotspot. Atmospheric Pollution Research p. 102567 (2025) 19. Guion, A., Gressent, A., Descombes, G., Janati, Y., Real, E., Ung, A., Meleux, F., Schucht, S., Colette, A.: High-resolution mapping of air quality across europe: an ensemble machine and deep learning framework integrating multi-scale spatial predictors (chromap v1. 0). EGUsphere 2026, 1–35 (2026) 20. Han, T., Guo, S., Ling, F., Chen, K., Gong, J., Luo, J., Gu, J., Dai, K., Ouyang, W., Bai, L.: Fengwu-ghr: Learning the kilometer-scale medium-range global weather forecasting. arXiv preprint arXiv:2402.00059 (2024) 21. Han, W., He, T.L., Jiang, Z., Zhu, R., Jones, D., Miyazaki, K., Shen, Y.: The capability of deep learning model to predict ozone across continents in china, the united states and europe. Geophysical Research Letters 50(24), e2023GL104928 (2023) 22. Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horanyi, A., Munoz Sabater, J., Nicolas, J., Peubey, C., Radu, R., Rozum, I., Schepers, D., Simmons, A., Soci, C., Dee, D., Thepaut, J.N.: Era5 hourly data on single levels from 1940 to present (2023). https://doi.org/10.24381/cds.adbb2d47, https://doi.org/10.24381/ cds.adbb2d47, generated using or contains modified Copernicus Climate Change Service information. Neither the European Commission nor ECMWF is responsible for any use that may be made of the Copernicus information or data it contains. 23. Hossen, M.K., Peng, Y.T., Shao, A., Chen, M.C.: An ode based neural network approach for pm2. 5 forecasting: Mk hossen et al. Scientific Reports 15(1), 24830 (2025) 24. Hsieh, W.W., Tang, B.: Applying neural network models to prediction and data analysis in meteorology and oceanography. Bulletin of the American Meteorological Society 79(9), 1855–1870 (1998) 25. Inness, A., Ades, M., Agustí-Panareda, A., Barré, J., Benedictow, A., Blechschmidt, A.M., Dominguez, J.J., Engelen, R., Eskes, H., Flemming, J., Huijnen, V., Jones, L., Kipling, Z., Massart, S., Parrington, M., Peuch, V.H., Razinger, M., Remy, S., Schulz, M., Suttie, M.: The cams reanalysis of atmospheric composition. Atmospheric Chemistry and Physics 19, 3515–3556 (2019). https://doi.org/10. 5194/acp-19-3515-2019, https://doi.org/10.5194/acp-19-3515-2019 26. Kolehmainen, M., Martikainen, H., Hiltunen, T., Ruuskanen, J.: Forecasting air quality parameters using hybrid neural network modelling. Environmental Monitoring and Assessment 65(1), 277–286 (2000) 27. Koo, J.S., Wang, K.H., Yun, H.Y., Kwon, H.Y., Koo, Y.S.: Development of pm2. 5 forecast model combining convlstm and dnn in seoul. Atmosphere 15(11), 1276 (2024)

18

G. Wang et al.

28. Kuligowski, R.J., Barros, A.P.: Experiments in short-term precipitation forecasting using artificial neural networks. Monthly weather review 126(2), 470–482 (1998) 29. Lam, R., Sanchez-Gonzalez, A., Willson, M., Wirnsberger, P., Fortunato, M., Alet, F., Ravuri, S., Ewalds, T., Eaton-Rosen, Z., Hu, W., et al.: Learning skillful medium-range global weather forecasting. Science 382(6677), 1416–1421 (2023) 30. Lee, Y., Park, J., Kim, J., Woo, J.H., Lee, J.H.: Conditional unet emulation of cmaq simulations for fine particulate matter concentration prediction. Scientific Reports 15(1), 38616 (2025) 31. Liu, D., Grimmond, C., Tan, J., Ao, X., Peng, J., Cui, L., Ma, B., Hu, Y., Du, M.: A new model to downscale urban and rural surface and air temperatures evaluated in shanghai, china. Journal of Applied Meteorology and Climatology 57(10), 2267– 2283 (2018) 32. Lv, B., Hu, Y., Chang, H.H., Russell, A.G., Cai, J., Xu, B., Bai, Y.: Daily estimation of ground-level pm2. 5 concentrations at 4 km resolution over beijingtianjin-hebei by fusing modis aod and ground observations. Science of the Total Environment 580, 235–244 (2017) 33. Lyapustin, A., Wang, Y.: Modis/terra+aqua land aerosol optical depth daily l2g global 1km sin grid v061 (2022). https://doi.org/10.5067/MODIS/MCD19A2.061, https://doi.org/10.5067/MODIS/MCD19A2.061 34. Marzban, C., Stumpf, G.J.: A neural network for tornado prediction based on doppler radar-derived attributes. Journal of Applied Meteorology and Climatology 35(5), 617–626 (1996) 35. McCann, D.W.: A neural network short-term forecast of significant thunderstorms. Weather and Forecasting 7(3), 525–534 (1992) 36. Oliver, M.A., Webster, R.: Kriging: a method of interpolation for geographical information systems. International Journal of Geographical Information System 4(3), 313–332 (1990) 37. OpenAI, R.: Gpt-4 technical report. arxiv 2303.08774. View in Article 2(5), 1 (2023) 38. Organization, W.H., et al.: WHO global air quality guidelines: particulate matter (PM2. 5 and PM10), ozone, nitrogen dioxide, sulfur dioxide and carbon monoxide. World Health Organization (2021) 39. Pesaresi, M., Politis, P.: Ghs-built-s r2023a - ghs built-up surface grid, derived from sentinel2 composite and landsat, multitemporal (1975–2030) (2023). https: //doi.org/10.2905/9F06F36F- 4B11- 47EC- ABB0- 4F8B7B1D72EA, http://data. europa.eu/89h/9f06f36f-4b11-47ec-abb0-4f8b7b1d72ea 40. Pesaresi, M., Politis, P.: Ghs-built-v r2023a - ghs built-up volume grids derived from joint assessment of sentinel2, landsat, and global dem data, multitemporal (1975–2030) (2023). https://doi.org/10.2905/AB2F107A- 03CD47A3-85E5-139D8EC63283, http://data.europa.eu/89h/ab2f107a-03cd-47a385e5 - 139d8ec63283, pID: http://data.europa.eu/89h/ab2f107a-03cd-47a3-85e5139d8ec63283 41. Pu, Q., Yoo, E.H.: Ground pm2. 5 prediction using imputed maiac aod with uncertainty quantification. Environmental Pollution 274, 116574 (2021) 42. Qiu, Y., Feng, J., Zhang, Z., Zhao, X., Li, Z., Ma, Z., Liu, R., Zhu, J.: Regional aerosol forecasts based on deep learning and numerical weather prediction. npj Climate and Atmospheric Science 6(1), 71 (2023) 43. Rautela, K.S., Goyal, M.K., Nagpure, A.S.: Unequal spatio-temporal distribution of population-weighted pollution extremes through deep learning. npj Climate and Atmospheric Science 8(1), 340 (2025)

Air Quality Downscaling with Station-Guided Pseudo-Supervision

19

44. Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: International Conference on Medical image computing and computer-assisted intervention. pp. 234–241. Springer (2015) 45. Rumelhart, D.E., Hinton, G.E., Williams, R.J.: Learning representations by backpropagating errors. nature 323(6088), 533–536 (1986) 46. Schiavina, M., Freire, S., Carioli, A., MacManus, K.: Ghs-pop r2023a - ghs population grid multitemporal (1975–2030) (2023). https : / / doi . org / 10 . 2905 / 2FF68A52 - 5B5B - 4A22 - 8F40 - C41DA8332CFE, http : / / data . europa . eu / 89h / 2ff68a52 - 5b5b - 4a22 - 8f40 - c41da8332cfe, pID: http://data.europa.eu/89h/2ff68a52-5b5b-4a22-8f40-c41da8332cfe 47. Shen, H., Tao, S., Chen, Y., Ciais, P., Güneralp, B., Ru, M., Zhong, Q., Yun, X., Zhu, X., Huang, T., et al.: Urbanization-induced population migration has reduced ambient pm2. 5 concentrations in china. Science Advances 3(7), e1700300 (2017) 48. Shen, H., Tao, S., Liu, J., Huang, Y., Chen, H., Li, W., Zhang, Y., Chen, Y., Su, S., Lin, N., et al.: Global lung cancer risk from pah exposure highly depends on emission sources and individual susceptibility. Scientific reports 4(1), 6561 (2014) 49. Shetty, S., Hamer, P.D., Stebel, K., Kylling, A., Hassani, A., Berntsen, T.K., Schneider, P.: Daily high-resolution surface pm2. 5 estimation over europe by mlbased downscaling of the cams regional forecast. Environmental Research 264, 120363 (2025) 50. Shi, Z., Song, C., Liu, B., Lu, G., Xu, J., Van Vu, T., Elliott, R.J., Li, W., Bloss, W.J., Harrison, R.M.: Abrupt but smaller than expected changes in surface air quality attributable to covid-19 lockdowns. Science advances 7(3), eabd6696 (2021) 51. Sørensen, J.H.: Sensitivity of the derma long-range gaussian dispersion model to meteorological input and diffusion parameters. Atmospheric Environment 32(24), 4195–4206 (1998) 52. Spellman, G.: An application of artificial neural networks to the prediction of surface ozone concentrations in the united kingdom. Applied Geography 19(2), 123–136 (1999) 53. Tangang, F., Hsieh, W., Tang, B.: Forecasting the equatorial pacific sea surface temperatures by neural network models. Climate Dynamics 13(2), 135–147 (1997) 54. Teng, M., Li, S., Xing, J., Fan, C., Yang, J., Wang, S., Song, G., Ding, Y., Dong, J., Wang, S.: 72-hour real-time forecasting of ambient pm2. 5 by hybrid graph deep neural network with aggregated neighborhood spatiotemporal information. Environment International 176, 107971 (2023) 55. Tørseth, K., Aas, W., Breivik, K., Fjæraa, A.M., Fiebig, M., Hjellbrekke, A.G., Lund Myhre, C., Solberg, S., Yttri, K.E.: Introduction to the european monitoring and evaluation programme (emep) and observed atmospheric composition change during 1972–2009. Atmospheric Chemistry and Physics 12(12), 5447–5481 (2012) 56. Varga-Balogh, A., Leelőssy, Á., Lagzi, I., Mészáros, R.: Time-dependent downscaling of pm2. 5 predictions from cams air quality models to urban monitoring sites in budapest. Atmosphere 11(6), 669 (2020) 57. Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, Ł., Polosukhin, I.: Attention is all you need. Advances in neural information processing systems 30 (2017) 58. Wang, F., Tian, D., Lowe, L., Kalin, L., Lehrter, J.: Deep learning for daily precipitation and temperature downscaling. Water Resources Research 57(4), e2020WR029308 (2021) 59. Wang, Y., Fernández-Godino, M.G., Gunawardena, N., Lucas, D.D., Yue, X.: Spatiotemporal predictions of toxic urban plumes using deep learning. PNAS nexus 4(6), pgaf198 (2025)

20

G. Wang et al.

60. Wolf, T., Debut, L., Sanh, V., Chaumond, J., Delangue, C., Moi, A., Cistac, P., Rault, T., Louf, R., Funtowicz, M., Davison, J., Shleifer, S., von Platen, P., Ma, C., Jernite, Y., Plu, J., Xu, C., Le Scao, T., Gugger, S., Drame, M., Lhoest, Q., Rush, A.: Transformers: State-of-the-art natural language processing. In: Liu, Q., Schlangen, D. (eds.) Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: System Demonstrations. pp. 38–45. Association for Computational Linguistics, Online (Oct 2020). https://doi.org/10.18653/v1/ 2020.emnlp-demos.6, https://aclanthology.org/2020.emnlp-demos.6/ 61. Xiao, Y., Wang, Y., Yuan, Q., He, J., Zhang, L.: Generating a long-term (20032020) hourly 0.25° global pm2. 5 dataset via spatiotemporal downscaling of cams with deep learning (deepcams). Science of The Total Environment 848, 157747 (2022) 62. Xie, E., Wang, W., Yu, Z., Anandkumar, A., Alvarez, J.M., Luo, P.: Segformer: Simple and efficient design for semantic segmentation with transformers. Advances in neural information processing systems 34, 12077–12090 (2021) 63. Xue, T., Zheng, Y., Tong, D., Zheng, B., Li, X., Zhu, T., Zhang, Q.: Spatiotemporal continuous estimates of pm2. 5 concentrations in china, 2000–2016: A machine learning method with inputs from satellites, chemical transport model, and ground observations. Environment international 123, 345–357 (2019) 64. Yang, Q., Yuan, Q., Li, T., Yue, L.: Mapping pm2. 5 concentration at high resolution using a cascade random forest based downscaling model: Evaluation and application. Journal of Cleaner Production 277, 123887 (2020)

Air Quality Downscaling with Station-Guided Pseudo-Supervision

21

Supplementary Material A

Additional Data Processing Details

A.1

Side Information Data Sources

The data used to train the model is obtained from various sources, including side information and ground-truth station measurements. A diverse set of geospatial and atmospheric datasets was employed to analyse fine particulate air pollution (PM2.5 ) and its driving factors. These datasets span atmospheric composition forecasts, human settlement and population data, land cover information, topography, satellite-derived aerosol measurements, and meteorological reanalysis. CAMS Global Atmospheric Composition Forecasts. The Copernicus Atmosphere Monitoring Service (CAMS) Global Atmospheric Composition Forecasts provide predictions of atmospheric pollutant concentrations worldwide. The CAMS forecasts are generated using an atmospheric model based on physics and chemistry, which assimilates satellite observations to ensure accuracy. CAMS issues forecasts twice daily, each providing hourly predictions of more than 50 chemical species and seven aerosol types. Notably, the CAMS forecasts include particulate matter outputs such as surface-level PM2.5 (particulate matter with a diameter smaller than 2.5 µm), expressed in units of mass concentration. The horizontal resolution of the global forecast is approximately 0.4◦ ×0.4◦ (≈40×40 km). We utilised the CAMS PM2.5 forecast fields as the air quality input, providing an initial estimate of pollutant concentrations across our study domain. This ensured that our PM2.5 analysis had a consistent global prior informed by both modelling and satellite observation integration, which can be combined with station observations. Global Human Settlement Layer (GHSL). The Global Human Settlement Layer is a collection of spatial datasets supported by the Joint Research Centre (JRC) and the Directorate-General for Regional and Urban Policy (DG REGIO) of the European Commission, together with the international partnership GEO Human Planet Initiative. GHSL produces global spatial information, evidencebased analytics, and knowledge describing the human presence on the planet. These data are useful for incorporating anthropogenic factors into environmental analyses, as they quantify the built environment and population distribution globally. In particular, we employed three GHSL products to capture humanrelated variables relevant to air pollution: – GHS-BUILT-S (Built-up Surface) GHS-BUILT-S depicts the distribution of built-up surfaces, expressed as the number of square metres. It describes the extent of built-up areas, providing the proportion of each cell that is covered by built-up structures. This built-up surface layer helps identify urbanised areas and infrastructure, which are often associated with higher emissions and concentrations of PM2.5 .

22

G. Wang et al.

– GHS-BUILT-V (Built-up Volume) GHS-BUILT-V provides estimates of the distribution of built-up volume (in cubic metres) within each grid cell. It combines information on built-up area and building height to quantify the three-dimensional volume of urban structures. This allows us to account not only for the horizontal footprint of urban areas but also for their vertical extent (e.g., differences between high-rise and low-rise buildings), which may influence local pollution emissions and dispersion. – GHS-POP (Population Grid) The GHSL population grid depicts the distribution of human population as the number of people per cell. This dataset enables analysis of the spatial relationship between population distribution and pollution sources, as well as the diffusion of pollution influenced by human presence. By incorporating population density, our assessment can identify areas where high PM2.5 concentrations coincide with large populations. Together, these GHSL layers enable a comprehensive characterisation of human activities. In our study, they are used to incorporate anthropogenic factors into the PM2.5 analysis. CORINE Land Cover. CORINE Land Cover (CLC) is a pan-European land cover inventory that provides land use and land cover information. CORINE provides a harmonised classification of land cover across Europe with 44 thematic classes, ranging from broad forested areas to individual vineyards. Land cover is an important factor in air quality studies, as it influences emission sources and affects pollutant dispersion—for instance, depending on whether an area is industrial, agricultural, or forested. By using CLC, we can distinguish urban areas from natural land. The inclusion of CLC therefore provides contextual information about the environment in which air pollution occurs, enhancing the robustness of our learning of PM2.5 patterns. EU-DEM (Digital Elevation Model). EU-DEM is a high-resolution digital elevation model covering Europe, representing ground elevations and topography. We included elevation data because terrain influences atmospheric flow and pollution distribution. For example, valleys can trap pollutants, leading to higher PM2.5 concentrations, while elevated or exposed areas may experience greater dispersion. By integrating EU-DEM, our analysis captures the effects of elevation by examining whether high pollution levels coincide with specific topographic features and by including elevation as a predictor in statistical models of PM2.5 . MODIS MCD19A2 Aerosol Optical Depth. To incorporate satellite observations of atmospheric aerosols, we used the MODIS MCD19A2 product. The MCD19A2 Version 6.1 dataset is a Moderate Resolution Imaging Spectroradiometer (MODIS) Terra and Aqua combined Multi-Angle Implementation

Air Quality Downscaling with Station-Guided Pseudo-Supervision

23

of Atmospheric Correction (MAIAC) Land Aerosol Optical Depth (AOD) gridded Level 2 product, produced daily at a 1 km pixel resolution. Specifically, MCD19A2 provides aerosol optical thickness measurements for the blue band (AOD at 0.47 µm) and the green band (AOD at 0.55 µm) for each day. By using this dataset, we incorporated an independent observational perspective on particulate pollution since PM2.5 constitutes a major component of atmospheric aerosols and is therefore correlated with AOD. This is particularly useful in regions or during periods where ground-based measurements are limited. Its daily resolution also captures day-to-day variability in aerosol levels, which is relevant for short-term air quality assessment. ERA5 Reanalysis (Wind Data). Meteorological conditions play a crucial role in air pollution dispersion. Therefore, we incorporated wind data from ERA5, the fifth-generation global atmospheric reanalysis produced by the European Centre for Medium-Range Weather Forecasts (ECMWF). ERA5 provides hourly estimates of a wide range of atmospheric variables at a spatial resolution of approximately 0.25◦ on a global scale. Specifically, we extracted near-surface wind fields (10 m wind components) from the ERA5 single-level dataset. The eastward (u10 ) and northward (v10 ) 10 m wind components describe the ventilation and transport of pollutants: stronger winds generally enhance dispersion, whereas stagnant conditions promote pollutant accumulation. The ERA5 fields provide reliable meteorological inputs corresponding to the temporal and spatial coverage of our PM2.5 data. In our study, we used ERA5 wind speed and direction to investigate how local and regional airflow patterns influence the spatiotemporal variability of PM2.5 . Including ERA5 meteorological data thus provides a dynamic physical context for interpreting PM2.5 patterns, helping to distinguish whether elevated pollutant concentrations result from emission sources or from meteorological stagnation. A.2

Coordinate and Resolution Conversions.

These datasets are provided by different organisations and therefore use different coordinate systems and spatial resolutions. In this study, we adopt a fine resolution of 0.01◦ . Consequently, we need a way to reconcile the differences among datasets. A practical solution is to reproject all the side information to our desired coordinate system and resolution beforehand and then load the standardised data during training. For the station data, OpenAQ provides point measurements that are tied to specific latitude–longitude locations. We also need to project each point onto a 2-D image representation (using the same coordinate system and resolution as the side information) to enable patch-to-patch training. For each patch, we must include the measurement values from all stations falling within that patch, along with the corresponding mask locations. However, including 10 years of data from stations, satellites, and other anthropogenic data sources across Europe is a large-scale data-processing challenge, and speeding up the data-generation process is essential. We therefore

24

G. Wang et al.

distribute the work across multiple CPU cores, assigning each worker a different hour within the same European region so that the workflow runs in parallel.

B

Additional Training Details

B.1

Input Construction and Label Projection

For each timestamp, we construct the complete model input by concatenating the hourly atmospheric variables with the slowly varying geospatial variables, which are typically updated annually. OpenAQ PM2.5 observations are projected onto the common 0.01◦ grid and matched to the corresponding hourly input timestamp; since the side information is at hourly resolution, minute-level differences in OpenAQ timestamps are discarded and only the hour is retained. In the rare case where multiple observations fall into the same grid cell at the same timestamp, their values are averaged and a single record is retained. This has negligible impact in practice, as such stations are usually within about 1 km of each other and report very similar pollution values. The station-level train/validation/test split is then applied to determine the observations available for supervision at each timestamp. B.2

Patch-Based Sampling

Processing the full European domain directly is computationally prohibitive, as loading the entire map into GPU memory is infeasible. We therefore adopt patch-based training with 64 × 64 input patches; at the 0.01◦ (≈1 km) target resolution, each patch covers an approximate spatial extent of 64 km × 64 km. For each timestamp, a station is selected as the patch centre and a 64 × 64 window is cropped around it, with random horizontal and vertical offsets within a controlled range to increase input diversity. All stations located inside the window are used as supervision signals for that patch. Because the stations available at a given timestamp may span a wide geographical area across Europe, we sample multiple station-centred patches from the same timestamp within each epoch, allowing the model to learn from all station locations. Repeatedly sampling the same station is also beneficial, since the random spatial offset introduces local variation and improves robustness. To preserve the large-scale spatial structure of the CAMS prior, we additionally sample background patches without station supervision. For each patch, the corresponding pseudo-label is constructed using different Gaussian kernel designs. Each sample is recorded together with its timestamp and, for station patches, the station coordinates, in the form (\text {timestamp\_id}, \text {station}, \text {station coordinates}) or (\text {timestamp\_id}, \text {background}). This ensures that every station can be sampled at least once as the centre of a training patch, so that all available stations are traversed.

Air Quality Downscaling with Station-Guided Pseudo-Supervision

25

Algorithm 1 Bucketed Timestamp Sampling Require: Timestamp set T , sample sets {St }t∈T , bucket size K Ensure: Ordered mini-batches for one training epoch 1: Randomly shuffle the timestamps in T 2: Partition the shuffled timestamps into buckets {B1 , . . . , BM }, where |Bm | ≤ K 3: for m = 1, . . . , M do 4: Collect all samples associated with the timestamps in the bucket: \mathcal {S}_m \leftarrow \bigcup _{t\in \mathcal {B}_m}\mathcal {S}_t

5: Randomly shuffle the samples in Sm 6: Split Sm into mini-batches 7: Yield the mini-batches from Sm sequentially 8: end for

B.3

Bucketed Timestamp Sampling

Our batch size is large and training is conducted on multiple GPUs, so many timestamp files may be opened simultaneously. Each file contains nine channels over the full European domain and exceeds 400 MB. Although the files are accessed through file handles rather than being fully loaded into memory, opening many large files at once causes heavy I/O pressure, leaving the GPUs underutilised while they wait for data loading. In the worst case, each training step may need to access approximately \text {batch size} \times \text {number of GPUs} different files, which significantly slows down training. Loading all patches from a single timestamp minimises the number of files opened and improves I/O locality, but it sacrifices temporal diversity: PM2.5 is highly time-sensitive and the side information can vary substantially across timestamps. A batch drawn from a single timestamp may cause the model to overfit to a specific time point and may produce gradients that vary significantly between steps, leading to unstable training. Conversely, fully random sampling provides high data diversity but causes extremely slow I/O. Our goal is to balance these two extremes. Inspired by bucket sampling in NLP, where samples of similar length are grouped together to reduce padding, we group patches from a small number of timestamps into the same bucket. Within a bucket, several patches reuse the same timestamp files, so the data only needs to be read from storage once, while later accesses benefit from the file cache. The complete procedure is summarised in Algorithm 1. At the beginning of every epoch, the timestamp order is reshuffled, so that different timestamps are grouped together across epochs. Samples are also shuffled within each bucket, preventing consecutive mini-batches from being dominated by a single station or location. Processing each bucket sequentially restricts data

26

G. Wang et al.

access to a small set of timestamp files, thereby improving I/O locality while retaining temporal variation within each bucket. The bucket size, defined as the number of timestamps in each bucket, controls this trade-off. A size of one maximises I/O efficiency but degenerates into singletimestamp sampling and loses diversity, whereas a very large size approaches fully random sampling. We use a bucket size of K = 4 in all experiments. Preliminary training runs revealed that this choice reduces the training time from around 10 days per epoch under fully shuffled sampling to approximately 1–2 days per epoch. Worked Example: DataLoader Construction. The final data loader follows a bucketed timestamp sampling procedure. Consider a simple example with 12 timestamps and a bucket size K = 4. First, all timestamps are placed into a list and randomly shuffled: [\text {ts}_0, \text {ts}_1, \ldots , \text {ts}_{11}] may become [\text {ts}_7, \text {ts}_2, \text {ts}_{11}, \text {ts}_5, \text {ts}_0, \text {ts}_9, \text {ts}_3, \text {ts}_8, \text {ts}_1, \text {ts}_{10}, \text {ts}_6, \text {ts}_4]. The shuffled timestamp list is then divided into consecutive chunks of size K = 4: \text {Group A} = [\text {ts}_7, \text {ts}_2, \text {ts}_{11}, \text {ts}_5], \text {Group B} = [\text {ts}_0, \text {ts}_9, \text {ts}_3, \text {ts}_8], \text {Group C} = [\text {ts}_1, \text {ts}_{10}, \text {ts}_6, \text {ts}_4]. Next, all samples from the timestamps within the same bucket are collected. For example, suppose each timestamp contains around 500 samples, including station patches and background patches. Group A may contain 520 + 480 + 510 + 490 = 2000 samples in total. These 2000 samples are then randomly shuffled within the group. Finally, all shuffled groups are concatenated sequentially: [\text {shuffled Group A} \mid \text {shuffled Group B} \mid \text {shuffled Group C}]. The PyTorch DataLoader then splits this final index sequence into mini-batches, for example with a batch size of 384. In short, the procedure first shuffles the timestamp order, then groups every K timestamps into one bucket, and finally shuffles all patches within each bucket. The timestamp grouping is different in every epoch, which avoids a fixed sampling pattern. Within each bucket, patch-level shuffling ensures that a batch is not dominated by samples from a single station. Processing buckets sequentially improves I/O locality because consecutive batches tend to access the same small set of timestamp files.

Air Quality Downscaling with Station-Guided Pseudo-Supervision

C

Additional Evaluation Details

C.1

Evaluation Metrics

27

The predicted values at station locations are compared against the corresponding ground-truth observations. The model is optimised using a root mean squared error (RMSE) loss, and the checkpoint with the lowest mean absolute error (MAE) on the spatiotemporally held-out OpenAQ validation set is selected as the best model. Based on the matched prediction–observation pairs, we also evaluate model performance using the coefficient of determination (R2 ). For MAE and RMSE, lower values indicate better performance, as they measure the magnitude of prediction errors. MAE is more robust to outliers, while RMSE penalises larger errors more heavily due to the squared term. The coefficient of determination R2 ranges from −∞ to 1. A value of R2 = 1 indicates perfect prediction, while R2 = 0 means the model performs no better than predicting the mean of the ground-truth values. Negative values indicate that the model performs worse than this baseline. Overall, a good model is expected to achieve low MAE and RMSE and high R2 . C.2

Visualisation Procedure

For the visualisations, we perform sliding-window inference with window size Ws and overlap O, using stride s = Ws − O. Let the input raster be x ∈ RC×H×W and the model be fθ , producing a local prediction for each window. For each window k with top-left corner (hk , wk ), we extract a tile x_k = x[:,\, h_k:h_k+a_k,\; w_k:w_k+b_k], (where ak × bk is the valid window size; boundary windows may have ak < Ws and/or bk < Ws and are reflect-padded to Ws × Ws before inference). The model outputs a local prediction ŷk = fθ (xk ) on the valid region Ωk = {0, . . . , ak − 1} × {0, . . . , bk − 1}. To suppress boundary artefacts, we apply a distance-based fading weight. Let (u, v) ∈ Ωk be local coordinates and define the minimum distance to window edges d(u,v)=\min \!\bigl (u,\; v,\; a_k-1-u,\; b_k-1-v\bigr ). With fading width F = O/2, we use a quadratic fade \omega _k(u,v)= \begin {cases} \left (\dfrac {d(u,v)}{F}\right )^2, & d(u,v)<F,\\ 1, & \text {otherwise}. \end {cases} Overlapping predictions are merged by normalised weighted averaging: \hat {y}(h,w)= \frac { \sum \limits _{k:(h,w)\in \mathcal {R}_k}\omega _k(h-h_k,w-w_k)\,\hat {y}_k(h-h_k,w-w_k) }{ \sum \limits _{k:(h,w)\in \mathcal {R}_k}\omega _k(h-h_k,w-w_k) },

28

G. Wang et al.

where Rk = {(h, w) : (h − hk , w − wk ) ∈ Ωk } is the footprint of window k in the global raster. This design ensures smooth transitions between neighbouring patches and avoids visible discontinuities at patch boundaries.

D

Channel Importance in the Italian Region

For an input patch X ∈ RC×H×W , where the C channels include CAMS PM2.5 and auxiliary geospatial and meteorological variables, a trained CNN predicts Ŷ = fθ (X). Since station observations are sparse, errors are evaluated only at observed pixels Pusing the valid mask M ∈ {0, 1}H×W . The masked MAE is E(Ŷ, Y; M) =

Mp |Ŷp −Yp | p∈Ω P p∈Ω Mp

, and the full-channel baseline error is E0 =

E(fθ (X), Y; M). We compute channel importance using three complementary methods. D.1

Leave-One-Channel-Out Ablation

For each channel c, we construct an ablated input X(−c) by setting the c-th channel to zero while keeping all other channels unchanged: X^{(-c)}_j = \begin {cases} 0, & j = c, \\ X_j, & j \neq c. \end {cases} The ablation importance score is defined as the increase in MAE after removing channel c: I^{\mathrm {abl}}_c = \mathcal {E}(f_{\theta }(\mathbf {X}^{(-c)}), \mathbf {Y}; \mathbf {M}) - \mathcal {E}_0. A positive value indicates that removing this channel degrades the model performance, suggesting that the channel provides useful information. A negative value means that removing the channel slightly improves the local MAE, which may indicate noise, redundancy, or local overfitting with the station-level correction. D.2

Gradient-Based Sensitivity

The second method measures the local sensitivity of the loss with respect to each input channel. We use the masked MSE loss \mathcal {L} = \frac { \sum _{p \in \Omega } M_p \left (f_{\theta }(\mathbf {X})_p - Y_p\right )^2 }{ \sum _{p \in \Omega } M_p }. For each channel c, we compute the mean absolute input gradient: G_c = \frac {1}{|\mathcal {D}|} \sum _{\mathbf {X} \in \mathcal {D}} \frac {1}{HW} \sum _{p \in \Omega } \left | \frac {\partial \mathcal {L}}{\partial X_{c,p}} \right |,

Air Quality Downscaling with Station-Guided Pseudo-Supervision

29

where D denotes the sampled validation patches. The gradient importance is then normalised across all channels: I^{\mathrm {grad}}_c = 100 \times \frac {G_c}{\sum _{j=1}^{C} G_j}. Unlike ablation and permutation, this score does not directly measure the effect on MAE. Instead, it measures how sensitive the model output is to local perturbations in each input channel. D.3

Permutation-Based Importance

The third method evaluates how much the model relies on the spatially and temporally aligned information in each channel. For channel c, we randomly permute this channel across validation samples while leaving the remaining channels unchanged: X^{\mathrm {perm}(c,r)}_{b,j} = \begin {cases} X_{\pi _r(b),c}, & j = c, \\ X_{b,j}, & j \neq c, \end {cases} where b indexes validation samples, πr is a random permutation for repeat r, and r = 1, . . . , R. The permutation importance score is the average MAE increase over R repeats: I^{\mathrm {perm}}_c = \frac {1}{R} \sum _{r=1}^{R} \left [ \mathcal {E} \left ( f_{\theta }(\mathbf {X}^{\mathrm {perm}(c,r)}), \mathbf {Y}; \mathbf {M} \right ) - \mathcal {E}_0 \right ]. The corresponding standard deviation across repeats is also recorded to quantify the variability introduced by random permutation.