A Multi-Temporal Analyses of Land Surface Temperature Using Landsat-8 Data and Open Source Software: The Case Study of Modena, Italy

The Urban Heat Island (UHI) phenomenon, namely urban areas where the atmospheric temperature is significantly higher than in the surrounding rural areas, is currently a very well-known topic both in the scientific community and in public debates. Growing urbanization is one of the anthropic causes of UHI. The UHI phenomenon has a negative impact on the life quality of the local population (thermal discomfort, summer thermal shock, etc.), thus investigations and analyses on this topic are really useful and important for correct and sustainable urban planning; this study is included in this context. A multi-temporal analysis was performed in the municipality of Modena (Italy) to identify and estimate the Surface Urban Heat Island (SUHI, strictly correlated to the UHI phenomenon) from 2014 to 2017. For this purpose, Landsat-8 satellite images were processed with Quantum Geographic Information System (QGIS) to obtain the Land Surface Temperature (LST) and the Normalized Difference Vegetation Index (NDVI). For every pixel, LST and NDVI values of three regions of interest (ROI, i.e., Countryside, Suburbs, and City Center) were extracted and their correlations were investigated. A maximum variation of 6.4 ◦C in the LST values between City Center and Countryside was highlighted, confirming the presence of the SUHI phenomenon even in a medium-sized municipality like Modena. The implemented procedure demonstrates that satellite data are suitable for SUHI identification and estimation, therefore it could be a useful tool for public administration for urban planning policies.


Background
About half of the world population lives in urban areas [1]. The global urbanization rate is expected to increase by 70% compared to the current world population [2], both because of the continued emergence of new urban areas [3] and because of the constant population migration from rural to urban and suburban areas [4,5]. It is not therefore surprising that the negative impacts of urbanization are an ever-growing global concern [6][7][8][9][10]. Urbanization has a negative impact on the environment, mainly due to pollution, changes in the physical and chemical properties of the atmosphere, and in the type of cover of the soil surface [11]. These phenomena lead to so-called Urban Heat Islands (UHI), namely urban areas where the atmospheric temperature is significantly higher than that in the surrounding rural areas [12]. The presence of UHIs is an increasing phenomenon studied by the international scientific community because of its dangerous and significant effects. In fact, the temperature increase has effects on the environment (higher temperatures cause higher energy consumption, photochemical smog, and worsening of the air quality), on the climate and on human health [13][14][15][16].
quite wet, especially in July and between the end of August and the beginning of September. These weather and climate conditions surely have an influence on the results obtained by remote sensing, because evapotranspiration influences the temperatures of urban and rural areas [57][58][59][60][61].
In this framework it is important to highlight that the municipality of Modena, in 2010, put in place a plan of action for sustainable energy (SEAP). This plan is included in the "good practices" for the European Green Infrastructure Strategies because it provides an increase in green infrastructures in agricultural areas, ecological zones, and public green. In particular, an increase in urban forest of 127.5 ha has been planned for implementation between 2011 and 2020 [62,63].
Thus, using the methodology developed in this paper, the city administrators will have a powerful tool to identify SUHIs and their intensity in different years. In this way, urban policies can be directed to those areas where the SUHI phenomenon is still present and creates discomfort in the population. Furthermore, the public administration could evaluate appropriate mitigation measures such as green roofs, cool materials, and green areas [64,65] that lead to a more sustainable life quality. In addition, using this methodology, it is possible to evaluate whether the policies entailed within the framework of the European Green Infrastructure Strategies are producing the desired results.

