An Application of the LCZ Approach in Surface Urban Heat Island Mapping in Sofia, Bulgaria

Abstract

This article presents the results of the thermal survey of the capital of Bulgaria (Sofia) carried out in August 2019, with the application of an unmanned aerial system (UAS). The study is based on the concept of local climate zones (LCZs), taking into account the influence of the features of land use/land cover and urban morphology on the urban climate. The basic spatial units used in the study are presented in the form of a regular grid consisting of 3299 cells with sides of 250 × 250 m. A total of 13 types of LCZs were identified, of which LCZs 6, 5, 8, 4, D, and A form the largest share. In the thermal imaging of the surface, a stratified sampling scheme was applied, which allowed us to select 74 cells, which are interpreted as representative of all cells belonging to the corresponding LCZ in the urban space. The performed statistical analysis of the thermal data allowed us to identify both the most thermally loaded zones (LCZs 9, 4, and 5) and the cells forming Urban Cool Islands (mainly in LCZs D and C). The average surface temperature in Sofia during the study period (in the time interval between 8:00 p.m. and 10:00 p.m.) was estimated at 20.9 °C, and between the different zones it varied in the range 17.2–25.1 °C. The highest maximum values of LST (27.9–30.6 °C) were registered in LCZ 4 and LCZ 5. The relation between the spatial structure of the urban thermal patterns and urban surface characteristics was also analyzed. Regression analysis confirmed the hypothesis that as the proportion of green areas increases, surface temperatures decrease, and, vice versa, as the proportion of built-up and impermeable areas increases, surface temperatures increase. A heat load map (via applying a z-transformation to standardize the temperature values), a map of the average surface temperature, and a map of the average intensity of the heat island on the surface were generated in the GIS environment. The results of the study adequately reflect the complex spatial model of the studied phenomenon, which gives grounds to conclude that the research approach used is applicable to similar studies in other cities.

1. Introduction

The urban heat island (UHI) is one of the most significant examples of the impact of cities on the environment. This phenomenon was first documented in 1818 for the city of London by the “father of urban climatology” Luke Howard [1,2,3].

An UHI is characterized by positive temperature differences between the city and its surroundings, which usually reach their maximum a few hours after sunset [4]. As a traditional measure of the intensity of the phenomenon, the difference between two measurement locations is used—one representing “urban” conditions and the other—“rural” [5]. The magnitude of the temperature difference (ΔT_u–r) correlates positively with the size of the city [6], but there are studies that do not establish such a relationship [7]. The study of UHIs is not only of scientific interest, but also has important practical implications for urban planning. The negative effects of an UHI are usually associated with a general deterioration of the urban environment, with increased health risks, increased consumption of energy for cooling, the frequent formation of smog, an increased concentration of ground-level ozone, and others [8,9,10,11]. The study of the genesis and spatial–temporal characteristics of this phenomenon is necessary in developing strategies to mitigate the negative impacts at the meso-, local, and micro-level and to take actions to reduce the anthropogenic heat load of cities.

Cities are characterized by a high level of spatial heterogeneity that results in varied microclimate effects [12]. Therefore, understanding the effects of urbanization on an UHI requires us to assess the urban structure (e.g., the height-to-width ratio of buildings and streets), surface cover (e.g., the proportion of built-up vs. vegetated surfaces per unit area), urban fabric (the materials used to construct cities), and urban metabolism (anthropogenic heat, water, and pollution), as they modify local energy budgets and water balances, forming heat islands [13,14,15,16]. The urban thermal environment varies not only from its rural surroundings but also within the urban area due to intra-urban differences in land-use and surface characteristics, which led to a rethinking of the traditional dichotomy “urban–rural” in the study of the phenomenon [2,17,18,19,20,21,22]. The UHI cross-section usually distinguishes between “cliff”, “plateau”, and “peak”, where the highest temperatures are recorded [23]. However, this general model of the spatial structure of an UHI is often disturbed due to the uniqueness of the spatial structure in the land use/land cover (LULC) of each city, creating local variations, including conditions for the formation of Urban Cool Islands (UCIs), usually in parts of cities saturated with Blue–Green Infrastructure (BGI) [24,25,26] or due to shading by tall buildings [27].

