The Urban Heat Island Analysis for the City of Zagreb in the Period 2013–2022 Utilizing Landsat 8 Satellite Imagery

Abstract

Urban Heat Island (UHI) is a phenomenon specific to urban areas where higher air temperatures manifest in the city area in relation to its surrounding rural landscape. Currently, UHI is one of the most dangerous environmental conditions for cities as well as their residents. It is expected that the intensity of UHI will increase with climate change. This work presents an analysis of the UHI phenomenon for the City of Zagreb, Croatia in the summertime period 2013–2022. In order to explore UHI, Land Surface Temperature (LST) was calculated using Landsat 8 (OLI TIRS sensor) satellite imagery. After the delineation of UHI, calculated temperatures were put in relation to NDVI (Normalised Difference Vegetation Index) and NDBI (Normalised Difference Built-Up Index) indices for the study area. Results show the similarity of mean temperatures over the observed period. However, the influence of external variables on UHI’s spatial expression was observed. Forest-covered areas and other green parts of the city’s infrastructure express the lowest temperatures, while built-up sites are the hottest points in cities. Results confirm the importance of urban green infrastructure for resilient cities and present the results of a long-term UHI observation in a Southeast European city.

1. Introduction

Urban Heat Island (UHI) is a phenomenon that occurs when urban areas experience higher temperatures compared to their surrounding non-urban areas [1,2]. UHI is one of the most expressed characteristics of the urban climate caused by human activity in the landscape [3]. The negative impacts of UHI pose significant ecological and social problems in cities, including health risks for citizens due to high temperatures, increased water and energy consumption, and economic loss [4]. Given that more than half of the world’s population today is considered urban and that urbanisation is still ongoing, UHI, with its negative consequences for human health and the physical world, needs to be addressed when planning for the long-term sustainable development of cities. Initially, the UHI phenomenon was expressed in urban–rural linearity, where the difference in temperature between an urban area and a rural area indicated the existence of the long-known and well-studied phenomenon in scientific research [5,6]. However, recently, the research focus shifted from large-scale landscapes to a city- (meso-scale) or micro-scale, emphasising residents’ thermal comfort and the UHI intensity [7]. Because of the mentioned negative effects, UHI is one of the most important topics in research related to urban climate. It is still relevant, especially with climate change and its interaction [8,9]. UHI has been observed and demonstrated in numerous cities worldwide [10,11,12,13,14,15,16]. It is often quantified through LST (Land Surface Temperature), representing one of several possible quantifications of the UHI phenomena [17]. Intensification of the published research exploring UHI within city boundaries can be noticed since 2014 at different spatial scales from meso- to micro-scale [8]. The diverse LST values within urban landscapes are expressed as the SUHI or the surface urban heat island and have far-reaching consequences for the city and its residents [6,15,18].

Even though the factors influencing the emergence of the UHI are complex and diverse, the main one relates to a difference in the land cover between rural and urban areas. There are also differences in their thermal properties or the land cover and its physical properties when it comes to cities. The UHI phenomenon is predominantly caused by human influence on the landscape, i.e., by construction, which creates the difference in the materials covering urban areas and those that cover natural areas [19]. Land cover and land use patterns in the city significantly influence the intensity and the emergence of the UHI [7,14,20]. The second important factor influencing the emergence and the intensity of the UHI relates to favourable meteorological conditions such as clear skies and calm weather [21]. Usually when exploring UHI, two main approaches include measurements of air and land surface temperatures. There are several techniques to study the thermal behaviour of a site, such as remote sensing [6,22], data from fixed meteorological stations, data from qualified entities, in situ campaigns with portable thermal cameras, etc. [8]. Some studies adopt more than one of the aforementioned data sources to complement the information at the atmospheric and surface levels (UHI and SUHI, respectively). This variety of sources allows access to diurnal and/or nocturnal data in different year seasons, although the effect is more intense in summer and winter [11,16,23,24,25,26]. An intra-city variability of the calculated LST is caused by the heterogeneous land cover and land use in a city landscape, which influenced the spatial emergence of the UHI [11]. There is also a time variation of the UHI effect in the city regarding the part of the day when it emerges [14]. Consequently, many studies monitor and model spatial–temporal variation of the phenomenon in a study area [27]. Furthermore, vertical and horizontal distribution and intensity of the UHI are explored within the city with regard to where the UHI occurs, differentiating between surface heat island, canopy layer heat island, and boundary layer heat island, expressing how serious and far-reaching thermal disbalance can occur in cities [8]. Finally, modelling is also one of the possible methods to capture UHI emergence, with developed models being manifold, usually in relation to the main aim of the study [27].

Remote sensing is one of the most often used methods in UHI phenomenon research employed on all spatial scales [8]. Some sensors used for remote sensing can record and measure land surface temperature, and this property enable successful research on UHI around the world [6]. At the same time, Landsat satellite imagery, with its medium spatial resolution, is proven appropriate for the research of particular city areas [28]. Landsat satellite imagery has been successfully used in research on small and large spatial scales for LST measurements and the delineation of UHI [12,29,30]. Landsat 8 has a two-channel infrared TIRS sensor (Thermal Infrared Sensor). TIRS sensor can note neighbouring parts of the spectrum ranging from 10 to 12 μm, which allows the precise calculation of the LST using the data about atmospheric throughput and reflection collected by the satellite. Based on the aforementioned, Landsat 8 allows the simple calculation of the LST, and, consequently, the UHI, environmental critical index (ECI), and other indicators related to the LST from the same satellite imagery [18]. Thermal remote sensing represents one of the essential tools and methods employed that results in high-precision measurements of thermal properties inside cities with a resolution ranging from 30 m to 1 km, depending on a chosen satellite [8]. At the same time, an available spatial resolution could present a limitation to using remote sensing because it does not provide enough detail to study the phenomena in great detail, such as the building–environment relationship [27].

