SATELLITE-DRIVEN ASSESSMENT OF SURFACE URBAN HEAT ISLANDS IN THE CITY OF ZAGREB, CROATIA

This study aims to assess surface urban heat islands (SUHIs) pattern over the city of Zagreb, Croatia, based on satellite (optical and thermal) remote sensing data. The spatio-temporal identification of SUHIs is analysed using the 12 sets of Landsat 8 imagery acquired during 2017 (in each month of the year). Vegetation cover within the city boundaries is extracted by using Principal Component Analysis (PCA) data fusion method on calculated three vegetation indices (VI): Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI) and Ratio Vegetation Index (RVI) for each set of bands. The first principal component was used to compute the land surface temperature (LST) and deductive Environmental Criticality Index (ECI). As expected, the relationship between LST and all VI scores shows a negative correlation and is most negative with RVI. The environmentally critical areas and the patterns of seasonal variations of the SUHIs in the city of Zagreb were identified based on the LST, ECI and vegetation cover. The city centre, an industrial area in the eastern part and an area with shopping centers and commercial buildings in the western part of the city were identified as the most critical areas. * Corresponding author


INTRODUCTION
The surface urban heat island (SUHI) is a phenomenon in urban environmental studies which represents the difference of land surface temperature (LST) between urban and neighbouring non-urban surfaces (Zhou et al. 2019). SUHI occurs in a different way from urban heat island (UHI) (Roth et al., 1989;Tran et al., 2006). UHI is mainly caused by the differences in radiative cooling between urban and rural areas during nighttime, while SUHI is mainly caused by the differences in radiative surface heating between urban and rural areas during the daytime (Choi et al., 2014). Furthermore, according to (Zhou et al. 2019), UHIs can be classified into two broad categories based on the way they are formed and on heights (Oke, 1982): "air" (or "atmospheric") and "surface" UHIs. The UHI is usually measured by in situ observed air temperature data from meteorological stations on the ground (Nichol et al., 2009) which can provide a high temporal frequency and accuracy. But it has limitations in regard to the spatial resolution that can be derived from the observing stations. Due to these limitations, the in situ measured air temperatures usually can"t provide enough spatial details for detection and extraction of UHI. This fact affects the quality of the estimation of spatio-temporal variability in UHIs and the spatial distribution of it. In contrast, the SUHI is primarily measured by satellite thermal remote sensing data, which offers the ability to study the urban thermal environment at various spatial (from local to global scales) and temporal (daily, seasonally, and interannually) scales (Pal et al., 2012, Zhou et al. 2019. With the increase in population and economic factors in cities, concrete areas increase and areas under urban vegetation decrease. These result in increased heat in cities. As a result, densely populated parts of the urban area become warmer than sparsely populated areas and areas with much vegetation. LST increases in the urban area compared to the surrounding rural landscape and forms the urban heat island effect (Oke, 1982). Authors Weng (2001), Badarinath et al. (2005), Mallick et al. (2008) have shown a relationship between LST, built-up, vegetation areas and UHI impact.
Multi-temporal Landsat data have proven to be useful in the study of LST variations. The thermal band of Landsat imagery helps in the calculation of LST on the region of interest (Weng, 2001, Chen et al., 2002. The high-LST regions correspond with high rising residential buildings and industrial areas with low vegetation coverage (Li et al. 2013). The vegetation cover inside the city boundaries affects the physical environment of the city. It helps in the selective absorption and reflection of incident radiation and also regulates the exchange of permanent and sensible heat (Gallo et al. 1995, Nichol 1996. Urban areas covered with vegetation reduce the possibility of creating SUHIs, while the removal of existing vegetation areas contributes to the formation of SUHIs (Foley et al., 2005;Wong and Yu, 2005). Consequently, the availability of vegetation cover is considered as an indicator of ecological sustainability in an urban community (Senanayake et al. 2013a). In this way, the maintenance and increase of vegetation surfaces in urban areas has become one of the most important activities in the effort of sustainable urban development (Bernatzky, 1982, Senanayake et al., 2013b. Temporal analysis of the relationship between LST and Normalised Differential Vegetation Index (NDVI) shows that the relationship amongst these varies temporally and NDVI arises as reliable tool for quantitative analysis of LST over different time in urbanised areas (Chen et al., 2006). The negative differences of NDVI between urban and rural was analysed by Choi et al. (2014). They indicated that the fractional coverage and conditions of vegetation at the rural points were greater than those at the urban points during the daytime. In addition, Grover and Singh (2015) analysed the UHI in relation to NDVI. They used all the pixels of LST and NDVI in the regression analysis, where NDVI was considered as the independent and LST as the dependent variable, as LST has been found to be strongly determined by vegetation health (Yue et al., 2007, Zhang andChen, 2010).
Due to the importance of vegetation cover in LST and derivation of Environmental Criticality Index (ECI) (Senanayake et al. 2013a), two more vegetation indices: Enhanced Vegetation Index (EVI), Ratio Vegetation Index (RVI), beside NDVI, are used to better estimate the emisivity, LST and ECI. EVI was used due to the fact that it was developed to improve the NDVI, and RVI is widely used for vegetation monitoring, specifically, at high density vegetation coverage. For better insight into the surface and vegetation conditions, data fusion was performed using the Principal Component Analysis (PCA) method with all three vegetation indices. Different characteristics of vegetation indices were extracted in this way. This study aims to analyse the association between built area, green cover and LST and ECI for the purpose of detecting and extracting SUHIs over the area of the City of Zagreb, Croatia.

Study Area: City of Zagreb, Croatia
Zagreb (45° 49' N, 15° 59' E) is a medium sized European city and capital of the of Republic of Croatia, located in the central part of the state on the plain bordered by the River Sava to the south and the Mt. Medvednica to the north (Figure 1). The urban area extends along the plain to Medvenica in the south and its slopes in the north. The gently sloping hills of the Medvednica Mountain (north of the Sava River) and the flat land between Medvednica and the Sava River make most of Zagreb. The climate of Zagreb is classified as an oceanic climate, but with significant continental influences and very closely bordering on a humid continental climate as well as a humid subtropical climate, with distinctive four seasons. Zagreb summers are usually warm and wet, but the temperatures can touch 32.2°C (90°F) during July and August at the peak of the season, particularly during a heat wave. Winters are cold, and the temperatures regularly drop below freezing with a moderate amount of snow and icy precipitation. Spring is mild and pleasant with unpredictable weather changes, which can push down the temperatures due to the prominent southerly winds. Autumn is temperate but rainy and foggy during the latter half of the season (Weather Atlas, 2020). A densely populated urban area with 792,875 citizens (according to the last population census in 2011), covering an area of 641.36 m², Zagreb is characterised by moderate urban expansion. The extent of the study area covers the entire municipality of Zagreb, while satellite images were processed within the administrative boundary of the City of Zagreb.

Satellite and Meteorological Data
USGS/NASA Landsat 8 satellite data (OLI and TIRS) covering Zagreb are used in this study for the calculation of LST during each season in 2017 (in total of 10 months). The data included Landsat 8 images obtained on January 20, February 14, March 25, April 3, April 10, May 28, June 22, July 31, August 25, October 12, October 19 and December 6 ( Table 1).
Since there were no usable visible and near infrared Landsat 8 images for the months of September and November (Zagreb area was covered with clouds), LST estimation was not made for these months. Bands 1-7 and 9 (with spectral range from 0.433 to 1.390 µm) are acquired in 30 m resolution while thermal bands, 10 and 11 (with spectral range from 10.6 to 12.5 µm), are acquired in 100 m resolution. Band 8 is the panchromatic band (with spectral range from 0.500 to 0.680 µm) which can get more light at once and it"s the sharpest of all   (Table   2). CMHS data were used to compare with the LST results, that is, to analyse the approximate accuracy of the calculated LSTs.
Grič meteorological station is located near the city center (urban environment), Puntjarka station on the northern slope of Medvednica (rural area) that is mostly covered by forest, and Maksimir station is located within a city park surrounded by buildings (mixed environment) (Figure 1). Land temperatures for the UTC time of 09:40 h when all the satellite images used in this study were collected, were linearly modelled relatively to the in situ measured temperatures at Maksimir station at 07:00 and 14:00 h, for the dates indicated.

Vegetation Indices and Principal Component Analysis
Because vegetation surfaces play an important role in minimising the possibilities for the creation of SUHI and calculation of LST, special attention was given to the calculation of the vegetation index (VI). A similar study was conducted and presented in Kayet et al. (2016), and a methodology has been proposed for the relationships between multiple vegetation indices and LST based on thermal remote sensing data. The author's conclusion is that the experimentation in this field needs to increase. In this study, we used multiple vegetation indices in the process of computing LST. Normalized Difference Vegetation Index (NDVI), the Enhanced Vegetation Index (EVI) and the Ratio Vegetation Index (RVI) are used to enhance the vegetation present in the observed area. NDVI (Kriegler et al., 1969) is the common vegetation index referring to VI. The ratio of the NIR and red band is used for the calculation because the absorption by chlorophyll of these two bands of the electromagnetic spectrum is the highest (Bindi et al. 2009). NDVI was calculated using the (Landsat 8) red (B4) and near-infrared (B5) bands with the Equation (1): EVI (Huete et al., 1997) was developed to improve the NDVI by optimizing the vegetation signal in Leaf Area Index (LAI) regions by using the blue reflectance to correct for soil background signals and reduce atmospheric influences. This VI is therefore most useful in LAI regions, where the NDVI may saturate. EVI was calculated using the (Landsat 8) blue (B2), red (B4) and near-infrared (B5) bands with the Equation (2): A common practice in remote sensing is the use of band ratios to eliminate various albedo effects. In case of using RVI (Jordan, 1969), the vegetation isoline converges at origin. And it is calculated with the Equation (3): The fusion of the NDVI, EVI and RVI was performed using the principal components analysis (PCA) to better evaluate and define vegetation cover. In this case, the input for LST calculation is the first principal component calculated from NDVI, EVI and RVI (PC1_NDVI-EVI-RVI = PC1_NER).

Land Surface Temperature Calculation from Landsat 8 Thermal Bands
Averaged values of both thermal bands (10 and 11) of the Landsat 8 images were used to calculate the LST distribution pattern of the city of Zagreb using remote sensing and image processing techniques (Figure 3). The spectral range of these thermal bands is from 10.60 to 11.12 µm and 10.50 to 12.51 µm respectively. Landsat 8 Data Users Handbook describes the retrieval method of LST from the thermal band of an image (USGS, 2019). In terms of atmospheric correction, the digital number (DN) values of all 12 visible and thermal bands were converted to spectral radiance values by using offset (bias) and gain values of the images as shown in the Equation (4)  where T S = the emissivity corrected LST in Kelvin (K), λ = the wavelength of emitted radiance, q = hc K -1 1 (1.438 * 10 -2 mK, h = Planck"s Constant (6.626 * 10 -34 J s -1 ), c = the velocity of light (2.998 * 10 8 m s -1 ), K = Boltzman constant (1.38 * 10 -23 J K -1 ), ε = surface emissivity.
The calculated temperature values in Kelvin were converted to Celsius degrees by using the Equation (7):

Emissivity Correction
In the Equation (6), only emissivity is an unknown value and needs to be predicted. It is mentioned in (Barane and Dwarakish, 2017) that the method based on LU/LC classified map is the simplest one, but the accuracy of LU/LC classification has significant influence on emissivity prediction, especially when using Landsat 8 satellite images with a spatial resolution of 30 m. In this research, the method which uses ratio values of vegetation and bare land was conducted. It is an easy to use method to predict emissivity which using NDVI image (in this case, the first principal component calculated from NDVI, EVI and RVI: PC1_NER) as given by Sobrino et al. (2004) with the Equation (8): where P v = proportion of vegetation obtained (Carlson et al., 1997), Equation (9).

Environmental Criticality Index Retrieval
In an analysis conducted in (Senanayake et al. 2013a), the authors concluded that LST and NDVI have direct and inversely proportional relationships, respectively, with the environmental criticality level of the analysed area. Based on this fact, they defined a deductive Environmental Criticality Index (ECI) as shown in the Equation (10)  The pixel values of the raster used for ECI computing are stretched from 1 to 255 to increase clarity and contrast in the results and to avoid infinite small values in NDVI. In this study, has been replaced with . ECI is calculated to show the impact of SUHIs on urban areas, i.e. to show where surface heat islands have the greatest impact with respect to land cover.

RESULTS AND DISCUSION
Atmospheric correction of the Landsat 8 bands was conducted using the Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) (ENVI, 2009) in the ENVI 5.0 software package applying the data provided in the included metadata files. Vegetation indices (VI) were also calculated in the same software package along with principal component analysis (PCA). The land surface temperatures (LST) and environmental criticality indices (ECI) were calculated within QGIS using the Remote Sensing & GIS plugin (Barane, 2017).

Spatial Distribution of Land Surface Temperature
In this research VI, LST and ECI are calculated over Zagreb area for each month in 2017. Greater attention has been given to calculations of VIs. For this purpose, NDVIs, EVIs and RVIs were calculated for each season (Table 1) in 2017, and the PCA fusion of their data was performed. The first principal component (PC1_NER) was used in the LST calculation, in addition to Landsat 8 thermal bands.

Correlation
NDVI EVI RVI LST -0,28436 -0,36028 -0,44548 Table 3. Correlations between LST and NDVI, EVI and RVI The warmer areas are residential areas north (densely populated downtown) and south of the Sava River (densely populated newer part of the city), as well as densely populated areas on the slopes of Medvednica. Colder places are the areas along the Sava River, parks, peripheral rural areas in the south and north and the sparsely populated northern slopes of Medvednica. The visual interpretation reveals that the hottest area is the strict center of the city, which lies in the plain between the Sava River in the south and the Medvednica Mountain in the north.
For the purpose of accuracy estimation of the derived LSTs ((5), (6) and (7)), the comparisons of land and air temperature were measured in situ and from retrieved LST (Figure 4) on corresponding pixels, noting that the in situ LST values and air temperature were linearly modelled for the purpose of estimating the temperature in the moment of the satellite passing over Zagreb. The difference between in situ and retrieved LST is the greatest in summer (in situ values are higher) and the smallest in winter (retrieved values are greater) ( Figure 5). On the other hand, the largest differences between in situ measured air temperatures and retrieved LST values are the greatest in spring and the least in summer ( Figure 5).
A correlation between in situ measured LST and air temperature, and derived LST values was performed. The high correlations (between in situ and derivative LST is 0.974 and between in situ air temperature and derivative LST is 0.977, Figure 6) justify the use of the methodology based on thermal remote sensing data (Landsat 8) for the calculation of LST.

Spatial Distribution of Environmental Criticality Index
ECI (Senanayake et al. 2013a) was calculated as shown in the Equation (10) to identify the level of combined environmental criticality in the city of Zagreb on the basis of LST and existing vegetation cover. One LST was made for each season. This was done by averaging the LST over all months of one season. The same was done with calculated VI (PC1_NER). This is the input to calculating ECI for each season. The number of classes within the LST and ECI rasters depends on the climate, i.e. air temperature. Due to the values of temperature and condition of vegetation in each season, spring, summer and autumn have one common classification, while winter has a different classification of LST values. Air and soil temperatures reach maximum values during summer days, and vegetation is at its peak, so the contrast within the LST is the strongest. Due to the low temperature values, winter LSTs were not taken into account when calculating the final ECI value because they would be too averaged (this can easily be observed in Figure 3). Therefore, the final ECI was calculated by averaging the LSTs for spring, summer, and autumn (Figure 7). ECIs and LSTs were used to detect SUHIs. Given that the area of interest is an urban area dominated by the areas with buildings with little vegetation cover, ECI is a function of (mostly) LST and has the highest values in the city centre area. A look at the final ECI shows the supposed locations of SUHI in the City of Zagreb.

Hot Spots Determination
The LST is always lower in areas covered with vegetation and water bodies, and in some instances creates sharp boundaries between regions with high LSTs. According to LST and ECI maps of Zagreb, most of the SUHIs identified are located away from water bodies (the River Sava, Jarun Lake) and large forested areas of the northern slopes of Medvednica. The study reveals that the areas north and south of the Sava River have great LST values mainly due to dense build up area. One of the largest high-temperature zones is in the city centre ( Garden (in poligons) with the lowest ECI mainly due to anthropogenic land use. Dense concrete blocks of buildings are intersecting with asphalt streets (with a few trees) where heavy traffic is taking place (Figure 7). In the City of Zagreb, beside the city centre, the most common locations of SUHIs are in Žitnjak (industrial zone in the east), Jakuševac (garbage dump), Prečko and Jankomir (shopping centers and commercial buildings in the west). Larger areas covered by low vegetation and parks (Zrinjevac and Botanical Garden downtown (Figure 7), Maksimir in the east, Bundek in the central and Jarun in the western part of the city) can also be observed within the urban area with lower LST and ECI values. The ECI map was overlaid on Google Maps (Figure 7) to detect sites of increased surface temperature, or places of potential SUHIs.
The initial objective of this research for the purpose of detecting potential SUHIs and creating the basis for further research involving experts in the fields of meteorology, agronomy, forestry, sociology, to create a detailed picture of SUHI and measures to reduce this problem in the City of Zagreb was accomplished by making LSTs and ECIs maps.

CONCLUSION
In this study, locations of SUHIs and environmentally critical areas based on LST and ECI distribution and vegetation cover were identified in the city of Zagreb, Croatia. This was done with the integration of the daily satellite remote sensing Landsat 8 data and in situ measured air and land temperature on three meteorological stations in every month in 2017 (except for September and November which were cloudy). A comparison of in situ and satellite data showed that LST and ECI derived from the Landsat 8 spectral and thermal data are a good base for SUHI detection. Retrieved LST can evaluate urban surface temperature by quantity and in spatial patterns. The PCA method of image fusion (vegetation indices) was implemented to reduce redundant data from multiple VIs (NDVI, EVI, RVI) and to create a product that contains the highest percentage of variance of the entire set. Vegetation indices comparable to the NDVI, like EVI and RVI are useful measures in getting better perspective to the condition of vegetation cover and the relation to built-up areas and water bodies.
The Zagreb city centre, Jankomir in the western part and Žitnjak in the eastern part of the city were identified as the environmentally most critical region based on LST and ECI. Concrete and asphalt surfaces, roofs of buildings with low albedo and parking lots were identified as the major sources of high values of LSTs. On the other hand, the areas with lower LST values are also located, such as the areas along the banks of the Sava River, Parks Maksimir, Zrinjevac, Bundek and the Botanical Garden that break down the texture of high value LST areas within city borders.
Considering that the use of satellite remote sensing data allows a time and cost-effective methodology for detecting the SUHIs, the results of this and similar studies can be used for future urban development and planning the projects intended for the increase in vegetation areas in urban areas and reducing the impact of SUHI on the urban population. The continuation of the research is aimed at analysing LST and ECI over a long period of time in order to identify trends in the movement of LST and to involve experts in the field of meteorology, agronomy, forestry, civil engineering, sociology for the purpose of analysing other parameters affecting the value of LST.