Depending on the place of formation, different types of UHI are distinguished, the observations of which are performed by different methods. An UHI can be observed on the surface of the city (surface heat islands), below it (subsurface heat islands), and above it (atmospheric heat islands) [2]. Atmospheric heat islands are most often studied by observing the air temperature (T_a) from stationary or mobile meteorological stations. These data have a good temporal resolution but have a poor spatial resolution due to the low density of the monitoring network in most cities. It also appears to be problematic to provide data from reference “rural” stations for comparison with data from “urban” stations because, due to urbanization, even stations from suburban airports are already in an urban environment [28]. Subsurface heat islands are less frequently studied, using data from temperature measurements at different depths in the soil/substrate or in the groundwater bodies of cities [29,30]. Surface heat islands (SUHI) as well as surface cool islands (SUCI) are studied using land surface temperature (LST) or radiant temperature data obtained by remote sensing of the radiation and thermodynamic properties of urban surfaces (horizontal, oblique, or vertical), recorded by thermal infrared (TIR) sensors mounted on satellite, airplane, or unmanned platforms or helicopters [31,32]. The low temporal resolution of the satellite data limits our ability to continuously estimate the intensity of the UHI. Usually, under the same meteorological conditions, the surface and ground air temperatures do not coincide because the remote sensor captures the long-wave radiation from the surface and not the air temperature read by thermometers, which forms the thermal environment in the inhabited atmospheric layer (the urban canopy layer, UCL). Due to the differences in the processes forming the energy budgets of the two types of media, LSTs vary in much wider ranges than T_a [33,34,35]. As the UHI is a multiscale phenomenon, combined methods are increasingly used in its study, which allows us to cover micro-, local, and mesoscale processes affecting the magnitude and spatial–temporal characteristics of different types of UHI [2].

