Sofia, the capital and largest city of Bulgaria, is located in the western part of the country and encompasses an urbanised area of 245.5 km². With a population of 1,286,965 in 2023, the city has undergone substantial urban expansion and industrialisation in recent decades. This development has led to environmental challenges, including air pollution [12], traffic congestion [13], and associated health risks [14].

To address these issues, the Sofia Municipality has launched several initiatives to enhance the city's infrastructure. For example, the Strategic Integrated Project aims to support local authorities in creating the necessary administrative, organisational, technical, and financial conditions for implementing Sustainable Urban Mobility Plans (SUMPs), facilitating the transition to a climate-neutral and sustainable urban environment [15]. As part of its commitment to climate action, the Sofia Municipal Council has joined the Covenant of Mayors for Climate and Energy for the 2021–2030 period [16]. The municipality has pledged to reduce greenhouse gas emissions by 40% compared to 2007 levels and implement adaptation strategies for climate change. This commitment aligns with national regulations on energy efficiency and renewable energy. The city's action plan integrates measures for climate change mitigation, energy conservation, and the promotion of renewable energy sources to support a more sustainable future. A satellite imagery of the study area, covering Sofia city and the surrounding area, is shown in Fig. 2.

This section presents data collection and preprocessing, including a description of the data sources, selection criteria, and preprocessing techniques used to structure the datasets for analysis. Land use and land cover data (see Fig. 3) are obtained from the Copernicus Urban Atlas (UA) and Corine Land Cover (CLC) for 2018. Aspect and slope images are derived from SRTM (30 m), and building heights are retrieved from the Global Human Settlement Layer - Building Height (GHS-BH, release R2020A). A 24-hour air surface temperature simulation at a 250 m resolution, performed using the WRF model and elevation data, is also employed throughout the study. Open data for the environment [17] and weather [18] were used to complete the input datasets.

UA, available in vector format, provides detailed land use information in Europe, categorised into 17 urban and 10 rural classes. CLC provides raster land-use data with 44 thematic classes in raster format. Using the Raster to Polygon tool from ArcGIS Pro, the CLC data is transformed into vector data.

High-resolution elevation data governs lapse rate cooling by day and cold air pooling at night (see Fig. 4a). GHS-BH, shown in Fig. 4b, is a value-added product obtained by linearly regressing a composite digital elevation model, ALOS World 3D 30 m (AW3D30, 2006–2011) and NASA SRTM 30 m (11–22 February 2000), against building-shadow indicators extracted from Sentinel2 image data composite for 2018 (GHS-composite-S2, R2020A). The result is a dataset that provides building height estimates at 100 m resolution in the Mollweide projection (EPSG:54009).

Derived slope and aspect raster images flag terrain-induced processes that pure elevation misses (see Fig. 5).

The WRF study area spans 59.25 × 59.25 km, subdivided into a 237 × 237 grid, with each cell measuring 250 × 250 meters. These grid cells are merged to form a polygon representing the WRF model's study area. The data collected from nine observation stations, shown in Fig. 6, are categorised by location and urban density (see Table 1). StringMeteo supplies meteorological data collected by Lozen, Orlandovtsi, and Buckton stations [19]. Sensor Community provides meteorological data for Kazichene, Ovcha Kupel 1, Bunkera, Gurmazovo, and Bustmantsi [17]. Data for the Sofia Airport meteorological station is available from the Weather Spark platform [20].

Crowd-sourced measurements from stations were validated for their accuracy and internal consistency from night to day. WRF temperatures were cross-checked against those of the stations for bias estimation, after which the WRF field became the sandbox “truth” on which optimisation experiments operated. Every dataset was reprojected to WGS-84 and resampled, so its cell centres coincide with WRF grid centroids.