Materials and Methods
The analysis of the SUHI phenomenon in these areas was carried out using an experimental plugin for QGIS [66], called the Semi-Automatic Classification Plugin (SCP) [67], that allows use of the free open source software QGIS as a remote sensing software. With this plugin, temperature maps can be obtained from "raw" remote sensing images acquired by the Landsat-8 satellite sensor.
The study area includes the municipality of Modena located in the Emilia Romagna region, Italy. The environmental context is the Po Valley, which is characterized by a high population density and intensive agriculture. The LST was computed using the equation suggested by Weng [23]. Then, additional processing was conducted in order to create three Regions of Interest (ROI) within each image: the City Center region, the Suburbs region, and the Countryside region. For each ROI, LST maps and Normalized Difference Vegetation Index (NDVI) maps were retrieved. Creating these ROIs also allowed the estimation of the temperature difference between Countryside and City Center, and thus allowed us to correlate it with the SUHI phenomenon in the studied period (2014-2017) [68][69][70]. Finally, NDVI maps were used to analyze the correlation between land cover and LST maps.

Study Area
The study area is the municipality of Modena, Italy ( Figure 1). Modena has 184,826 inhabitants and an area of 183.19 km 2 [71]. Modena is located in the north of Italy, along the "via Aemilia", an ancient Roman road running from Rimini to Piacenza on the river Po ( Figure 2). The Po Valley is characterized by a high population density and by processes of industrialization and intensive agriculture (it is among the most productive agricultural areas within Europe) [72]. The climate of this region is partially continental: summers are hot with intense heat waves while winters are cold and wet, usually with atmospheric stability conditions and fog. Normally precipitation is concentrated during autumn and spring. Summer and winter are the driest seasons [57,58,72,73]. In order to have a clear climate situation of the study area, Table 1 reports some climate information. The data are divided into two series: an historical series from 1971-2000 and another time series from 2001 to 2017. For each series, Table 1 shows values of seasonal mean temperature, seasonal mean of the maximum daily temperatures, seasonal mean of the minimum daily temperatures, and seasonal mean of precipitation. The data referred to the City Center region were acquired by the Weather Station "Osservatorio Geofisico di Modena", located in the urban area of Modena (Long: 10 • 55 47.2 E, Lat: 44 • 38 52.9 N) [74], while the data referred to the Countryside region were acquired by the ARPAE Weather Station of Albareto (Long: 10 • 57 24 E, Lat: 44 • 42 7 N). In Table 1, all temperature values of the historical series increase in the time series in every season. The increase of the values vary between 1.4 • C and 2 • C. Also, the values of precipitation increase in every season except summer. In Modena, starting from the 1980s, climate change has exhibited the same trend as the global climate change. The changes cause extreme events such as jumps in temperature or increases in the distribution and amount of precipitation [57].
The area of Modena was chosen as a first area of interest for the proposed methodology for its characteristics in terms of:

•
Dimensions. Modena is a medium-sized municipality and the scientific literature does not investigate SUHIs of these kind of cities with remote sensing data; • Climate situation; • Geographical position. Modena is located in one of the most industrialized areas of Europe, as well as one with the highest levels of air pollutants [75]; • Active urban policies. Analyzing three years of data allowed us to investigate whether urban policies on green infrastructures are effective or not.
These characteristics render Modena an ideal test site for the developed methodology.

Landsat-8 Data
The Landsat-8 satellite was successfully launched on 11 February 2013 and deployed into orbit with two instruments on board: (1) the Operational Land Imager (OLI), whose spatial resolution is 30 m, with nine bands in the visible (VIS), the near infrared (NIR), and the short-wave infrared (SWIR) spectral regions; and (2) the Thermal Infrared Sensor (TIRS) with two spectral bands in the long-wave infrared (LWIR) region. The spatial resolution of TIRS data is 100 m with a revisit time of 16 days [75][76][77][78][79] (Table 2). The data used for this study are the "raw" remote sensing images acquired from the Landsat-8 satellite sensors, provided free of charge by the United States Geological Survey (USGS) ( Table 3). Images were divided into three seasons depending on mean LST values retrieved by satellite data: hot season (LST > 30 • C); mid season (10 • C < LST < 30 • C); cold season (LST < 10 • C).
As reported in Table 3, there are 15 daily images covering about the hot, mid, and cold seasons of four years. The scene center time is about 10:00 UTC. Only images with a clear-sky condition were selected.
There are only two images from 2015 because of the cloud cover, but winter 2015 can be represented by the image taken on 12 December 2014. Autumn 2015 can be represented by the image taken on 10 September 2015, which presents typical autumnal climatic conditions. Therefore, the year 2015 can be considered completely covered, except for spring.
For the entire procedure, the Semi-Automatic Classification Plugin (SCP) on the open source software QGIS was used. The SCP plugin allows one to use remote sensing functions in QGIS [22,24,67]. The methodology, from input data to obtained results, is shown in the conceptual diagram in Figure 3. Each step is described in detail in the following subsections.