The idea of dividing cities into “climate zones” has a long tradition in urban climatology [36]. Since UHIs depend mainly on the specific features of each city, spatial differentiation (zoning) of urban space is usually performed according to some classification scheme so that the intensity of the UHI can be related to certain factors and types of urban structures (intra-urban variations of LULC, urban morphology, degree/type of human activities, and other criteria) [37,38,39,40]. Based on previous classification schemes [15,41,42,43], Stewart and Oke developed a new system of local climatic zones (LCZs) [44,45,46], which is the most common scheme for spatial differentiation of urban areas for climatic purposes. The classification divides urban and rural landscapes into 17 standard classes, each defined by ten structural or “built-up” (LCZs numbered from 1 to 10) and seven land cover or “natural” (LCZ named with letters from A to G) climate-relevant surface properties that influence the air temperature. Thus, according to the building types, LCZs from 1 to 3 are compact high-rise, compact midrise, and compact low-rise; LCZs from 4 to 6 are open high-rise, open midrise, and open low-rise; and LCZs from 7 to 10 are lightweight low-rise, large low-rise, sparsely built-up, and heavy industry. Based on seasonal fluctuations, it is possible to differentiate additional variable subclasses of LCZs (b—bare trees, s—snow cover, d—dry ground, and w—wet ground). LCZs are regions of uniform surface cover, structure, material, and human activity spanning hundreds of meters to several kilometers on a horizontal scale. In the context of the LCZ classification system, the urban heat island intensity is not an urban–rural near-surface air temperature difference (ΔTu—r), but an air temperature difference between pairs of LCZ types (ΔT_{LCZ X–Y}; that is, an inter-zone temperature difference [47].

Available spatial databases on urban structures are used as primary sources of information. This includes digital cadastre data, remotely received LULC and LST data, orthophotos, DEM/DTM data and building geometry (via LiDAR- and SAR-based technologies), data from ground stations for T_a, and freely available data from OpenStreetMap, Google Earth Engine, Google Street View, etc. Primary methods for LCZ mapping include in situ measurements, geographic information system (GIS)-based and remote-sensing-image-based calculations (the WUDAPT method), or a combination of them, which are presented in detail in many publications [48,49,50,51,52,53,54,55,56]. To improve the accuracy, some scientists offer an alternative GIS-based workflow to translate Copernicus datasets from Urban Atlas and Corine Land Cover shapefiles and Raster GeoTiff into LCZ maps, which, however, limits the possibilities for practical application of the method outside major European cities [57]. Testing of the methodology for five southern European cities (Athens, Barcelona, Lisbon, Marseille, and Naples) shows promising results [58]. At the same time, methods for GIS-based spatial analysis and UHI map generation are being improved [59], as well as the interpolation techniques for cartographic representation [60]. Correlation analysis and regression modeling are commonly used to detect the influence of LULC and urban morphology on the UHI [61,62]. Thanks to advances in geospatial technologies, UHI/SUHI maps based on the LCZ concept have been created for hundreds of cities around the world today. There are also examples of the application of the LCZ system in studies for which it was not originally intended (e.g., for mapping and assessment of urban ecosystem conditions and services [63] and urban ecohydrological studies [64]).

Although remote thermal studies do not provide direct observations of T_a in the lower layer of the urban atmosphere, they give an idea of SUHI, which is the main engine of atmospheric UHI. The various methods for remotely acquiring LST data have both advantages and disadvantages. The main advantage of using data from satellite sensors is that LSTs can be observed simultaneously over a large area and temperature variations can be analyzed depending on the influence of certain factors (LULC and urban morphology). Using appropriate satellite data depends on the individual sensors and their spectral and time resolution. Still, the optimal time to receive data may not coincide with the time the satellite passes over the target object because the temporal resolution is not controlled by the user [65].

Unlike satellite data, data from airborne thermal remote sensing (from airplanes and helicopters) have a much higher level of detail, which allows for the study of microscale variations in LST between the individual elements of the urban landscape. Another significant advantage is that you can choose the most appropriate time of day to receive thermal data. However, the disadvantages are the high costs, the complex organization, and the strict rules of flight management in a country. This type of observation is usually combined with ground temperature measurements from stationary and/or mobile meteorological stations [66,67,68,69].

An alternative approach is the use of unmanned aerial vehicle systems (UAVs, or drones), which are increasingly used in various scientific fields [70]. These platforms can be equipped not only with thermal, but also with multispectral and hyperspectral sensors, as well as with LiDAR, which allows for obtaining very detailed information about the ground cover and urban morphology needed to generate 2-D and 3-D models of the urban environment. This is important for the assessment of the heat load of cities, as the total area of the different types of heat-emitting surfaces in them is much larger than in the suburbs [71]. In recent years, there have been many examples of successful use of UAV-based thermal studies in different parts of the world. For instance, in several studies, Gaitani et al. [72,73,74] presented the results of a thermal survey in a suburb of Athens, and Naughton & McDonald [75] used a UAV quadcopter drone and thermal camera at two university campuses in Milwaukee (Wisconsin) and El Paso (Texas). However, we are not aware of any UAV-based thermal studies of LCZs that have been performed so far throughout the whole city. The main advantage of thermal studies with an UAV is obtaining high-quality data for LST with a very high spatial resolution in a flexible, easy, and accessible way. On the other side, the application of UAVs often is affected by a number of limitations related to urban safety regulations and privacy [76].

Collecting air temperature data with an appropriate spatial resolution for an entire city like Sofia is a complex task that requires expensive equipment and labor-intensive organization. A possible approach is to use remote sensing to provide a snapshot of the LST in time as a proxy for air temperature. This publication presents some of the results of a more extensive study using GIS-based spatial analysis methods, traditional satellite data sources, and high spatial resolution data from unmanned aerial vehicle systems (UAVs) [77]. The paper is organized as follows. After reviewing the literature on the topic (this Section), we present the approach used in the development of the map of LCZs in Sofia (Section 2.2) and UAV-based thermal surveys performed using a stratified sampling scheme (Section 2.3). In Section 3.1, we present the statistical characteristics of LST by LCZs, in Section 3.2 we perform spatial interpretation and analysis of the results obtained from the thermal imaging of the surface, and in Section 3.3 we demonstrate the map of the intensity of SUHI in Sofia. In Section 4, we discuss the results of the research, and in Section 5 we draw some conclusions about future research.

2. Materials and Methods

2.1. Study Area

The present study is focused on the urbanized area of Sofia—the capital of Bulgaria. The city is located in the largest valley in Bulgaria—Sofia (about 1200 sq. km) at 550 m above sea level. The relief of the valley floor is flat, with a slight slope up to the north. To the north, it is surrounded by the Balkan Mountains and to the south—by the Vitosha Mountain. The climate is temperate continental with four clearly distinguished seasons. According to the climatic classification of Köppen–Geiger, the study area has a climatic index Cfb (C—warm temperate, f—fully humid, b—warm summer), [78,79]. According to the Sofia Station, which has been in operation from 1881 to the present, the average annual temperature is 10.3 °C, and the average yearly rainfall is 612 mm, with a maximum in May–June. The average temperature of the coldest month is –1.7 °C and the warmest month is 21.2 °C. The measured absolute minimum air temperature is –27.5 °C, and the absolute maximum temperature is 37.4 °C (Figure S1). Sofia’s urban area is a complex combination of different functional zones (administrative, residential, industrial, transport, commercial, parks, etc.), which have been developed at different times, over a very long period.

2.2. Local Climate Zones (LCZs)

For the study, an integrated spatial model of local climate zones (LCZs) was created, which is based on the classical concept developed by Oke and Stewart. In the present study, preliminary research and fieldwork were performed to determine the location of specific areas within a given LCZ. Metadata sources include information from aerial imagery with an UAV, orthophotos, a DEM, LULC maps, cadastral maps, street maps, Copernicus datasets from Urban Atlas, and data from Google Earth and Google Street View, which were integrated and organized into a target geodatabase. This, in turn, made it possible to determine the general characteristics of the urban landscape and to specify the spatial distribution of different types of local surfaces (the impervious areas and the Blue–Green Infrastructure elements) and the peculiarities of the urban morphology (e.g., street center-lines and width on the streets, building footprints, the density and height of construction).

Usually, the basic spatial units for which LCZ information is collected are regular grids (with a resolution of 100–500 m) [80,81], lot area polygons [82], urban blocks [83], or urban function zones (UFZs) [84]. For this study, a regular grid was used, consisting of 3299 cells with sides of 250 × 250 m, each cell covering an area of 6.25 ha. In this way, both the unification of the range of the basic units, which are the subject of research, is achieved, and the application of the geostatistical and spatial-interpolation methods in the further analysis is facilitated. Based on similarities in land use and morphological characteristics, most of these cells were grouped into different LCZ classes. A basic rule in determining the affiliation of each cell to a particular LCZ is that of the predominant area of its surface (majority rule), i.e., if different morphological and land-use types fall within the scope of a cell, the cell is determined by the predominant type. This approach, although leading to a greater degree of generalization, allows for a better adaptation of spatial data for the purposes of stratified sampling in subsequent geostatistical processing. Based on orthophoto images and three-dimensional models of the surface in the city of Sofia provided in the Google Earth Pro platform, as well as with the help of Google Street View, a detailed analysis of the more complex areas (cells) and validation of already defined group cells were performed. In determining the height of the building, information on the height of buildings with an Urban atlas source with a resolution of 10 m was used, as the height information is based on IRS-P5 stereo images and is presented as a digital surface model (Figure S2). With the results of this analysis, a map of LCZs was prepared, and, subsequently, some previously unconfirmed areas were validated by the field team. Thus, 13 types of LCZ were identified within the urban area (Figure 1). Built classes cover around 71.1% of the total study area, compared with 28.9% of the non-built classes. Most classes are in the LCZ types 6, 5, 8, 4, D, and A. The share distribution for each type of LCZ is presented in

Table 1. Local climate zones (LCZs).

LCZ Class	Number of Cells	Surface (ha)	% of Study Area	Sample Cells Per Stratum (n)
LCZ_3	26	162.5	0.8	1
LCZ_4	346	2162.5	10.5	7
LCZ_5	636	3975.0	19.3	15
LCZ_6	723	4518.8	21.9	12
LCZ_8	547	3418.8	16.6	14
LCZ_9	34	212.5	1.0	2
LCZ_10	34	212.5	1.0	2
LCZ_A	216	1350.0	6.5	4
LCZ_B	117	731.3	3.5	2
LCZ_C	194	1212.5	5.9	2
LCZ_D	298	1862.5	9.0	10
LCZ_E	105	656.3	3.2	2
LCZ_G	23	143.8	0.7	1
Total	3299	20,618.8	100	74

.

2.3. UAV-Based Thermal Survey of LCZs

LST data within the urbanized space of Sofia were collected in the field for selected representative cells in the grid using an UAV with an installed Thermal Infrared (TIR) sensor. In this case, the Albris platform of the Swiss company Sensefly was used, which is a V-shaped inspection quadcopter with five sensors and three camcorders, allowing for three types of imaging with image detail below 1 cm depending on the subject of study. The TIR camera has parameters of 80 × 60 pixels, with a 50-degree angle and a fixed measuring point, positioned on a three-axis head (gimble). The thermal imager is an uncooled microbolometer type with a thermal sensitivity of 150 mK and an accuracy of ±3 °C or ±5% of reading. The unmanned aerial system enables remote capture of thermal images and video with a dynamically adaptable scale in real-time (Figure 2).

Stratified sampling (STR) was applied to the selection of representative cells, in which the population comprising all 3299 cells was first divided into non-overlapping subpopulations (13 LCZ classes), called “strata”, and sampling was performed independently in each of the strata. As the individual areas of LCZs vary widely, the weights of the areas whose total sum is equal to 1 were also determined; for example, LCZ 6 received a weight of 0.219, as it occupies the most significant space of all zones, and LCZ G recevied a weight of 0.007 because it has the smallest area compared with all other sites. When transmitting the required number of cells (samples) in each zone, the sample size was determined by the method of irreversible selection [85,86,87]. An essential advantage of the stratified sample, compared with other sampling schemes (simple random, systematic, clustered, etc.), is that with the same sample size the results have less stochastic error and higher accuracy due to less scattering in the strata compared with the total scattering in the population [88,89,90]. As a result of the statistical approach used, a total of 74 cells belonging to different local climatic zones were identified from all 3299 cells (Table 1). The locations of these cells play the role of “virtual meteorological stations”, and the measured LSTs in the respective cells are interpreted as representative of all cells belonging to a given LCZ in the urban space.

The most appropriate study period of the year is August, which is the month with the highest average temperature of the Earth’s surface, which is a fundamental prerequisite for the most visible signs of the UHI phenomenon. Its effect is most clearly observed during the early evening hours, immediately after the active sunshine, with a relatively windless and cloudless sky in the situation of anticyclonic synoptic conditions. In this part of the day, there is the best opportunity to maximize the potential for precise identification of micro and local climatic differences between urbanized territories and their response to the UHI. Regarding this, the time interval between 8:00 p.m. and 10:00 p.m. was selected for thermal imaging observations for the 74 predefined grid cells. The UAV-based thermal imaging campaign was conducted over eight days on calendar dates 7, 8, 12, 13,14, 21, 22, and 23 August 2019.

In each of the studied cells, there are different types of land cover and types of construction, which are in different proportions. Therefore, in the range of each cell, the areas of four main categories of surfaces were quantified: (1) Sealed surfaces; (2) Buildings; (3) Vegetation; and (4) Water bodies, that were used in the calculation of the weighted average LST value for the given thermally examined cell.

3. Results

As a result of the research, a general geoinformation model of LCZs of the urbanized space of Sofia was created in the form of a standardized grid of size 250 × 250 m and various types of attributive information obtained through the integrated application, in situ observations and surveys, mapping of the territory via an UAV, and classical remote sensing and geostatistical and geospatial analytical operations.

3.1. LST Statistics by LCZs

For each LCZ, the statistics of the LSTs measured in them were processed, with the exception of LCZ 3 and LCZ G, which are represented by only one sampling cell and contain only one temperature value. Of all 74 cells, 58 cells covered 77.3% of the total area of all cells examined. These are LCZ 4 (7 cells), LCZ 5 (15 cells), LCZ 6 (12 cells), LCZ 8 (14 cells), and LCZ D (10 cells).

Table 2. LST statistics by LCZs.

LCZ Class	Min	Max	Mean	Median	Stand. dev	25 Prcntil	75 Prcntil	Skewness	Kurtosis	Coeff. var	Std. Error
LCZ 3			20.6
LCZ 4	18.8	30.6	23.9	23.8	4.3	19.1	26.6	0.3	−0.6	17.8	1.6
LCZ 5	19.0	27.9	22.5	22.5	2.4	20.7	23.7	0.7	0.4	10.5	0.6
LCZ 6	17.1	23.3	20.1	20.2	2.2	17.8	21.9	0.1	−1.4	10.8	0.6
LCZ 8	17.1	24.8	21.3	22.0	2.6	18.4	23.9	−0.3	−1.5	12.4	0.7
LCZ 9	23.7	26.4	25.1	25.1	1.9	17.8	19.8	0.0	−2.8	7.6	1.3
LCZ 10	21.0	23.5	22.3	22.3	1.8	15.7	17.7	0.0	−2.8	8.2	1.3
LCZ A	17.0	20.2	19.0	19.4	1.5	17.5	20.2	−1.0	−0.7	7.9	0.8
LCZ B	21.2	21.3	21.3	21.3	0.1	15.9	16.0	0.0	−2.8	0.4	0.1
LCZ C	16.8	17.6	17.2	17.2	0.6	12.6	13.2	0.0	−2.8	3.6	3.6
LCZ D	14.4	22.6	17.8	18.1	2.7	15.0	19.7	0.3	−0.7	14.9	0.8
LCZ E	19.6	21.9	20.7	20.7	1.7	14.7	16.4	0.0	−2.8	8.1	1.2
LCZ G			17.8
All cells	14.4	30.6	20.9	20.7	3.1	18.6	23.1	0.4	0.4	14.9	0.4

presents the statistical characteristics of the LST of the studied cells.

The highest mean and median values of LST (25.1 °C) were registered in LCZ 9, but due to the insignificant relative share of the zone (~1.0%) its contribution to the total heat load of the city is minor. On the other hand, the highest maximum values of LST (27.9–30.6 °C) were observed in LCZ 4 and LCZ 5, which occupy a significant share (about 30%) of the urban area. As expected, the lowest temperatures are measured in the “natural” classes of climatic zones (LCZ A—LCZ G). Intra-LCZ variability in the measured LSTs is best seen in the Box-Plot diagrams of Figure 3. It can be seen that the most extensive range of LST is in LCZ 4 (11.8 °C) and in LCZ 5 (8.9 °C), which can be explained by the more complex and diverse internal spatial structure of these areas—the presence of complex combinations of anthropogenic elements (buildings with different stores and geometries and differences in the type and area of sealed spaces) and natural features (green spaces with diverse vegetation and/or bare permeable soil substrates).

3.2. Spatial Interpretation and Analysis of Field Research Data

In this part of the study, we present the results of the performed z-transformation of the temperature values, as well as the generated map based on this transformation, showing the spatial distribution of the different degrees of heat load in Sofia. Then, we present the constructed regression model, revealing the relationship between the main types of land cover and the measured temperatures in them. Finally, we present the map of the spatial distribution of temperature values.

The data for LST were used for the preliminary estimation of the heat load in Sofia. The heat load can be categorized into four stages (most favorable, favorable, less favorable, and most unfavorable) on the basis of the standardized temperature values by the applied z-transformation according to the following formula [91]:

z = (X − µ)/σ,

(1)

where X is the mean temperature of a single test cell, µ is the mean temperature of all cells tested, and σ is the standard deviation of the mean LST of all cells.

Table 3. Thermal load of urban areas.

Heat Load Categories	z-Scores	Area (%)	Population (%)
1 (most favourable)	<−1.00	16.22	20.3
2 (favourable)	from −0.99 to 0.00	35.14	38.1
3 (less favourable)	from 0.01 to 0.99	36.49	30.2
4 (most unfavourable)	>1.00	12.16	11.4

and Figure 4 present the results of the spatial z-transformation of the temperature values.

The areas with values of z <−1 are distinguished by the lowest levels of heat load. They cover 16.22% of the urbanized area and include mainly forested areas, areas with grass or shrubs and surface water bodies (LCZ A, LCZ C, LCZ D, and LCZ G), as well as low-rise areas with yards and gardens, with trees scattered between the buildings (LCZ 6). The lowest value of the z-score in this category is registered in LCZ D (z-score = −2.10). Category 4 covers 12.16% of the urbanized area and shows the highest degree of heat load (hot island) with a value of z of over 1 (i.e., areas with a LST higher than the mean + one std). These are mainly areas with high and medium–high construction (LCZ 4 and LCZ 5), as well as large low-rise buildings with extensive sealed spaces between them (LCZ 8). Even single cells with z-score = 3.11 (LCZ 4) and z-score = 2.24 (LCZ 5) fall into this category. The other categories (2 and 3) occupy the largest share of the urbanized territory (71.62%).

The results from the in-situ UAV-based thermal surveys were used to construct linear regression dependencies through the least-squares method between the areas of the main land-use types (as an independent quantity) and the average measured temperatures by land-use types (as a dependent value). The general conclusions confirm the hypothesis that the greater the proportion of green space, the lower are the measured temperature values, and vice versa—the higher the proportion of built-up and impervious areas, the higher are the surface temperatures (Figure 5).

Using the created regression dependencies and the spatial distribution of the measurements made by the UAV, a geostatistical model of the spatial distribution of the surface temperature within the city of Sofia was created. For this purpose, an approach for deterministic local interpolation was applied whereby the predicted values are generated within local groups of adjacent points based on the network of measured values and their locations via the UAV. The result of spatial interpolation clearly shows that within the urban space in Sofia a poly-structural (polycentric) urban heat island has been formed (Figure 6).

The great diversity of the territories involved in the different LCZs, as well as the high spatial variability in their spatial structure, determines to a large extent the complex and poly-structural nature of the UHI phenomenon. For example, only one section of 1000 square meters exhibits significant temperature contrasts. The figure below shows a fragment of a field measurement of the surface temperature at 10:00 p.m., carried out in the area of a large hypermarket in Sofia. The image clearly shows that within this confined space, substantial differences in surface temperature are observed—in the grassland, it is within 12.9–13.1 °C, in the parking places made of concrete rosette panels it is 17.1 °C, and on the asphalt road it is 23.6 °C (Figure 7).

3.3. Map of SUHI Intensity

The observed temperature differences characterize the complexity of the UHI and represent one of its most important parameters—its intensity. The intensity of SUHI represents the temperature difference (ΔT) between two local climatic zones x and y (ΔT_{LCZx – LCZy}). Inter-LCZ differences in measured LST in Sofia are presented in

Table 4. Temperature differences between LCZs necessary for calculation and mapping of SUHI magnitude.

LCZ		3	4	5	6	8	9	10	A	B	C	D	E	G
	Average LSTs (°C)	20.6	23.9	22.5	20.1	21.3	25.1	22.3	19.0	21.3	17.2	17.8	20.7	17.8
3	20.6		−3.28	−1.94	0.53	−0.71	−4.46	−1.66	1.56	−0.67	3.39	2.77	−0.15	2.77
4	23.9	3.28		1.34	3.81	2.56	−1.16	1.62	4.83	2.61	6.67	6.05	3.13	6.05
5	22.5	1.94	−1.34		2.47	1.22	−2.56	0.28	3.49	1.27	5.33	4.71	1.79	4.71
6	20.1	−0.53	−3.81	−2.47		−1.24	−4.96	−2.19	1.03	−1.2	2.86	2.24	−0.68	2.24
8	21.3	0.71	−2.56	−1.22	1.24		−3.76	−0.95	2.27	0.04	4.1	3.48	0.56	3.48
9	25.1	4.46	1.16	2.56	4.96	3.76		2.76	6.06	3.76	7.86	7.26	4.36	7.26
10	22.3	1.66	−1.62	−0.28	2.19	0.95	−2.76		3.22	0.99	5.05	4.43	1.51	4.43
A	19.0	−1.56	−4.83	−3.49	−1.03	−2.27	−6.06	−3.22		−2.23	1.83	1.21	−1.71	1.21
B	21.3	0.67	−2.61	−1.27	1.2	−0.04	−3.76	−0.99	2.23		4.06	3.44	0.52	3.44
C	17.2	−3.39	−6.67	−5.33	−2.86	−4.1	−7.86	−5.05	−1.83	−4.06		−0.62	−3.54	−0.62
D	17.8	−2.77	−6.05	−4.71	−2.24	−3.48	−7.26	−4.43	−1.21	−3.44	0.62		−2.92	0
E	20.7	0.15	−3.13	−1.79	0.68	−0.56	−4.36	−1.51	1.71	−0.52	3.54	2.92		2.92
G	17.8	−2.77	−6.05	−4.71	−2.24	−3.48	−7.26	−4.43	−1.21	−3.44	0.62	0	−2.92

.

Based on the data from the matrix in the table, a statistical surface was created using local polynomial interpolation, representing in a continuous form the UHI magnitude in the urban space of Sofia (Figure 8).

4. Discussion

The study clearly showed that the “artificial” LCZ classes have significantly higher surface temperatures than the “natural” LCZ classes, which are relatively cooler. It was also found that there are significant temperature differences not only between the individual LCZs, but also within the zones themselves. Therefore, future research should focus on revealing the micro-scale effects caused by the influences of the type of construction, the orientation of the buildings, the ventilation possibilities, the urban canyons, the used building materials, the peculiarities of the green infrastructure, the meteorological conditions and several other factors affecting the climate in Sofia.

The range and spatial configuration of SUHI show significant variability over relatively short distances due to the complex combination of anthropogenic and natural elements forming the urban landscape. Therefore, the SUHI looks more like an “archipelago” than a clearly defined single thermal island. However, the warmest places generally coincide with the residential areas, which also affects the climatic comfort of the people. Therefore, future micro-scale studies should focus mainly on these urban areas.

5. Conclusions

The results of the study adequately reflect the complex spatial model of the studied phenomenon, which gives grounds to conclude that the research approach used is applicable to similar studies in other cities. The developed geoinformation model of LCZs and the related integrated database constitute a suitable basis for further research related to mapping and assessing the effects of the urban heat island. Because of advancing urbanization and deepening climate change, urban areas will inevitably become the focus of such research. The approaches, methods, and technological solutions used in the present study were fully consistent with the challenges that objectively exist in the study of UHIs, as well as with the complex geospatial nature of this phenomenon. We are taking into account all advantages as well as the disadvantages of the used UAV-based solution, whose major problem is the relatively low resolution of the thermal imager. Therefore, currently, we are moving to a more advanced solution based on the integration of thermal radiometric infrared (FLIR) and and a photogrammetric camera. This is the Duet T sensor of the Sensefly Ebee X platform (https://www.sensefly.com/camera/sensefly-duet-t-thermal-mapping-camera/, accessed on 15 September 2021), which has provided very promising results in our preliminary research experiments.

In our opinion, because of the need for data with an extremely high spatial resolution for the changing urban environment, as well as based on the requirements for the time in which the field mapping activities of the thermodynamic characteristics of urban areas should be carried out, modern technological means, such as unmanned aerial systems, provide a high level of flexibility and operability while providing opportunities for accurate collection of high-resolution geospatial data, objectively characterizing various aspects of the territory, including land cover, land use, urban morphology, and local climatic conditions.

The latest technological innovations related to thermal photogrammetry and hyperspectral sensors installed on UAVs open up new possibilities for detailed studies on the urban climate and the factors that determine it. The integration of these new geospatial technologies with traditional geographic information systems and tools will contribute to a better understanding of the mechanism of manifestation and spatial distribution of UHI effects, which is an essential prerequisite for mitigating its adverse effects on urban areas and their adaptation to climate change.