For UA, sections of it falling inside the WRF polygon were clipped (intersection), yielding the Intersectioned UA (IUA). Gaps remain where UA provides no coverage. Then CLC was intersected with the same polygon to create the Intersectioned CLC (ICLC), which spans the entire domain. Using a disjoint operation, ICLC parts that do not overlap with the IUA were isolated and subsequently merged with the IUA, producing a seamless land use layer for the entire study area. All land use classes were harmonised and assigned numeric codes. Urban zones were further split by urban fraction ranges into Low Residential (0.01–0.50), Medium Residential (0.51–0.80), and High Residential (0.81–1.00).

The slope and aspect images are processed using OpenTopography’s clipping service. The bounding box of the research area is specified, and only the pixels inside that rectangular grid are retrieved. The portal supplies slope (in degrees) and aspect (0–360°, clockwise from north), which were calculated using a standard 3 × 3 Horn gradient kernel. Outputs were saved as GeoTIFFs. These layers’ Coordinate Reference System is WGS84, which coincides with the WRF model’s one, so no reprojecting is required.

Each clipped layer was then sampled at the WRF cell centroids and exported as a point table (longitude, latitude, data value). The terrain layers were at 30 m, so they were first aggregated to the 250 m WRF grid. The longitude and latitude coordinates of every WRF grid cell centroid were read from a CSV file containing them. Around each centroid, a square polygon measuring 250 m × 250 m was constructed. Half-widths of 125 m were converted to degrees based on Eq. (1) and Eq. (2).

 (latitudinal) (1)

Each 30 m slope/aspect point whose centre fell within a WRF square was identified with a point-in-polygon test. For every WRF cell, the mean of all contained 30 m values was calculated, after discarding no data codes (if there were any). The resulting 250 m aligned fields were written to two new CSV files, preserving the original WRF centroid order to guarantee one-to-one correspondence with the model grid.

Two tiles of the GHS-BH dataset, which included the research zone (R4_C20 and R4_C21), were merged and reprojected to WGS 84 geographic coordinates (EPSG: 4326), aligning horizontally with the WRF model grid. A polygon representing the study area was supplied to the clipping routine, so that only pixels lying inside the research boundary are retained. The clipped raster was saved as a GeoTIFF file. The 100 m building height layer had to be spatially aggregated so that it aligns with the WRF grid resolution of 250 m. Using a centre-in-polygon test, the 100m GHS-BH grid was sampled within each 250 m WRF cell. No data values were excluded, and the arithmetic mean of sampled pixels defined the building height value for that cell. The final result is saved in a CSV file, preserving the centroid order of the WRF grid.

Each land cell is given one row of features, consisting of elevation, slope, aspect expressed as sine and cosine, merged land cover one-hots, a residential intensity index, mean building height, water proximity weights, neighbourhood residential fraction in a ~ 2.7 km radius and two cyclic hour of day variables added at model runtime. Crowd-sourced in-situ observations serve three roles. First, they provide an external bias check on the WRF run before it is adopted as the OSSE truth field. Second, their hourly temperatures are fed into the hybrid Random Forest–Kriging model, which generates a temperature field. Third, the 32 coordinates of these real sites define the geometry of the CoLoc-32 network. By sampling bias-free temperatures directly from the WRF grid for those same points, a synthetic counterpart is obtained, allowing to quantify how much sensor measurement bias influences the results.

Each land use category is one-hot encoded and supplied to the Random Forest, allowing the model to learn, for example, that densely built surfaces retain heat at night, while vegetated pixels cool rapidly. The same classes underpin the initial sensor seeding step. K-means clustering is performed in a mixed feature space that deliberately includes the five urban fractions (from UA) under a single composite variable representing residential area, ensuring that at least one sensor falls within every dominant land use regime. Additionally, a binary water body mask is derived from the data, which removes lakes and rivers from the modelling grid and yields a distance-to-water predictor that captures shoreline thermal inertia. All candidate station coordinates must fall on land pixels.