Data Processing
a. Conversion to Top of Atmosphere (TOA) Radiance OLI and TIRS band data can be converted to top of atmosphere (TOA) spectral radiance using radiance rescaling factors provided in the metadata file [80]: where: • L λ = top of atmosphere spectral radiance (TOA) (W/(m 2 srad µm)). • M L = band-specific multiplicative rescaling factor from the image metadata (W/(m 2 srad µm)). • A L = band-specific additive rescaling factor from the image metadata (W/(m 2 srad µm)).
b. Conversion to TOA Reflectance OLI band data were converted to TOA planetary reflectance using reflectance rescaling coefficients provided in the image metadata files [80]. The following equation is used to convert DN values to TOA reflectance for OLI data: where: • ρλ ' = TOA reflectance. • M ρ = band-specific multiplicative rescaling factor from the metadata. • A ρ = band-specific additive rescaling factor from the metadata.
c. Conversion to At-Satellite Brightness Temperature TIRS band data were converted from spectral radiance to brightness temperature using the thermal constants provided in the image the metadata file [80]: where: • T = at-satellite brightness temperature (K). • L λ = TOA spectral radiance (W/(m 2 srad µm)). • K 1 = band-specific thermal conversion constant from the metadata (W/(m 2 srad µm)). This index is computed using spectral reflectance in the near infrared band and in the red band. It provides a rapid estimation of the presence of vegetation. NDVI values ranges from −1 to 1. Higher NDVI values indicate dense vegetation while lower values (typically from 0 to 0.2) identify light or dark soils [24]. NDVI is computed using the following equation: where: NDVI was used to calculate FVC, an index that estimates the proportion of an area covered by a set of predefined type of vegetation or soil cover [24].
The following equation is used to calculate FVC: where: • NDVI s = NDVI value for bare soils. • NDVI v = NDVI value for fully vegetated soils.
The NDVI s value was set equal to 0.1, while the NDVI v value was set equal to 0.65 [22].
g. Calculation of Land Surface Emissivity (LSE), which is a measurement of the capacity of a material to radiate energy [24] where: • ε s = typical soil emissivity (0.93).
h. Estimation of LST LST can be calculated from the at-satellite brightness temperature T B [23]. LST was first computed in K, then converted to • C. where: • λ = wavelength of the emitted radiance (11.5 µm for band 10 in Landsat-8 OLI).