Scientific research on LST and UHI showed that the relationship between them is somewhat complex and that surface temperature response is a function of, among others, different land cover [13,14,31]. The said has led to further research on the relationship between LST and land cover, particularly vegetation abundance [22,32,33,34]. Remote sensing techniques can be used to obtain various vegetation indices and assess vegetation cover [35]. The Normalised Difference Vegetation Index (NDVI) has been widely used to extract vegetation information. Higher NDVI values indicate a higher vegetation coverage in a pixel. The Normalised Difference Built-up Index (NDBI) has been widely used for urban area extraction. Building a connection between land cover and LST can be valuable for urban climate studies [28]. Green infrastructure, primarily urban and peri-urban forests, is considered an effective tool for regulating urban climate and reducing the negative impact of climate change [9,12,20,36]. This property of the green infrastructure is further highlighted in an ecosystem services concept as part of regulating ecosystem services [37]. Urban green infrastructure, with its numerous ecosystem services, is an invaluable part of the city’s infrastructure [38]. The complex mechanism behind the green infrastructure’s cooling effect is evapotranspiration [39] and shading provided by the canopy cover [40].

For the City of Zagreb, the study area in this research, data on the relationship between LST (SUHI), and several vegetation indices have already been presented for 2017. The results showed a significant contribution of vegetation to urban climate regulation in a local context [18]. Furthermore, the results of this study also confirm the suitability of the remote sensing methodology in calculating LST for the City of Zagreb, with high correlation values found between in situ measurements and calculated LST during the day time [18]. In the study area, more research has been conducted in climatology and meteorology using more in situ measurements to describe and explain long-term climate change in the City of Zagreb [41,42,43]. On a European level, it has been calculated that vegetation provides a cooling effect of, on average, 1–1.5 °C, with only a small percentage (less than 10%) of an area showing no cooling effect, for the Functional Urban Area of Zagreb [12]. However, longer-term research of the UHI’s spatial distribution in the City of Zagreb using remote sensing was identified as lacking in the current understanding of Zagreb’s climate. Climate change causing a significant increase in the temperature worldwide is a relevant factor to be considered when planning for a city’s sustainable development in the long term [13]. Understanding current trends and factors affecting UHI can help the City of Zagreb prepare for future environmental challenges. Since satellite imagery is a widely available tool, this paper fills the gap by utilising satellite data and presenting long-term monitoring of the SUHI effect in the City of Zagreb over nine years. This data can be beneficial as a basis for future research in the context of the City of Zagreb and in developing sustainable urban development practices.

In this study, the effects of SUHI and the relationship between LST, NDVI, and NDBI have been studied for the City of Zagreb in the period 2013–2022. The aims of this study are as follows:

• To present long-term monitoring of SUHI manifestation in the City of Zagreb in a summertime period;

• To explore the relationship between calculated LST and NDVI, and NDBI values for the same period;

• To describe the behaviour of SUHI manifestation in the city of Zagreb in relation to calculated indices and detect other possible influences.

2. Materials and Methods

2.1. Study Area

The City of Zagreb is the capital of the Republic of Croatia, located in the continental part of Croatia. It is a mid-sized Southeastern European city with a climate categorised as Cfb, a moderately warm, humid climate with warm summers [44]. Recent research shows changes in the city’s climatic characteristics, with more expressed arid features in some months and extreme temperature regimes at some meteorological stations in the city [41,42]. Two large natural areas in the city that contribute to local climate regulation and serve as a solution to mitigate the negative effects of climate changes are Medvednica Mountain and Sava River. On the one hand, Medvednica is a nature park in the north of the city, covering approximately 8500 ha of the city’s area. On the other, the Sava River runs through the city in the direction of west–east and holds ecological and historical significance. In addition to these two areas, other smaller patches of green spaces, parks, larger urban forests, and water features also contribute to the city’s climate regulation. The study area is presented in Figure 1.

2.2. Satellite Imagery

Six Landsat 8 OLI-TIRS images were obtained from the USGS (United States Geological Survey) website. The images covered the City of Zagreb (path/row 190/28, WGS84, UTM33). The criteria for selecting the images for the analysis were (1) taken in August of a particular year and (2) with minimal cloud coverage. August was chosen as it has the highest number of sunny days in Zagreb, thus increasing the chances of meeting the cloud requirements. Only one image per year was selected to represent the year. This research is exploratory and aims to observe the UHI intensity for a period of years, not conduct a time series study, thus we believe that one image in the hottest period of the year is sufficient to meet the aims of this paper. Only the necessary bands for the analyses were downloaded, hence bands 4, 5, 6, and 10.

2.3. Time Period

The selected period for collecting and analysing satellite imagery was from 2013 to 2022. Landsat 8 imagery was chosen as the main data source, thus the earliest year included in this research was 2013. By examining available data before download, it was found that imagery that met the proposed requirements was available for the years 2013, 2015, 2017, 2019, 2020, and 2022. Therefore, imagery from these years was obtained and used. Imagery from other years did not meet the cloud criteria for August. The dates when the chosen imagery was taken were: 5 August 2013; 11 August 2015; 25 August 2017; 6 August 2019; 8 August 2020; and 6 August 2022.

2.4. Preparation and Calibration of the Raw Temperature Data