Elevation appears three times: as a direct predictor in the Random Forest (RF) trend, as the vertical dimension in the Gaussian Process kernel, and as the basis for the low, mid and high elevation strata used when the first twelve hotspot stations are seeded. In addition, steeper slopes signal potential drainage corridors, whereas aspect, which is stored as sine and cosine components to avoid the 360° wrap-around problem, modulates solar exposure and hence daytime warming.

One-hot encoding was applied to the land use data, and the following composite variable is formed (from UA’s urban fabric classes). A Residential Urban Fabric Index (RUFI) summarises residential built-up intensity per grid cell for five Urban Atlas residential classes: continuous urban fabric (> 80%) and discontinuous urban fabric bands (50–80%, 30–50%, 10–30%, < 10%). Let denote the indicator that the cell belongs to the class , with exactly one class active among the five. Assigning monotone weights , the RUFI is defined according to Eq. 3 as follows:

which equals the weight of the active class for that cell. Next, the five land use components that form the upper variable are dropped. Additionally, a second composite variable, Green Vegetation Indicator (GVI), merges the principal vegetation classes, forests, natural grassland, and wetlands, into a single flag, enabling the temperature generation model to recognise cells that are predominantly natural cover. A 500 m circular convolution of the RUFI l variable gives a neighbourhood residential fraction, which is a smooth measure of built-up intensity.

Because lake and river cells are masked in the model, we introduce an auxiliary predictor that indicates the proximity of each land cell to open water. First, the binary water mask contains a total of 744 grid cell centroids whose centres fall inside mapped water bodies. The shortest great circle distance to the nearest water centroid is computed for every land cell centroid using the Haversine formula on a sphere of radius 6371 km. The distances are converted to a unitless proximity weight based on Eq. (4).

with a cut-off .

The sensor network is generated in two separate stages. A deterministic initial set of stations is generated using stratified hotspot sampling combined with multi-feature k-means clustering. This is followed by a stochastic refinement stage, in which any station may be replaced if the substitution increases the overall R², or, with a low probability, if it results in only a minor decrease, facilitating exploration of the solution space. Parameterisation is area-dependent, and no parameter is held strictly constant. Prior to implementing the sampling algorithm, we distinguish two types of station sets used throughout: OBS, the fixed set of existing, quality-controlled observation stations in the domain, and F-OBS, fictive observation station sets produced either by a rule or by the siting algorithm at specified densities. In the OSSE setting, temperatures at F-OBS locations are sampled from the WRF field and used as support for the hybrid RF-GP (Random Forest-Gaussian Process) mapper. Reconstructions are evaluated over the masked grid.

The RF-GP model follows a 3-step procedure. First, the RF model, with hyperparameters optimised for the specific set of OBS/F-OBS (including the number of trees, the number of features considered at each split, the minimum leaf size, and the maximum tree depth), was trained on the engineered predictors to produce a broad-scale estimate (Eq. (5)).

Next, GP residual interpolation is performed. Residuals at the station locations are modelled as a zero-mean GP with covariance (Eq. (6)).

where is the Matern kernel .

Forty random hyperparameter draws are scored on an external validation subset, yielding the best set of horizontal length scales, Matern kernel, and vertical scale.

Finally, the lapse rate is corrected. The grid-wide residual field is added to the RF trend and adjusted back to actual elevation with a standard lapse rate (Eq. (7)).

To ensure baseline representation of the principal thermal regimes, the study area, comprising 55,425 grid cells (with water bodies masked), was stratified into three elevation bands (0–700 m, 700–1200 m, and > 1200 m) and four nocturnal temperature clusters. The clusters were derived using k-means clustering applied to 01:00 AM temperature data, as this variable yielded the most distinct thermal groupings among all hours and mean temperature tested. Within each of the resulting twelve strata, the grid cell exhibiting the highest local nighttime temperature variance (quantified as the standard deviation within a ~ 2.7 km neighbourhood) was selected. Subsequently, a multi-feature k-means procedure was employed to expand the network to the predefined target size 𝑁 (i.e., the desired number of F-OBS stations within the study region), by embedding the remaining grid cells in an eight-dimensional feature space.