LST Maps
The methodology was applied to the dataset in order to obtain LST maps of Modena for different climate seasons. In this section only a few significant images are shown of the 15 processed. The statistical results on the next pages were, however, retrieved using all of the images listed in Table 3. Figure 4 shows the LST map of the study area on 7 April 2014. The minimum value is 17.5 • C, the maximum value is 29 • C. Thanks to the spatial resolution of the TIRS sensor, inside the city in the industrial areas (LST values between 28 • C and 29 • C) it is possible to identify dense areas (most with LST values between 25.5 • C and 27.8 • C) and vegetation areas (LST values between 17.5 • C to 22 • C). Also, small residential areas and bare soils (LST values from 23 • C to 29 • C) could identified. Furthermore, main roads are clearly visible (LST values of about 23 • C). In the main urban area, LSTs are on average higher than in the countryside (the maximum temperature difference between the urban area and rural area is about 11 • C), confirming the presence of the SUHI phenomenon. In the following section, the difference between the urban LST and the rural LST will be further analyzed.
In winter, the LST maps show a more homogeneous distribution of temperatures and sometimes the urban LST values (from 0 • C to 2 • C) are slightly lower than the rural LST values (highly variable values that average between 0 • C and 5 • C). The map in Figure 5, for example, shows the study area on 9 January 2017. The situation is significantly different from that of Figure 4. LST values are very similar both in urban and in rural zones, thus in this period the SUHI phenomenon is not visible from this kind of data. In accordance with the scientific literature [61,81], during the winter period an inversion of temperature trends could occur (urban areas have lower temperatures than rural areas), especially in medium-sized municipalities like Modena. In order to explain this map, an investigation of the presence of water in the study area is required. When humidity levels are very high, the water evaporation and the plant transpiration saturate the atmosphere, thus countryside and city temperatures tend to flatten, out becoming much more similar to each other [61,81]. For this reason, in humid periods the SUHI is barely identifiable using remote sensing techniques, which however does not mean that the SUHI does not exist. For example, a temperature increase of just one or two degrees with humidity levels exceeding 90% does not allow the definition of the UHI phenomenon on an urban scale, but the population feels health discomforts nonetheless [82].

Region of Interest
To analyze temperature variations between urban areas and rural areas, three ROIs were created. A shapefile built thanks to municipality information was overlaid on each image in order to divide the municipality of Modena in three areas: City Center, Suburbs, and Countryside ( Figure 6). These ROIs correspond to three specific Local Climate Zones (LCZs) defined by Stewart and Oke as "regions of uniform surfaces cover, structure, material, and human activity that span hundreds of meters to several kilometers in horizontal scale" [83]. City Center corresponds to LCZ1 (Compact Midrise) characterized by a building surface fraction between 40% and 70%. Suburbs corresponds to LCZ5 (Open Midrise) characterized by a building surface fraction between 20% and 40%. Countryside corresponds to LCZD (Low Plants) characterized by a building surface fraction less than 10% [83]. The building surface fraction represents the ratio of building plan area to total plan area (%). For Modena's ROIs, the building surface fraction is equal to 46% for City Center, 20% for Suburbs, and 2% for Countryside. LCZs are defined in a univocal way based on specific parameters. The link between ROIs and LCZs allows the repeatability of the present methodology as well as comparison with other study areas.  Table 4 shows the average LST value of each ROI for every image and the difference (∆LST) between the City Center LST values (LST CC ) and the Countryside LST values (LST CS ). These values allow to have a first evaluation of the intensity of the SUHI phenomenon. Suburbs LST values are not considered for this estimation because they are often very similar to City Center LST values. Looking at Table 4, it can be observed that the images of cold seasons are usually associated with low LST variations (≤2 • C). Accordingly, with Figure 4 the SUHI phenomenon is not visible in these seasons with Landsat-8 data and thus it is not possible to correlate SUHIs with the UHI phenomenon.
Otherwise, images of hot seasons (spring and summer) are mostly associated with high LST variations (>2 • C). In this case it is possible to clearly identify the SUHI phenomenon and estimate its intensity. In Table 4 it is possible to observe the maximum and the minimum LST variation in the three seasons of interest. For the hot season, this value reached 6.4 • C, confirming the presence of the SUHI phenomenon and thus of the UHI phenomenon.
An exception of this analysis is represented by the image acquired on 9 January 2017. LST values for different regions are very close and the ∆LST between City Center and Countryside is lower than 2 • C. This anomaly is due to particular weather conditions of the days prior to the image acquisition; 2 and 3 July 2017 were in fact cloudy days with temperatures below the usual seasonal means. Figure 7 shows a graphical representation of mean LST values for each acquired image and for each ROI-during hot seasons, the SUHI phenomenon is clearly visible. The results highlighted in Table 5 and visually shown in Figure 7 could be compared to scientific studies on UHI even if the present study deals with SUHI. In Asian large cities, the UHI intensity ranges from 0.4 to 11 • C. These values are largely influenced by the monsoon circulation [40]. In Europe, in particular in the Mediterranean area, it is possible to find UHI intensity values from 3 • C to 7 • C for cities like Athens and Salonicco [44,45]. Focusing on small cities, two interesting studies show the presence of UHI even in cities with less than 100,000 inhabitants [84,85]. These studies revealed a mean summer UHI intensity from 2.6 to 8 • C.   Table 6 show the UHI intensity retrieved for Italian cities in the scientific literature during the hot seasons. Even if these values are related to UHI and not to SUHI and these cities exhibit some differences from Modena (size, position, and inhabitants), it is possible to notice that the results obtained from this study are comparable with the values in Table 6. In particular, Zauli Sajani et al. [51] analyzed a large area that includes Bologna as well as Modena: the UHI intensity values measured by the meteorological station during the summer period are coherent with the results obtained from the methodology presented in this paper. In Figure 7, the air temperature values are also reported. In particular, air temperature data for City Center and Countryside were retrieved from two weather stations provided by ARPAE Emilia Romagna, TCC and TCS, respectively.
These air temperature values are highly correlated to satellite LST values as it is possible to observe both in Figures 7 and 8. The Pearson coefficients (r) shown in Figure 8 and in Table 7 have values closer to 1, which means there is an excellent correlation between LST values and air temperature values.  The results shown in Figure 8 are very important because they prove the effectiveness and reliability of satellite measurements for the study of SUHIs and consequently of UHIs. The high correlation with air temperatures measured by weather stations demonstrates that the methodology is repeatable and can be transferred to other areas similar to the city of Modena (medium-sized municipalities). In urban areas with different dimensions (large cities or small municipalities), the methodology can still be applied, only requiring the foresight to choose carefully the ROIs. In particular, the number and size of these regions must be appropriate to the study area. In this framework, the reference to LCZs [83] allows, for cities of different dimensions, the identification of the correct ROI number and size to study the SUHI phenomenon. Weather stations must also be chosen accurately as they are representative of the investigated ROI. Possible influences from traffic or emissions from industrial activities must be avoided if possible.