Calibration and transformation of the raw data are prepared and performed following the Landsat 8 (L8) Data User Handbook [45]. The software used for the data analysis was QGIS, v 3.6.14. Hannover (QGIS Geographic Information System. QGIS Association. http://www.qgis.org) and R v 4.1. with RStudio (R Core Team (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria).

The complete satellite image was cropped to cover only the City of Zagreb using the city’s boundary vector data to reduce the time needed for the algorithm to run. This step was performed for each image in the analysis. Once the image size was reduced, the calibration was performed following the Handbook. Even though two Landsat 8 bands carried thermal data, band 10 and band 11, due to higher precision, it was advisable to use only band 10 data for calculating LST.

The calibration procedure was as follows:

Digital Number (DN) transformation into Top of Atmosphere emission using an Equation (1). $M_L$ and $Q_{cal}$ values were extracted from the imagery’s metadata (RADIANCE_MULT_BAND_10 and RADIANCE_ADD_BAND_10, respectively) for each image, while $A_L$ is the value of the pixel expressed as DN.

(1)$L_{\lambda }=M_L{Q}_{cal}+A_L$

Furthermore, the Top of Atmosphere number was converted into reflection temperature in Kelvin using Equation (2). $K_1$ and $K_2$ are conversion constants extracted from the metadata file (K1_CONSTANT_BAND_10 and K2_CONSTANT_BAND_10, respectively) for each processed image.

(2)$T_K=\frac{K_2}{\mathrm{l}\mathrm{n} \left(\frac{K_1}{L_{\lambda }}+1\right)}$

The temperature in Kelvin was further converted into temperature in Celsius using Equation (3):

(3)$T_{°\mathrm{C}}=T_K-273.15$

Furthermore, the fraction of vegetation in a pixel was determined using an Equation (4) following the protocol by Guha et al. [46], where NDVI_min value is set to 0.2 and represents the bare soil, and NDVI_max value is set to 0.5, representing pixels considered as healthy vegetation:

(4)$F_{\upsilon =}\left(\frac{NDVI-{NDVI}_{min}}{{NDVI}_{max}-{NDVI}_{min}}\right)$

where NDVI is calculated using a common Equation (5):

(5)$NDVI=\frac{band5-band4}{band5+band4}$

For this purpose, L8 bands 4 and 5 were downloaded, cropped to the area of the City of Zagreb using city boundaries vector data, and used as such.

To calculate LST, it is important to also calculate ground emission beforehand due to the influence of ground characteristics on the surface temperature. Ground emission is calculated using Equation (6).

(6)$\epsilon =0.004 \ast F_{\upsilon }+0.986$

Lastly, LST was calculated using Equation (7).

(7)$LST=\frac{T_{°\mathrm{C}}}{1+\left(\frac{0.00115 \ast T_{°\mathrm{C}}}{1.4388}\right) \ast \mathrm{l}\mathrm{n} \left(\epsilon \right)}$

Urban Heat Island threshold value was calculated following the equation:

(8)$UHI>\mu +\frac{\sigma }{2}$

where μ represents the mean calculated LST value for the study area and σ is the standard deviation of the calculated LST for the study area. Following, UHI in the City of Zagreb was spatially delineated. For each observed year, a map was created highlighting areas in the city expressing UHI and its intensity expressed in °C. Everything above the threshold value calculated for each year can be considered UHI. Spatial delineation of the UHI for the City of Zagreb was then conducted based on individual calculated values for each year.

2.5. Calculating Spatial Indices

Commonly assessed indices were calculated for the City of Zagreb study area. These indices are NDVI and NDBI.

NDVI was calculated for the study area using Equation (5). While NDBI was calculated for the study area using Equation (9).

(9)$NDBI=\frac{band6-band5}{band6+band5}$

For the simplicity of the analysis, LST, NDVI, and NDBI were calculated on a city level. Descriptive statistics for the City of Zagreb were calculated, and the mean values were used to interpret the relationship between LST and spatial indices. Additionally, more detailed descriptive statistics were calculated on a city district level, and the calculated values were plotted. A linear regression line was fitted to explain the trend and relationship between mean calculated LST, NDVI, and NDBI values for the first and the last year explored between city districts in the City of Zagreb. The first and last years were chosen to detect any changes in an urban environment that may influence the UHI’s distribution, intensity, and possible mitigation strategies.

3. Results and Discussion

Cartographic Representation of the Delineated UHI in the City of Zagreb

Figure 2 shows the delineation of the UHI in the City of Zagreb for each year analysed. The blue areas on the maps indicate parts of the City of Zagreb where there is no UHI effect, meaning the temperature is lower than the calculated UHI threshold value. Other parts of the map present those areas in the City of Zagreb that are warmer than the UHI threshold value calculated, hence areas that can be interpreted as SUHI. Thresholds values calculated using Equation (8) were 29.57 °C (2013), 28.80 °C (2015), 30.59 °C (2017), 25.34 °C (2019), 25.17 °C (2020), and 29.73 °C (2022). The visual interpretation of the maps resulted in several conclusions, i.e., UHI in the City of Zagreb is usually connected to built-up areas and remains relatively constant in its spatial distribution over the years. The city’s southwest part, which has a rural character (see Figure 1), confirms this observation as it exhibited lower temperatures than urban areas. Since the UHI definition states that the difference in the temperature between urban and rural landscape indicates the existence of the phenomena, observing differences between the urban and rural part of the city can provide confirmation without looking above the current city borders. The distribution of the UHI also depends on air temperature and local climate expressed in the summertime of each year, as evidenced by the maps with the years 2013, 2017, and 2022 being hotter and having larger UHI in extent and intensity than the years 2019 and 2020 (see Table 1). Specific hotspots, the hottest points in the City of Zagreb, coincide with extensive facilities such as shopping centres and industrial zones that consume and emit a lot of energy and residential areas. The same conclusion was presented in a year-long study of Zagreb’s UHI, confirming the results’ reliability [18]. Industrial zones in the western and eastern parts of the city are continuously represented as UHI hotspots on the maps. These areas are also the most critical environmental points in the city and should be monitored in the future [18]. The results also show the connection between forest areas and the mitigation of the UHI in the city. Indeed, it was confirmed that forests are one of the best solutions available to cities to reduce the UHI effect and other negative consequences of climate change that can harm citizens and reduce their quality of life [47]. Maps also show the influence of the Sava River in the middle of the city on UHI manifestation, where the river stream can be clearly distinguished as a non-UHI area on the maps. However, the relationship between the water features and UHI manifestation is complex, with water sometimes intensifying the UHI effect, especially during the night time [39]. On the other hand, as is the case with this study also, water features can sometimes help mitigate UHI intensity, which is also shown in other research [7]. The visual interpretation also shows that the vegetation around the river expressed UHI characteristics compared to years without the influence of drought and heat waves during heatwave periods (2019 and 2020).

With a more detailed interpretation of the maps, it can be noticed that other hotspot points in the city are usually large concrete or asphalt areas, along with the dump hill that emerged as the hotspot. All of the mentioned highlight the usability of the method, and these results, presented in visually clear form, can easily be communicated to stakeholders and the interested public.

The LST was calculated for the City of Zagreb using remote sensing technology and satellite data. Although there are some differences in temperatures measured on-site and off-site, the advantage of remote sensing and satellite imagery is their ability to cover large areas and produce detailed data at low or no cost. Studies have shown that remote sensing methods are precise enough to be used in research and practice, despite their differences [9]. Furthermore, extracting LST values from remote sensing data has become one of the most crucial components of current urban climate studies that enables long-term spatial assessments of the phenomenon [8].

Descriptive statistics for the LST were calculated for the entire City of Zagreb as a single unit and are presented in Table 1. Mean temperature ranges from 24 °C to 29 °C, and minimum temperatures range from 14 °C to 21 °C while maximal temperatures are measured in the range from 37 °C to 45 °C. The year 2017 was the hottest in the observed period, with a severe heat wave throughout Croatia and Zagreb. The previously said was confirmed by exploring meteorological data for July/August 2017, which showed numerous cities breaking historical temperature records in Croatia, including Zagreb, with the highest calculated UHI threshold value. The cartographic representation of the UHI for the same year shows a large UHI spatial distribution and more severe temperature hotspots. The second hottest year was 2022 when an extreme heat wave and drought affected most of Europe, including Croatia. The combined effect of the heat wave and drought can also be seen in cartographic representation, where a large area indicating UHI can be detected. Alongside built-up areas, bare land and agricultural areas were also part of UHI’s spatial distribution. This suggests that air temperature and other external influences can impact UHI spatial expression in the landscape. The third hottest year was 2013, when a heat wave was also present, according to historical meteorological data. However, by comparing UHI spatial distribution in 2013, 2017, and 2022, it can be observed that UHI was smaller in spatial extent in 2013, likely due to already mentioned external influences that were not controlled in this study but could influence UHI manifestation and magnitudes such as drought [48]. The scientific literature confirms the relationship between heat waves and the amplification of UHI intensity, making them an even more serious threat to the urban population [4,30]. The spatial distribution of the UHI in the City of Zagreb under the heat waves was similar to those in other cities, such as Cluj-Napoca in Romania, where agricultural areas without vegetation during a severe heat wave became part of the detected UHI [4]. The lowest minimum temperature was recorded in 2022. It should be stressed here that this temperature was due to the small number of clouds present in satellite imagery despite efforts to use cloudless imagery. The issue of clouds is well-known in remote sensing, and their presence can sometimes hinder the analysis and interpretation of satellite imagery. Because of that, one of the requirements for including satellite imagery in the analysis was as few clouds as possible to get the most reliable results, even though with these requirements we chose the hottest days. Other years that were included in the analysis, such as 2019 and 2020, showed lower temperatures and smaller UHI spatial distribution throughout the City of Zagreb for a similar period.

Compared to the temperatures measured using satellite sensors presented in this paper, meteorological tools at official measurement sites in Zagreb showed lower maximal temperatures for the City of Zagreb than those presented here for daytime temperatures [41]. However, the advantage of the presented approach is its spatial extent and pixel-based measurements that provide high precision in detecting spatially explicit hotspots in the city. For the City of Zagreb, these micro-locations, such as industrial sites and shopping facilities, increase the mean temperature and hold the highest LST interpreted here as the maximal temperature.

NDVI is a widely used spatial index, which is relatively easy to calculate and interpret. NDVI can be calculated from satellite imagery, like LST, and compared to LST using the same dataset (imagery). The descriptive statistics of the NDVI for the City of Zagreb are presented in Table 2. It can be seen in the table that the City of Zagreb, overall, has a satisfactory level of vegetation, with mean values being around 0.37, which is relatively high compared to the mean NDVI values in cities like Naples and Florence in Italy [46]. The difference between the cities emerges from different local climates and the particular position of the City of Zagreb, which is on the slopes of the Medvednica mountain. The Medvednica mountain is a large, forested area in the north of the city. Therefore, its significance is acknowledged on the city level and beyond it and is essential not only as part of the green infrastructure of the City of Zagreb but also as a provider of different ecosystem services, including climate regulation. Alongside Medvednica, patches of these forests remained in the urban landscape and formed several urban forests in the City of Zagreb. Finally, private- and state-owned forests in the city contribute to a high NDVI along with a typical urban green infrastructure such as parks, tree alleys, and others. The contribution of trees and forests in cities is present through the provision of many ecosystem services and a significant cooling effect, which is the most expressed in the summertime [9]. Minimal and maximal values are similar to those in Naples and Florence [46].

The second commonly used spatial index that can be derived from remote sensing imagery is the NDBI index. The built-up counterpart to NDVI indicates how built-up the observed area is. Table 3 presents the descriptive statistics of the NDBI values for the City of Zagreb for each year. The mean values for the City of Zagreb are higher than the reported values for the cities of Naples and Florence [46] but similar to those reported for the Lagos metropolis in Nigeria [49].

Finally, the mean LST was calculated for each city district in the City of Zagreb using OpenStreetMap’s city district boundaries vector data. The LST and mean NDVI and NDBI were also extracted for each district. This was done for 2013 and 2022 as the first and last years in the selected time frame. Two plots were produced and presented in Figure 3, one exploring the relationship between mean temperature in each city district and its mean NDVI, and the second one exploring the relationship between mean temperature and mean NDBI in each city district in the City of Zagreb. Consistent with previous research, we found a strong negative relationship between LST and NDVI, meaning that those city districts with a high share of vegetation had a lower LST. This trend was observed both in 2013 and 2022, indicating that this association is universal. In a large-scale European study, Zagreb was found to have a cooling effect from the vegetation of 1–1.5 °C, placing it among cities with sufficient greenery to provide benefits to residents [12]. Here, it should be stressed that one-third of the City of Zagreb is allocated to forests, one-third to agricultural areas, and one-third to built-up areas. Thus, forests significantly impact the city’s cooling capacity, with larger forested areas having a greater extent of cooling than smaller green spaces, which is confirmed in other studies [20].

Contrary to the relationship between NDVI and LST, NDBI and LST showed an opposite trend. Higher NDBI values were found to be related to higher LST. This is due to the nature of UHI, where impervious surfaces emit energy and warm up the surrounding air. Several city districts in the City of Zagreb have high LST due to being heavily built up. Consequently, these city districts are also those with low NDVI values. This result confirms the usefulness of the statistical approach applied. The use of remote sensing provides a detailed view of UHI in the city and identifies those areas that require special attention. Furthermore, regression analysis, a commonly used statistical approach in many studies, can help identify those factors and indices that influence the distribution and intensity of UHI, thus improving our understanding of the UHI phenomenon [8].

Regression analysis showed interesting trends in the City of Zagreb, particularly regarding NDBI. The results indicate that NDBI values increased in 2022, suggesting an increase in new constructions throughout the city. On the other hand, NDVI values decreased in 2022 compared to 2013. This trend can be attributed to converting agricultural land into built-up areas, as seen in a comparison of land use/land cover categories in Zagreb between 1968 and 2012 [43]. The same research shows that heat load can be observed even in the forest land cover in a long-term period, but to a lesser extent than other land covers [43]. Since the Republic of Croatia was a transition economy, the trend of new construction was expected, given that the City of Zagreb is the capital of the Republic of Croatia and an attractive location for investment. Long-term planning and management should prioritise sustainable development and citizens’ well-being, highlighting urban green infrastructure and its ecosystem services as a potential strategy for mitigating the UHI effect.

4. Limitations

Some limitations in this study can be addressed. Firstly, this study only derives temperature from satellite imagery and not in situ measurements using meteorological tools. Research suggests combining both methods provides a more comprehensive understanding of the UHI phenomena [50]. Secondly, this study only uses one clear image from August for each observed year, making it difficult to draw general conclusions about the LST in the summertime in Zagreb. However, differences in the LST temperatures and spatial distributions can serve as a foundation for future, more detailed research on UHI in the City of Zagreb, particularly during the night time as the imagery used in this research was captured during the daytime (morning) hours, due to satellite trajectory.

5. Conclusions

This paper presents an assessment of the LST for the City of Zagreb over a longer term and the spatial distribution of UHI. The study also explores the relationship between LST and commonly used spatial indices (NDVI, NDBI). The results are shown as a cartographic representation of UHI, which highlights the parts of the city with higher temperatures in the summertime and where they can be harmful to the residents or visitors, such as shopping centres and open, large, concrete areas. Cartographic representation helps inform stakeholders and the public of the potential risks associated with prolonged exposure to UHI. It also serves as a foundation for future research and city planning focused on combating the negative impacts of climate change.

A negative relationship between NDVI and calculated LST in city districts with a large proportion of green areas indicates the critical role of forests and other green areas in providing thermal comfort in the city. Conversely, a positive relationship was observed between higher NDBI values and LST, highlighting UHI and the relationship between highly built-up areas and thermal discomfort for city residents. This study also observed other factors influencing UHI, such as long-term drought. This opens up new scientific questions and avenues for further research. This paper presents one of the first assessments of UHI intensity and spatial distribution over a long period in the City of Zagreb, the capital and largest city in Croatia. The southern part of Europe, including Croatia, is recognised as a concerning area regarding heat waves and climate change, therefore, it is necessary to establish long-term monitoring of UHI to develop successful mitigation strategies. The findings of this study could provide a foundation for future research. The year 2022 is emphasised as concerning since severe drought and a heat wave strongly influenced the largest UHI delineated in the City of Zagreb, expressing the adverse effects of climate change on the environment and citizens alike. Future research could explore the effects of heatwaves on night time temperatures, their impact on citizens’ health, and the role of vegetation in mitigating UHI in cities under heatwave pressure, among other topics.

Author Contributions

Conceptualization: A.S. and M.K.; Data curation: A.S. and M.K.; Formal analysis: A.S. and M.K.; Visualization: M.K.; Supervision: A.S. and R.P.; Writing—original draft: A.S. and M.K.; Writing—review and editing: M.A., J.K. and R.P. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data will be provided upon request.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. The study area. CLC (Corine Land Cover) for the city of Zagreb with city district boundaries (OpenStreetMap data) (a); and the study area in a spatial context of Europe (b).

Figure 2. Spatial distribution of the UHI in the city of Zagreb in time period 2013–2022.

Figure 3. Relationship between LST (°C) and NDVI (a) and NDBI (b) index calculated for each city district in the City of Zagreb with regression line for the year 2013 (blue) and 2022 (red).

Table 1

Descriptive statistics of the LST (°C) for the city of Zagreb in each year.

Year	MIN	MAX	MEAN	SD
2013	21.27	41.95	28.19	2.77
2015	18.72	42.56	27.50	2.62
2017	20.73	44.94	29.31	2.56
2019	17.42	37.38	24.20	2.28
2020	18.11	38.38	24.11	2.12
2022	13.97	38.58	28.46	2.53

Table 2

Descriptive statistics of NDVI values for the city of Zagreb in each year.

Year	MIN	MAX	MEAN	SD
2013	−0.11	0.61	0.36	0.11
2015	−0.13	0.61	0.38	0.11
2017	−0.15	0.60	0.36	0.12
2019	−0.15	0.63	0.39	0.12
2020	−0.14	0.63	0.40	0.11
2022	−0.15	0.62	0.34	0.12

Table 3

Descriptive statistics of NDBI values for the city of Zagreb in each year.

Year	MIN	MAX	MEAN	SD
2013	−0.54	0.63	−0.16	0.1
2015	−0.40	0.61	−0.18	0.09
2017	−0.44	0.48	−0.15	0.1
2019	−0.44	0.42	−0.19	0.1
2020	−0.21	0.59	−0.21	0.09
2022	−0.52	0.39	−0.13	0.11

References

1. Rizwan, A.M.; Dennis, L.Y.C.; Liu, C. A Review on the Generation, Determination and Mitigation of Urban Heat Island. J. Environ. Sci.; 2008; 20, pp. 120-128. [DOI: https://dx.doi.org/10.1016/S1001-0742(08)60019-4]

2. Oke, T.R.; Mills, G.; Christen, A.; Voogt, J.A. Urban Climates; Cambridge University Press: Cambridge, UK, 2017.

3. Yang, L.; Qian, F.; Song, D.-X.; Zheng, K.-J. Research on Urban Heat-Island Effect. Procedia Eng.; 2016; 169, pp. 11-18. [DOI: https://dx.doi.org/10.1016/j.proeng.2016.10.002]

4. Herbel, I.; Croitoru, A.-E.; Rus, A.V.; Roşca, C.F.; Harpa, G.V.; Ciupertea, A.-F.; Rus, I. The Impact of Heat Waves on Surface Urban Heat Island and Local Economy in Cluj-Napoca City, Romania. Theor. Appl. Climatol.; 2018; 133, pp. 681-695. [DOI: https://dx.doi.org/10.1007/s00704-017-2196-4]

5. Ünal, Y.S.; Sonuç, C.Y.; Incecik, S.; Topcu, H.S.; Diren-Üstün, D.H.; Temizöz, H.P. Investigating Urban Heat Island Intensity in Istanbul. Theor. Appl. Climatol.; 2020; 139, pp. 175-190. [DOI: https://dx.doi.org/10.1007/s00704-019-02953-2]

6. Voogt, J.A.; Oke, T.R. Thermal Remote Sensing of Urban Climates. Remote Sens. Environ.; 2003; 86, pp. 370-384. [DOI: https://dx.doi.org/10.1016/S0034-4257(03)00079-8]

7. Naikoo, M.W.; Islam, A.R.M.T.; Mallick, J.; Rahman, A. Land Use/Land Cover Change and Its Impact on Surface Urban Heat Island and Urban Thermal Comfort in a Metropolitan City. Urban Clim.; 2022; 41, 101052. [DOI: https://dx.doi.org/10.1016/j.uclim.2021.101052]

8. Kim, S.W.; Brown, R.D. Urban Heat Island (UHI) Variations within a City Boundary: A Systematic Literature Review. Renew. Sustain. Energy Rev.; 2021; 148, 111256. [DOI: https://dx.doi.org/10.1016/j.rser.2021.111256]

9. Lafortezza, R.; Chen, J. The Provision of Ecosystem Services in Response to Global Change: Evidences and Applications. Environ. Res.; 2016; 147, pp. 576-579. [DOI: https://dx.doi.org/10.1016/j.envres.2016.02.018]

10. Santamouris, M. Analyzing the Heat Island Magnitude and Characteristics in One Hundred Asian and Australian Cities and Regions. Sci. Total Environ.; 2015; 512–513, pp. 582-598. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2015.01.060]

11. Chakraborty, T.; Hsu, A.; Manya, D.; Sheriff, G. A Spatially Explicit Surface Urban Heat Island Database for the United States: Characterization, Uncertainties, and Possible Applications. ISPRS J. Photogramm. Remote Sens.; 2020; 168, pp. 74-88. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2020.07.021]

12. Marando, F.; Heris, M.P.; Zulian, G.; Udías, A.; Mentaschi, L.; Chrysoulakis, N.; Parastatidis, D.; Maes, J. Urban Heat Island Mitigation by Green Infrastructure in European Functional Urban Areas. Sustain. Cities Soc.; 2022; 77, 103564. [DOI: https://dx.doi.org/10.1016/j.scs.2021.103564]

13. Schwarz, N.; Manceur, A.M. Analyzing the Influence of Urban Forms on Surface Urban Heat Islands in Europe. J. Urban Plan. Dev.; 2015; 141, A4014003. [DOI: https://dx.doi.org/10.1061/(ASCE)UP.1943-5444.0000263]

14. Dugord, P.-A.; Lauf, S.; Schuster, C.; Kleinschmit, B. Land Use Patterns, Temperature Distribution, and Potential Heat Stress Risk—The Case Study Berlin, Germany. Comput. Environ. Urban Syst.; 2014; 48, pp. 86-98. [DOI: https://dx.doi.org/10.1016/j.compenvurbsys.2014.07.005]

15. Geletič, J.; Lehnert, M.; Savić, S.; Milošević, D. Inter-/Intra-Zonal Seasonal Variability of the Surface Urban Heat Island Based on Local Climate Zones in Three Central European Cities. Build. Environ.; 2019; 156, pp. 21-32. [DOI: https://dx.doi.org/10.1016/j.buildenv.2019.04.011]

16. Zhou, B.; Rybski, D.; Kropp, J.P. On the Statistics of Urban Heat Island Intensity. Geophys. Res. Lett.; 2013; 40, pp. 5486-5491. [DOI: https://dx.doi.org/10.1002/2013GL057320]

17. Marando, F.; Salvatori, E.; Sebastiani, A.; Fusaro, L.; Manes, F. Regulating Ecosystem Services and Green Infrastructure: Assessment of Urban Heat Island Effect Mitigation in the Municipality of Rome, Italy. Ecol. Model.; 2019; 392, pp. 92-102. [DOI: https://dx.doi.org/10.1016/j.ecolmodel.2018.11.011]

18. Krtalić, A.; Divjak, A.K.; Čmrlec, K. Satellite-driven assessment of surface urban heat islands in the city of Zagreb, Croatia. ISPRS Ann. Photogramm. Remote Sens. Spat. Inf. Sci.; 2020; 3, pp. 757-764. [DOI: https://dx.doi.org/10.5194/isprs-annals-V-3-2020-757-2020]

19. Žgela, M. Urbana Klimatologija—Primjer Toplinskog Otoka Grada Zagreba. Geogr. Horiz.; 2018; 64, pp. 31-40.

20. Jaganmohan, M.; Knapp, S.; Buchmann, C.M.; Schwarz, N. The Bigger, the Better? The Influence of Urban Green Space Design on Cooling Effects for Residential Areas. J. Environ. Qual.; 2016; 45, pp. 134-145. [DOI: https://dx.doi.org/10.2134/jeq2015.01.0062]

21. László, E.; Bottyán, Z.; Szegedi, S. Long-Term Changes of Meteorological Conditions of Urban Heat Island Development in the Region of Debrecen, Hungary. Theor. Appl. Climatol.; 2016; 124, pp. 365-373. [DOI: https://dx.doi.org/10.1007/s00704-015-1427-9]

22. Weng, Q.; Lu, D.; Schubring, J. Estimation of Land Surface Temperature–Vegetation Abundance Relationship for Urban Heat Island Studies. Remote Sens. Environ.; 2004; 89, pp. 467-483. [DOI: https://dx.doi.org/10.1016/j.rse.2003.11.005]

23. Santamouris, M.; Cartalis, C.; Synnefa, A. Local Urban Warming, Possible Impacts and a Resilience Plan to Climate Change for the Historical Center of Athens, Greece. Sustain. Cities Soc.; 2015; 19, pp. 281-291. [DOI: https://dx.doi.org/10.1016/j.scs.2015.02.001]

24. Santamouris, M. On the Energy Impact of Urban Heat Island and Global Warming on Buildings. Energy Build.; 2014; 82, pp. 100-113. [DOI: https://dx.doi.org/10.1016/j.enbuild.2014.07.022]

25. Tan, J.; Zheng, Y.; Tang, X.; Guo, C.; Li, L.; Song, G.; Zhen, X.; Yuan, D.; Kalkstein, A.J.; Li, F. et al. The Urban Heat Island and Its Impact on Heat Waves and Human Health in Shanghai. Int. J. Biometeorol.; 2010; 54, pp. 75-84. [DOI: https://dx.doi.org/10.1007/s00484-009-0256-x]

26. Keikhosravi, Q. The Effect of Heat Waves on the Intensification of the Heat Island of Iran’s Metropolises (Tehran, Mashhad, Tabriz, Ahvaz). Urban Clim.; 2019; 28, 100453. [DOI: https://dx.doi.org/10.1016/j.uclim.2019.100453]

27. Mirzaei, P.A. Recent Challenges in Modeling of Urban Heat Island. Sustain. Cities Soc.; 2015; 19, pp. 200-206. [DOI: https://dx.doi.org/10.1016/j.scs.2015.04.001]

28. Kaplan, G.; Avdan, U.; Avdan, Z.Y. Urban Heat Island Analysis Using the Landsat 8 Satellite Data: A Case Study in Skopje, Macedonia. Proceedings; 2018; 2, 358. [DOI: https://dx.doi.org/10.3390/ecrs-2-05171]

29. Schmidt, P.; Lawrence, B.T. Association between Land Surface Temperature and Green Volume in Bochum, Germany. Sustainability; 2022; 14, 14642. [DOI: https://dx.doi.org/10.3390/su142114642]

30. Founda, D.; Santamouris, M. Synergies between Urban Heat Island and Heat Waves in Athens (Greece), during an Extremely Hot Summer (2012). Sci. Rep.; 2017; 7, 10973. [DOI: https://dx.doi.org/10.1038/s41598-017-11407-6]

31. Owen, T.W.; Carlson, T.N.; Gillies, R.R. An Assessment of Satellite Remotely-Sensed Land Cover Parameters in Quantitatively Describing the Climatic Effect of Urbanization. Int. J. Remote Sens.; 1998; 19, pp. 1663-1681. [DOI: https://dx.doi.org/10.1080/014311698215171]

32. Chen, X.-L.; Zhao, H.-M.; Li, P.-X.; Yin, Z.-Y. Remote Sensing Image-Based Analysis of the Relationship between Urban Heat Island and Land Use/Cover Changes. Remote Sens. Environ.; 2006; 104, pp. 133-146. [DOI: https://dx.doi.org/10.1016/j.rse.2005.11.016]

33. Gallo, K.P.; Owen, T.W. Satellite-Based Adjustments for the Urban Heat Island Temperature Bias. J. Appl. Meteorol.; 1999; 38, pp. 806-813. [DOI: https://dx.doi.org/10.1175/1520-0450(1999)038<0806:SBAFTU>2.0.CO;2]

34. Weng, Q. A Remote Sensing?GIS Evaluation of Urban Expansion and Its Impact on Surface Temperature in the Zhujiang Delta, China. Int. J. Remote Sens.; 2001; 22, pp. 1999-2014. [DOI: https://dx.doi.org/10.1080/713860788]

35. Florim, I.; Albert, B.; Shpejtim, B. Measuring Uhi Using Landsat 8 Oli and Tirs Data with Ndvi and Ndbi in Municipality of Prishtina. Disaster Adv.; 2021; 14, pp. 25-36. [DOI: https://dx.doi.org/10.25303/1411da2536]

36. Gill, S.E.; Handley, J.F.; Ennos, A.R.; Pauleit, S. Adapting Cities for Climate Change: The Role of the Green Infrastructure. Built Environ.; 2007; 33, pp. 115-133. [DOI: https://dx.doi.org/10.2148/benv.33.1.115]

37. Arico, S.; Bridgewater, P.; El-beltagy, A.; Finlayson, M.; Harms, E.; Program, S.; Hepworth, R.; Leitner, K.; Oteng-yeboah, A.; Ramos, M.A. et al. Millennium Ecosystem Assessment, 2005—Ecosystems and Human Well-Being: Synthesis; Island Press: Washington, DC, USA, 2005; ISBN 1-59726-040-1

38. Livesley, S.J.; McPherson, E.G.; Calfapietra, C. The Urban Forest and Ecosystem Services: Impacts on Urban Water, Heat, and Pollution Cycles at the Tree, Street, and City Scale. J. Environ. Qual.; 2016; 45, pp. 119-124. [DOI: https://dx.doi.org/10.2134/jeq2015.11.0567]

39. Gunawardena, K.R.; Wells, M.J.; Kershaw, T. Utilising Green and Bluespace to Mitigate Urban Heat Island Intensity. Sci. Total Environ.; 2017; 584–585, pp. 1040-1055. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2017.01.158]

40. Oke, T.R. The Micrometeorology of the Urban Forest. Philos. Trans. R. Soc. Lond. B; 1989; 324, pp. 335-349. [DOI: https://dx.doi.org/10.1098/rstb.1989.0051]

41. Ugarković, D.; Matijević, M.; Tikvić, I.; Popić, K. Neka Obilježja Klime i Klimatskih Elemenata Na Području Grada Zagreba. Šumarski List; 2021; 145, pp. 479-488. [DOI: https://dx.doi.org/10.31298/sl.145.9-10.6]

42. Nimac, I.; Herceg-Bulić, I.; Cindrić Kalin, K.; Perčec Tadić, M. Changes in Extreme Air Temperatures in the Mid-Sized European City Situated on Southern Base of a Mountain (Zagreb, Croatia). Theor. Appl. Climatol.; 2021; 146, pp. 429-441. [DOI: https://dx.doi.org/10.1007/s00704-021-03689-8]

43. Nimac, I.; Herceg-Bulić, I.; Žuvela-Aloise, M. The Contribution of Urbanisation and Climate Conditions to Increased Urban Heat Load in Zagreb (Croatia) since the 1960s. Urban Clim.; 2022; 46, 101343. [DOI: https://dx.doi.org/10.1016/j.uclim.2022.101343]

44. Šegota, T.; Filipčić, A. Köppenova Podjela Klima i Hrvatsko Nazivlje. Geoadria; 2003; 8, pp. 17-37. [DOI: https://dx.doi.org/10.15291/geoadria.93]

45. USGS. Landsat 8 (L8) Data Users Handbook; Earth Resources Observation and Science (EROS) Center: Garretson, SD, USA, 2019.

46. Guha, S.; Govil, H.; Dey, A.; Gill, N. Analytical Study of Land Surface Temperature with NDVI and NDBI Using Landsat 8 OLI and TIRS Data in Florence and Naples City, Italy. Eur. J. Remote Sens.; 2018; 51, pp. 667-678. [DOI: https://dx.doi.org/10.1080/22797254.2018.1474494]

47. Nastran, M.; Kobal, M.; Eler, K. Urban Heat Islands in Relation to Green Land Use in European Cities. Urban For. Urban Green.; 2019; 37, pp. 33-41. [DOI: https://dx.doi.org/10.1016/j.ufug.2018.01.008]

48. Kong, J.; Zhao, Y.; Carmeliet, J.; Lei, C. Urban Heat Island and Its Interaction with Heatwaves: A Review of Studies on Mesoscale. Sustainability; 2021; 13, 10923. [DOI: https://dx.doi.org/10.3390/su131910923]

49. Alademomi, A.S.; Okolie, C.J.; Daramola, O.E.; Akinnusi, S.A.; Adediran, E.; Olanrewaju, H.O.; Alabi, A.O.; Salami, T.J.; Odumosu, J. The Interrelationship between LST, NDVI, NDBI, and Land Cover Change in a Section of Lagos Metropolis, Nigeria. Appl. Geomat.; 2022; 14, pp. 299-314. [DOI: https://dx.doi.org/10.1007/s12518-022-00434-2]

50. Elmarakby, E.; Khalifa, M.; Elshater, A.; Afifi, S. Tailored Methods for Mapping Urban Heat Islands in Greater Cairo Region. Ain Shams Eng. J.; 2022; 13, 101545. [DOI: https://dx.doi.org/10.1016/j.asej.2021.06.030]