where longitude λ, latitude ϕ and height h describe location,

R is the relief over a ~ 2.7km radius (standard deviation of elevation in the neighbourhood),

C is the range of neighbour means in day and nighttime windows,

is the maximum eight-neighbour temperature difference over 24 h,

is the weighted residential fraction index, so built-up areas could be recognised by the model.

A k-means clustering with 𝑘=𝑁−12 was performed to generate the final set of F-OBS, producing the near-optimal network of stations required. Figure 7 illustrates the spatial stratification of the study area into discrete strata defined by elevation bands and nighttime temperature classes, excluding water bodies.

Simulated annealing refinement is performed, starting from the k-means seeds, and a one-for-one swap procedure searches the discrete space of layouts. At each iteration, one station index is replaced by a random non-member, the hybrid RF-GP model is retrained, and the mean hourly validation R² over the full WRF field is re-computed. Improvements are accepted deterministically: declines are accepted with probability . The temperature T decays geometrically each step. The performance of the algorithm on a representative dataset is shown in Fig. 8.

For each candidate network, the hybrid RF-GP model is trained only on that network’s stations and validated on the entire WRF grid, excluding those stations and the water body cells. The primary metric used for validation is the mean hourly coefficient of determination, shown in Eq. (10).

Firstly, the WRF model simulation has been rigorously validated against data from nine observation stations located in the city centre, the suburbs of Sofia, and the surrounding rural areas. The analysis indicates that the WRF model accurately simulates surface air temperatures on August 22, 2018, compared to measurements from observation stations, with biases across the extensive study area ranging from − 0.5 to 0.2°C, except for stations Bukston and Busmantsi, where the bias is approximately 1.0°C. The temperature distribution is illustrated in Fig. 9.

The visualisation of the hourly observations from the OBS network confirms the presence of the UHI between 18:00 and 20:00 LST, which is also captured in the measured data (Fig. 10).

Secondly, the WRF model data, which are treated as reference (“truth”), were compared with the temperature fields generated by the RF-GP model for the first five predefined OBS/F-OBS networks, as well as the three F-OBS networks produced by the sampling algorithm. For clarity, visual comparisons are presented only for three representative hours: 07:00 (coldest), 14:00 (hottest), and 19:00 (when the UHI effect is observed). In addition, scatter plots are provided for each network at each hour to illustrate the correspondence between observed and predicted values. The overall performance across all networks is summarised in terms of and presented in Table 2.

Figure 11 shows the temperature obtained by simulation with the WRF model, while Fig. 12 shows the temperature obtained from Uniform-113.

Table 2 presents the R² scores for the eight defined networks. Across all 24 hours, the performance hierarchy of the networks is consistent. The two 32-station layouts occupy the lowest tier, the three uniform networks perform at an intermediate level, and the three optimised layouts consistently achieve the highest scores across all time slices. The OBS-32 network, which is sparse and concentrated in urban areas, achieves values of 0.70–0.76 between 11:00 and 17:00, when variance is largely explained by slope, aspect, and land-use. However, it collapses to a negative before sunrise due to the near absence of sensors in valleys and ridge crests, where cold air drainage and radiation inversions are most pronounced.

The CoLoc-32 network, which mitigates instrument bias by sampling WRF values at the same sites, exhibits even lower nighttime performance, indicating that spatial placement, rather than measurement error, is the primary limiting factor. Increasing the density of F-OBS in the Uniform-113 network improves nighttime from ~ 0.1 to 0.65 and elevates the 11:00–17:00 window above 0.90. Adding another 104 stations (Uniform-217) provides only a modest mean increase (0.799 → 0.837) and negligible improvement after noon, with the performance curve flattening beyond 217 stations. The Uniform-313 network exhibits a similar daytime skill to Uniform-217, with a marginal improvement of + 0.05 at night.