NDVI Maps
NDVI maps were retrieved in order to partially identify surface land cover and correlate it with LST values [87]. Figure 9 represents the NDVI map of the "Modena 7 April 2014" Landsat image. High values of NDVI (0.6 to 0.9) suggest a high density of vegetation, and low values of NDVI (from 0 to 0.3) suggest the presence of bare soil or artificial surfaces. Obviously, urban areas present low NDVI values and rural areas have high NDVI values. From the NDVI map presented in Figure 9, it is possible to recognize built-up areas (or bare soil areas) and vegetation areas. NDVI maps are consistent with scientific literature data [28,88]. As reported by Yuan et al. [15], lower LSTs usually are found in areas with high NDVI. This negative correlation between NDVI and LST is valuable for urban climate studies. However, NDVI measurements are subject to seasonal variations, which may influence the results of SUHI studies.

LST-NDVI Correlation
For every pixel of the images from 7 April 2014, 15 July 2015, and 9 January 2017, the LST and NDVI values were extracted to make a correlation analysis.
Focusing on the image from 7 April 2014 and the image from 15 July 2015, Figure 10 shows that high values of NDVI correspond to low values of LST and vice versa. Thus, the presence of vegetation corresponds to lower LST values.
Regarding the 9 January 2017 image, it is difficult to make observations by just looking at the maps; it is necessary to analyze Table 8. This table reports the Pearson correlation coefficient between LST and NDVI computed for each ROI in the three images analyzed in Figure 10  For the images of April 2014 and July 2015, r values are close to 1, suggesting a strong correlation between NDVI and LST. Instead, for the image of January 2017, the r value is close to 0, thus a poor correlation between NDVI and LST is highlighted.  Compared to "Modena 7 April 2014" and to "Modena 15 July 2015", the January 2017 image has two main differences: the low LST-NDVI correlation and lower LST values in the city Center. This is probably due to the difference in thermal inertia of the materials mainly present in each ROI. When, the night before the image acquisition day, the air temperature is really low and there is a clear sky condition, the surface of the city becomes very cold. Thus, on the image acquisition day, more time is needed by the urban surfaces (mainly concrete roofs and asphalt) to increase their temperature compared to the countryside [59][60][61]. The countryside has more vegetation and bare soils compared to the city, therefore the temperature increase is faster here than in urban areas.