In contrast, the Optim-113 network, using the same number of stations as Uniform-113, matches or exceeds the performance of Uniform-217 for every hour, demonstrating that feature-guided k-means seeding combined with annealing swaps recovers more thermal structure than uniform densification. Optim-217 provides the best overall results in terms of station count and reaching 0.79 during the most inversion-prone hour (03:00) and 0.94 at the afternoon peak, yielding a 24-hour mean of 0.863 while using one-third fewer F-OBS than Uniform-313. Adding 96 additional F-OBS in Optim-313 yields only marginal improvements (+ 0.01–0.02 at night), highlighting the principle that, once critical gradients are instrumented, further densification offers diminishing returns.

These results suggest that a city or province could potentially reduce the number of sensors by nearly half through strategic placement. A real-world example demonstrated that in the Beijing-Tianjin-Hebei region, removing 194 of 481 weather stations, primarily those clustered on plains, and focusing on underserved mountainous and urban-edge areas resulted in negligible loss of forecast accuracy [21]. The refined 287-station network achieved 97–99% of the improvement provided by the full 481-station network for air quality forecasts. This operational OSSE exemplifies how a sparser, yet optimally distributed network, can rival a much denser one. Our study provides a concrete temperature mapping analogue, showing that Optim-217 ≈ Uniform-313, Optim-113 ≈ Uniform-217, in terms of network performance.

Figure 13 shows the distribution of the mean absolute error (24 hours) by strata for each network.

The temperature results from the two highest-density F-OBS networks, Uniform-313 and Optim-313, are compared and visualised in Fig. 14. Similar observations are made for the other 2 pairs of F-OBS datasets.

Across all elevation bands, the largest discrepancies occur during nighttime hours, particularly in the early morning (01:00–06:00), when atmospheric processes such as radiative cooling, humidity-driven inversions, and wind have a strong influence on temperatures. The absence of dynamic meteorological predictors limits the accuracy of the RF-GP model during these hours, as it relies solely on static spatial data and temperature observations. By contrast, errors are markedly reduced during daylight (11:00–19:00), with MAE dropping to ~ 0.5°C and consistently exceeding 0.9, demonstrating that under clear-sky daytime conditions, static predictors such as elevation, land use, and aspect provide sufficient explanatory power for spatial temperature variability.

The comparison between the Uniform-313 and Optim-313 networks highlights the benefits of strategic sensor placement, particularly under conditions where temporal meteorological factors dominate:

To evaluate how accurately the networks reproduce the WRF reference, the mean temperature difference (WRF minus generated) is mapped for four key phases of the diurnal cycle. The analysis focuses on the network with the highest , namely Optim-313. Figure 15a corresponds to the deep-night period (01:00–06:00), when strong inversions and cold-air drainage dominate. During this interval, the largest errors, reaching up to ± 3°C, occur in valley floors and on exposed ridges, where the predictors lack information on local humidity and wind structure. Figure 15b captures the sunrise transition period (07:00–09:00), during which cold-air pool biases diminish as atmospheric mixing begins, although isolated hollows and sheltered areas remain slightly cooler than observed. Figure 15c represents the broad daylight interval (10:00–19:00), when a well-mixed boundary layer prevails. During this period, the network operates optimally, with temperature differences generally within ± 1°C, except in areas such as Vitosha and north of Sofia, where the Iskar Valley meets the first Balkan ridges. Minor anomalies are also apparent over the urban core and along forested slopes. During the evening transition (20:00–22:00), cold and warm biases reappear in downslope areas and within the urban core as the nocturnal inversion begins to reestablish (see Fig. 15d).

The Optim-313 layout demonstrates that the hybrid model reproduces daytime temperatures accurately in regions of relatively uniform terrain, where the optimised sensor network effectively captures the limited meteorological contrasts. In areas with pronounced relief, however, residual errors indicate that even a strategically placed network cannot fully resolve fine-scale variations in wind, humidity, and radiation that are not captured by the predictors. Thus, while Optim-313 consistently outperforms a uniform grid, further improvements in complex terrain would require either additional targeted sensors or higher-resolution meteorological inputs.

In summary, the optimised F-OBS networks clearly demonstrate the value of strategic sensor placement. Stations are preferentially located in areas that capture the signatures of key processes, such as valley cold pools, slope-driven wind zones, and ridge tops, thereby resolving temperature extremes and gradients that a uniform grid often fails to represent. Consequently, an optimised layout achieves lower errors with fewer stations, particularly during nighttime conditions, when spatial temperature gradients are strongly influenced by unresolved meteorological processes. These results support the central hypothesis that intelligent siting can partially compensate for the absence of high-resolution meteorological data. By embedding sensors in critical locations, the static plus hybrid model can indirectly account for phenomena such as humidity-driven inversions and cold-air drainage, which are not explicitly included in the predictors. In contrast, simply increasing the density of uniformly spaced sensors yields diminishing returns, as many additional stations are placed in already well-characterised areas, leaving important microsites inadequately observed.

This study combined a WRF-based sandbox, a static feature–driven sensor optimisation algorithm, and a hybrid RF-GP interpolation model to assess the potential of strategically placed sensor stations for urban temperature monitoring. The results demonstrate that terrain and land-use proxies capture the majority of spatial temperature variability during daytime and a substantial fraction at night. However, the strategic placement of sensor stations remains essential to resolve the remaining gradients. An optimised network of approximately 200 stations, selected using our algorithm, achieved near-equivalent mapping accuracy to a naïvely dense uniform grid of over 300 sensor stations, highlighting the efficiency and value of information-centric network design. The largest improvements from optimisation occurred in complex terrain and stable nocturnal conditions, where uniform networks exhibited errors of 2–3°C; targeted placement in valleys, transitional slopes, and ridge tops reduced these errors by up to half.

These findings reinforce previous research across diverse climates showing that optimal observing networks frequently require non-uniform layouts. Strategic siting can act as a partial substitute for dense measurements: a carefully designed sensor network effectively “learns” sub-grid meteorological patterns that would otherwise require high-resolution humidity, wind, or radiation data. Nevertheless, there are inherent limits. Even the Optim-313 network, while outperforming a uniform grid, could not fully resolve all sub-grid processes, particularly in the coldest valley regions, suggesting that further improvements would require additional sensors or the inclusion of fine-scale meteorological predictors. From a practical perspective, this study provides city planners with evidence that optimised sensor station deployments can produce higher-quality temperature data for the same investment. Even under resource constraints, carefully placed stations guided by terrain and land-cover analytics can capture critical microclimates, enabling the generation of high-fidelity urban heat maps for public health, energy management, and planning applications. The results show that an optimised layout, focused on key problem areas, can achieve most of the benefits of a much denser network. Open-source workflows, such as the one presented here, combine static GIS data, modelling, and optimisation to support cities in designing cost-effective, data-driven, and strategically placed observation networks, enhancing urban resilience to heat hazards prior to deployment.

Potential improvements could be achieved by enriching the sandbox with downscaled meteorological data. By incorporating kilometre-scale humidity, wind, and radiation fields that are statistically or dynamically downscaled to the 250 m grid, the hybrid RF-GP model could be retrained with these additional predictors, and station placement re-optimised. This approach would clarify whether the inclusion of high-resolution atmospheric information reduces the total number of sensors required or primarily shifts their optimal locations.

Another direction for future work is the deployment of a cost-effective, synthetic sensor layout at a reduced scale. After a full season of operation, measured temperatures could be compared with sandbox predictions (for days where weather conditions correspond to the WRF model, assuming no expansion of input data) and with denser virtual networks. Such an experiment would provide real-world validation of the optimisation framework, accounting for measurement noise and unmodeled variability, and would test whether active learning approaches, which involve adding or relocating a small subset of sensors based on initial error maps, can further enhance coverage without incurring substantial additional costs.