Conclusions
In this work the presence, extension, and intensity of the SUHI in the municipality of Modena were studied using Landsat-8 data.
The results showed that with satellite data, the SUHI is clearly visible during hot and mid seasons, while an opposite phenomenon is observed during cold seasons.
In particular, during hot seasons, Modena records a significant SUHI phenomenon, with a difference between City Center LST values and Countryside LST values of up to 6.4 • C.
Observing the entire studied period (from 2014 to 2017), the following assessments can be made: • Mean LST values are on average higher from 2014 to 2015 than from 2016 to 2017; • The SUHI phenomenon is more evident during hot seasons, especially in the years 2014 and 2015.
These assessments are certainly consistent with the policies implemented by the municipality of Modena in recent years concerning European Green Infrastructures Strategies. Thus, the application of the methodology presented in this paper highlights the mitigation effect on the SUHI (and therefore on the UHI) of the Modena SEAP agreement.
NDVI maps were useful to identify the vegetation cover and to correlate it with LST maps. This study confirmed that, for hot and mid seasons, high values of LST correspond to low values of NDVI (bare soil or artificial covers). For example, for the image of 15 July 2015, the Pearson coefficient (r) shows a correlation between LST values and NDVI values equal to 0.91 for City Center, 0.79 for Suburbs, and 0.77 for Countryside. For the mid season, the correlation between NDVI and LST is still strong. For example, for the image from 7 April 2014, the Pearson coefficient values are between 0.85 and 0.92. Otherwise, during cold seasons, this study revealed a null or opposite trend of the phenomenon and a low correlation between NDVI values and LST values.
Remote sensing LST values were also compared with air temperature values measured by two weather stations (one for the City Center region and one for the Countryside region). These values showed a high correlation with the Pearson coefficient equal to 0.98 for City Center and to 0.97 for Countryside. Remote sensing data are therefore representative of real surface temperatures and thus they can be used for the identification of the SUHI and the estimation of its intensity.
In conclusion, this study showed that the SUHI is a phenomenon that can also be observed in medium-sized municipalities like Modena and that can be investigated with a free and accessible methodology. The proposed methodology could be easily used not only for the identification of SUHI but also for its intensity estimation thanks to the good correlation between LST values and air temperature values. This study could thus represent a powerful tool for public administration for sustainable planning policies based on UHI mitigation strategies. Using the proposed methodology for several images acquired in different years allows one to monitor the UHI mitigation actions and verify their effect on the territory. Furthermore, the proposed methodology could be used not only for cities similar to Modena but also for different kind of cities, from large metropolises to small municipalities. For the transferability of the methodology, the choice of ROIs is important: each ROI has to be compared with LCZs definitions in order to set universal benchmarks to compare obtained results.
This study is not concluded. More images will be analyzed, in particular during hot seasons, and, for each image, the Albedo parameter will be calculated to correlate the SUHI to surface reflectance. Moreover, complete meteorological data will be collected in order to better understand the SUHI phenomenon, taking the local climate into consideration. Additionally, data such as the daily precipitation of the three days before the image acquisition time, the average daily wind speed, and the main wind direction will provide a complete climatic and environmental characterization of the acquisition time for each image.

Conflicts of Interest:
The authors declare no conflict of interest.