Satellite Images and Gaussian Parameterization for an Extensive Analysis of Urban Heat Islands in Thailand

Abstract

For the first time, an extensive study of the surface urban heat island (SUHI) in Thailand’s six major cities is reported, using 728 MODIS (MODerate Resolution Imaging Spectroradiometer) images for each city. The SUHI analysis was performed at three timescales—diurnal, seasonal, and multiyear. The diurnal variation is represented by the four MODIS passages (10:00, 14:00, 22:00, and 02:00 local time) and the seasonal variation by summer and winter maps, with images covering a 14-year interval (2003–2016). Also, 126 Landsat scenes were processed to classify and map land cover changes for each city. To analyze and compare the SUHI patterns, a least-square Gaussian fitting method has been applied and the corresponding empirical metrics quantified. Such an approach represents, when applicable, an efficient quantitative tool to perform comparisons that a visual inspection of a great number of maps would not allow. Results point out that SUHI does not show significant seasonality differences, while SUHI in the daytime is a more evident phenomenon with respect to nighttime, mainly due to solar forcing and intense human activities and traffic. Across the 14 years, the biggest city, Bangkok, shows the highest SUHI maximum intensities during daytime, with values ranging between 4 °C and 6 °C; during nighttime, the intensities are rather similar for all the six cities, between 1 °C and 2 °C. However, these maximum intensities are not correlated with the urban growth over the years. For each city, the SUHI spatial extension represented by the Gaussian footprint is generally not affected by the urban area sprawl across the years, except for Bangkok and Chiang Mai, whose daytime SUHI footprints show a slight increase over the years. Orientation angle and central location of the fitted surface also provide information on the SUHI layout in relation to the land use of the urban texture.

1. Introduction

Urban heating is usually quantified by means of the urban heat island (UHI), a well-known anthropogenic thermal modification whose effects have implications for human comfort and health, ecosystem function, local weather, air pollution, urban planning, and energy management [1,2,3,4]. The UHI phenomenon is linked to urban land transformations, especially when vegetated and rural zones are substituted by impervious surfaces, causing an increase in sensible heat flux and a decrease in latent heat in urban areas [5]. Urban heating is mainly ascribed to the trapping and absorption of solar irradiation in built-up areas connected to the thermal capacity of construction materials and the urban canyon phenomenon [6]. As a result, a difference in air and surface temperature between urban area and rural surroundings arises [7], quantified as UHI intensity or magnitude.

Many studies have estimated the UHI intensity by using ground-based weather stations observing air temperature in urban and rural locations [8], however, the sparsely distributed stations do not allow a global view over the urban area. Air temperature has a diurnal cycle with a nighttime UHI formation, while surface temperatures typically exhibit a surface urban heat island (SUHI) throughout the diurnal cycle, with a maximum in the daytime [9,10,11].

Advances in space-based remote sensing and satellite data availability have contributed to the growing number of studies on SUHI in several cities around the world, exploiting the global spatial coverage and the revisit time of satellite remote sensing missions [12]. SUHI map images were extensively retrieved by using Landsat data, providing spatiotemporal analysis and allowing the modeling of the urban thermal patterns with the land cover classes [13,14,15]. The MODerate Resolution Imaging Spectroradiometer (MODIS), aboard Terra and Aqua satellites, allows monitoring SUHI daily variations but with a lower spatial resolution (1 km) [16,17,18].

SUHI mapping can also satisfy the growing demand for integrating urban planning and landscape ecosystems, supplying scientific support for urban development policies [19]. In this context, remotely sensed imagery is often used to map impervious surfaces together with SUHI, the former widely recognized as a significant urbanization indicator affecting the thermal pattern. Instead, it is well known that the increase of vegetation and green spaces in the urban texture is one of the main strategies to mitigate SUHI intensity [20,21].

In this work, we provide a SUHI insight on different Thailand urban areas. Even though a great number of papers have explored the SUHI maps retrieved by remote sensing data in several cities worldwide, a complete study in Thailand is missing. For instance, different works were carried out with ground-based weather stations measuring air temperature, including: investigation of the UHI intensity in three big cities in Thailand (Bangkok, Chiang Mai and Songkhla) by using meteorological weather stations for the period 2004–2008 [22]; analysis of the air temperature data from four weather stations for a UHI insight in Bangkok area [23,24]. Spaceborne imagery in Thailand was used with limited spatial and temporal coverage: study of the land surface temperature (LST) pattern in Chiang Mai city (northern Thailand) with only two Landsat Enhanced Thematic Mapper Plus (ETM+) images (one in 2000 and one in 2006) [25]; Kantarawichai district analysis in the northeastern Thailand with only two Landsat images (one in 2006 and one in 2009) [26]; SUHI monitoring on Bangkok with few MODIS images in the period 2001–2003 [27]. Recently, a SUHI insight over the land use zoning plan of Bangkok using seven Landsat 8 scenes was proposed [28].

To bridge this gap, this work provides, for the first time, an extensive temporal and spatial investigation in Thailand by examining the SUHI of six major cities located in different geographic regions during a long-time interval (2003–2016) using 728 MODIS images, allowing daytime and nighttime monitoring. Also, 126 Landsat scenes were processed to classify and map land-cover changes in each city. Since Thailand has experienced significant changes in urban development during the last 20 years, this study points out the spatiotemporal SUHI trend for each city and intercity (trend among cities) and performs a comparison with the built-up surface growth derived from Landsat imagery.

To examine the SUHI patterns and their relevance in terms of temporal and spatial trends, empirical metrics were retrieved to quantify the phenomenon and to perform comparative analyses. The study was carried out by means of a two-dimensional least-square Gaussian fitting technique, whose parameterization and visualization proves to be useful in SUHI pattern insight.

2. Study Area

Six major cities in Thailand were selected for this study (Figure 1), including: (1) Bangkok (BG), the capital of Thailand, located in the central region; (2) Chiang Mai (CM), the largest city in the northern region, and four large cities located in northeastern region: (3) Nakhon Ratchasima (NR); (4) Ubon Ratchathani (UR); (5) Khon Kaen (KK); and (6) Udon Thani (UT).

Bangkok has a population of about 6.3 million within municipality and the largest urban extension in Thailand. It is located in the Chao Phraya River delta in Thailand central plains. The area is flat and low-laying, with a mean elevation of 1.5 m above sea level. It has a tropical wet and dry climate with three seasons—namely, summer, rainy season, and winter. The air temperature ranges from an average value of 22.0 °C in December to 35.4 °C in April. The rainy season begins around mid-May until October, with the dry and cool northeast monsoon occurring until February, when the summer season starts until mid-May.

Chiang Mai (136,700 inhabitants within municipality) is the capital of Chiang Mai Province, situated at about 300 m above sea level, 696 km north of Bangkok. About 70% of the Chiang Mai province consists of mountains covered with forests. The coolest months are December and January. The temperature throughout the year varies between 14 and 30 °C, with an average value of 26 °C.

Finally, four of the largest cities located in the northeastern Thailand were selected: Nakhon Ratchasima (135,300 inhabitants within municipality), Ubon Ratchathani (82,800), Khon Kaen (110,000), and Udon Thani (135,500). They are located mostly in flat plains between latitudes 14°0′E and 18°30′N and longitudes 101°0′E and 105°0′E. Average temperatures range from 19.6 °C to 30.2 °C. Northeastern Thailand is the focal point of the Greater Mekong Subregion, which operates as an economic center with Laos, China (Yunnan), and Vietnam. Development projects at large scales have sprung up, including construction of highways to provide access to the northeastern Thailand area for trade increasing. All of these aspects have expanded the northeastern economy and fostered rapid growth of the urban territories.

Figure 2 shows the urban surface evolution for the six cities in the years 2004, 2008, 2012, and 2016. These maps have been obtained from Landsat images as described in Section 3.1. The overview of Figure 2 shows the urban expansion across the years, which makes the study of the surface thermal impact interesting in terms of spatial and temporal evolution. Concerning the land use of urban areas for the six cities, the city core always includes high- and medium-density residential areas, commercial zones, and government institutes, while outskirts are mainly made up of medium- and low-density residential areas and industrial zones.

3. Data and Methods

3.1. Satellite Data

Land surface temperature (LST) images retrieved from observations of the MODIS sensors on board Terra and Aqua satellites were used. LST data were downloaded from ORNL DAAC archive [29], selecting MOD11A2 (Terra) and MYD11A2 (Aqua) 8-day LST products with 1 km resolution. These LST products are composed of daily 1 km LST values averaged during an 8-day period; cloudy pixels for each image were masked out. Terra satellite passes over Thailand at around 10:00 am and 22:00 (local time), whilst Aqua satellite passes at around 14:00 and 2:00 am; local time is +7:00 from UTC time. In this work, observations from 2003 until the end of 2016 were selected.

Since Thailand is in the tropical zone and has very few cloud-free images during the rainy/monsoon season, in this study, satellite LST images during two intervals within each year were acquired: three months during winter (November–January) and summer (February–April), excluding the six rainy months (May–October). Overall, 336 daytime and 392 nighttime MODIS-Terra/Aqua LST images were processed for a total of 728 raster LST maps of 1 km spatial resolution for each city.

Therefore, the whole dataset will allow an analysis of the SUHI temporal variation at three timescales: diurnal, seasonal, and multiyear. The diurnal variation is represented by the four times that MODIS passes (daytime: 10:00 and 14:00; nighttime: 22:00 and 02:00), while the seasonal variation is examined by comparing the summer and winter maps. Overall, all the selected images (24 in summer and 28 in winter for each city and for each year) were averaged to get seasonal LST mean values (i.e., winter and summer separately) at the four diurnal times, from which the SUHI was extracted.

To identify the rural buffer for the SUHI computation, MODIS Land Cover Type Product (MCD12Q1) with 500 m spatial resolution was used. The MCD12Q1 products yearly categorize the land using combined Terra and Aqua MODIS observations to identify different land cover classes.

Land use and land cover (LULC) maps for each city were also derived from cloud-free Landsat-4/Landsat-5 Thematic Mapper (TM), Landsat-7 Enhanced Thematic Mapper Plus (ETM+), and Landsat-8 Operational Land Imager (OLI) images with a spatial resolution of 30 m downloaded from USGS archive [30]. A total of 126 scenes from 2003 to 2016 were used to extract the urban area extent, exploiting the higher spatial resolution with respect to the MCD12Q1 products. The land covers were classified into four broad types (i.e., built-up land, water body, agriculture, and forest, Figure 2) applying a supervised maximum likelihood classification approach [31,32]. The user and producer accuracy [33] for urban and agriculture classes—evaluated by high-resolution images incorporated in GoogleEarth Pro^® and reported as mean value over all the Landsat maps—are as follows: user (urban) = 0.80; user (agriculture) = 0.94; producer (urban) = 0.95; producer (agriculture) = 0.84.

3.2. Estimation of the SUHI

The SUHI intensity quantifies the heating effects in an urban area, and we computed it as the surface temperature difference between urban (LST_URB) and rural (LST_RUR) pixels:

SUHI = LST_URB − LST_RUR.

(1)

The reference LST_RUR represents the mean of LST values associated with the rural pixels in a buffer area, surrounding the urban area, of about 3 km. In such buffer rural areas, water pixels were discarded. The choice of the 3 km buffer is a trade-off that allows, for all the cities, to discard high-altitude rural zones (e.g., Chiang Mai forest in Figure 2) that would lower the LST rural bias. For the computation of Equation (1), to discern built-up and rural pixels at the LST spatial resolution (1 km), urban and rural categories were extracted from the MCD12Q1 product. Figure 3 shows an example of these reference maps for the six cities for 2013, resampled at 1 km pixel size from the native 500 m, used for the computation of the reference LST_RUR.

The thermal characterization of an urban area by means of SUHI maps is also appropriate for comparative analysis among different cities because SUHI intensity is not an absolute temperature, but a temperature difference between rural and urban areas. Therefore, the effects of local weather conditions and other error sources are dampened.

To perform SUHI pattern comparisons, a technique based on a two-dimensional least-square Gaussian fit of the heat island was applied. This SUHI fitting technique was proposed in [34] and applied in [10,13,27,35,36,37]. The 2D Gaussian surface can be a useful representation of the SUHI pattern when an urban area shows a clear heat island with respect to the rural surroundings, especially if clustered SUHI values are prone to decrease towards the map boundary. The Gaussian fitting of the spatially distributed SUHI(x, y) (°C) obtained from MODIS data is described as follows:

$$\mathrm{SUHI}\left(\mathrm{x},\mathrm{y}\right)={\mathrm{a}}_0\times \exp \left[-\frac{{\left(\left(\mathrm{x}-{\mathrm{x}}_0\right)\cos \mathsf{\varphi }+\left(\mathrm{y}-{\mathrm{y}}_0\right)\sin \mathsf{\varphi }\right)}^2}{0.5{\mathrm{a}}_{\mathrm{x}}^2}-\frac{{\left(\left(\mathrm{y}-{\mathrm{y}}_0\right)\cos \mathsf{\varphi }-\left(\mathrm{x}-{\mathrm{x}}_0\right)\sin \mathsf{\varphi }\right)}^2}{0.5{\mathrm{a}}_{\mathrm{y}}^2}\right],$$

(2)

where (x, y) (km) is the pixel location in the map, with the origin (0, 0) in the center of the map. From Equation (2), a nonlinear data-fitting in least-squares sense was performed [35] in order to estimate:

-

the magnitude a0 (°C) of the SUHI maximum intensity of the study area;

-

the spatial extent (a_x, a_y) (km), central location (x₀, y₀) (km), and orientation (φ) (degree) of the SUHI pattern.

A correlation coefficient (r) between the SUHI map from MODIS data and the fitted Gaussian surface was computed to assess the capability of the fitting surface to follow the actual pattern. The two parameters, a_x and a_y (spatial extent), were also combined to define the overall area of the SUHI. Such an area or footprint is an ellipse obtained by horizontally cross-sectioning the Gaussian surface. The ellipse axes are a_x and a_y, φ is the orientation angle, with the maximum SUHI value located in the footprint center x₀, y₀.

The ellipse area defined by the distance from the center, where SUHI decreases to 61% (e^−1/2) with respect to the maximum value a0, will be defined as Ar61% (km²). Similarly, a footprint area where SUHI is greater than a fixed threshold can be defined—for example, greater than 1 °C (Ar1C) [10,34].

This fitting technique works well for urban areas with regular shapes. The method is suitable for large-area analysis by satellite data of spatial resolution like MODIS. Since the fitting provides a smoothed surface, the greater variability of higher spatial resolution images (e.g., Landsat maps) reduces the correlation coefficient and the representativeness of the Gaussian footprints [35].

4. Results and Discussion

The spatial trend of the SUHI pattern was analyzed at three timescales (diurnal, seasonal, and multiyear), considering the trend for each city and comparing the trend for all the six cities (intercity analysis). These SUHI patterns were examined during a 14-year interval (2003–2016): for example, the SUHI during 2013 for the six cities is shown in Figure 4 considering the summer season at 22:00.

The example of Figure 4 shows that the nighttime urban heating is not too pronounced—that is, few degrees with respect to rural area. In the next subsections, the use of the Gaussian parameterization will allow us to have a SUHI trend overview and to perform comparisons that a simple visual inspection of a great number of maps would not allow.

4.1. Gaussian Fitting of SUHI Maps

As discussed in Section 3.2, the 2D Gaussian surface can be a useful representation of the SUHI pattern when an urban area shows a clear and clustered heat island decreasing towards the rural area. Some situations (anomalous spotted rural heating or negligible temperature differences between LST_URB and LST_RUR) can lead to a low correlation coefficient r between the SUHI maps from MODIS data and the fitted Gaussian surface.

Table 1. Correlation coefficient r between SUHI from MODIS (MODerate Resolution Imaging Spectroradiometer) and Gaussian surface for the six cities. The mean r value over the years 2003–2016 is reported in brackets. x = reliable Gaussian fit; - = no Gaussian fit.

City	Season	10:00	14:00	22:00	02:00
Bangkok	Summer	x (0.70)	x (0.70)	x (0.75)	x (0.75)
	Winter	x (0.74)	x (0.73)	x (0.66)	x (0.64)
Chiang May	Summer	x (0.73)	x (0.73)	x (0.90)	x (0.87)
	Winter	x (0.71)	x (0.69)	x (0.86)	x (0.82)
Nakhon	Summer	-	-	x (0.72)	x (0.82)
Ratchasima	Winter	-	-	x (0.73)	x (0.79)
Ubon	Summer	x (0.60)	x (0.63)	x (0.87)	x (0.83)
Ratchathani	Winter	x (0.73)	x (0.69)	-	-
Khon Kaen	Summer	-	-	x (0.81)	x (0.77)
	Winter	x (0.73)	x (0.68)	x (0.82)	x (0.82)
Udon Thani	Summer	x (0.71)	x (0.70)	x (0.81)	x (0.84)
	Winter	x (0.88)	x (0.88)	x (0.67)	x (0.70)

describes the capability of the fitting surface to follow the actual SUHI pattern, reporting the r mean value for the years 2003–2016. We have considered the fitting reliable when is greater than 0.6 and the Gaussian footprint covers the urban area.

Table 1 shows that some combinations of city–season–hour do not have a robust fit. An explanation will be provided in the successive Figure 5 and Figure 6 by examining the correspondent SUHI patterns. The correlation is always better for the summer nighttime fitting with respect to the summer daytime, whereas the same is not confirmed for the winter. Considering the whole daily cycle, Chiang May has the best performance in terms of correlation.

Figure 5 and Figure 6 show the SUHI maps and the correspondent Gaussian fitting for the six cities during 2014 both in summer and winter, reporting also the diurnal variation at 10:00, 14:00, 22:00, and 02:00. Each SUHI(x, y) pattern from MODIS data—with x and y representing the latitude and longitude—were fitted with a Gaussian surface and the descriptive parameters collected. In Figure 5 and Figure 6, the black continuous ellipse superimposed on the maps represents the footprint area Ar61%, while the black dashed ellipse represents Ar1C. The red and green straight lines correspond to the directions of the ellipse’s major and minor axes, the former also defining the orientation angle φ. The intersection point of the two lines (i.e., the footprint center) is the central location (x₀, y₀), expressed as shift from the center of the map (0, 0). For example, with reference to the first panel of Figure 5 (BG Summer 10:00), Ar61% = 755.6 km², Ar1C = 1815.7 km², φ = 72.3° (i.e., along the north–east and south–west (NE–SW) quarters) and (x₀, y₀) = (7.3, −5.5) km; the estimated magnitude a0 is 3.4 °C. The Gaussian fitting procedure is not generally conditioned by localized surface temperature peaks within the pattern, being prone to smooth such isolated thermal discontinuities.

As a general overview of Figure 5 and Figure 6, the daytime SUHI is a more evident phenomenon with respect to the nighttime one, while summer and winter SUHIs at the same hour are quite similar. At nighttime periods, due to the absence of solar forcing and the reduced working activities and traffic, the urban heat island intensity decreases, since rural and urban temperatures tend to get closer. The exceptions are NR and KK, since their daytime SUHI intensities are similar to those at night in the urban area. Particularly, KK, one of the six cities with the smallest urban size, always exhibits a very weak urban heating with respect to rural surroundings. Ar61% and Ar1C, clearly distinct during daytime, are almost the same at night. In some cases, Ar1C is inside Ar61% for NR and KK, indicating a weaker SUHI with respect to the other cities. The biggest city, BG, shows the highest SUHI intensities. In the next Section 4.2 and Section 4.3, these aspects will be analyzed in depth for each city across the years.

Figure 5 and Figure 6 also provide an explanation of the lack of Gaussian fitting in some cases, as reported in Table 1. Concerning the daytime NR SUHI pattern (Figure 5), the urban fitting is not possible due to anomalous spotted rural heating in the south part of the map. Specifically, in the southeast zone, paddy rice areas are present, while in southwest areas, there are corn, sugarcane, and cassava croplands. Rice cropland experiences very different surface temperatures during the year that affect the LST of the rural area, since their thermal conditions depend strongly on site-specific activities and on the irrigation practice. Rice irrigation usually takes place during rainy season (May–October), with harvest finishing in December. Therefore, paddy rice areas are hot bare soils in summer during daytime, while in winter, the harvesting season, some areas start to appear warm. In the nighttime, such rural warming disappears, and a moderate urban heating prevails. Concerning the summer daytime KK SUHI pattern (Figure 6), the southwest rural surroundings are mostly paddy rice areas, providing the same rural warm effect of NR, even if it is not recognizable in winter daytime.

The daytime and nighttime behavior of the UR SUHI pattern is very different, since the urban area is divided in two parts by the river and its rural surroundings (see Ubon Ratchathani in Figure 2), as revealed by the cold strip in the middle of the daytime pattern of Figure 6. During daytime, the urban heating of the upper and wider part of the city is predominant. During nighttime, the urban pixels and the rural/river strip pixels between the two city parts have similar LST values, resulting in an overall wider city area with a more uniform LST pattern. This effect causes a nighttime SUHI footprint (with intensities less than 2 °C) that is wider and has a different location. In winter nighttime, the pattern does not allow a reliable Gaussian fitting. A similar nighttime enlargement of the SUHI footprint happens also to KK during winter.

For a more complete overview of the processed maps, an additional 48 maps for the six cities considering a daytime and nighttime summer map during the years 2004, 2008, 2012, and 2016 are reported in Supplementary Materials.

4.2. Fitting Parameter Comparison

The fitting parameters allow a quantitative intercity comparison of different SUHI features. We will show the footprint area Ar61%, magnitude a0, orientation angle φ, and central location (x₀, y₀) as mean values and standard deviations (std) computed over the years 2003–2016. In the figures below, the seasonal behavior (summer/winter) and the diurnal variation (10:00, 14:00, 22:00, and 02:00) are reported for each city.

Figure 7 shows Ar61%. For BG and CM, the footprint area during daytime is greater than during nighttime, without significant seasonality differences. In addition, the two daytime values (10:00 and 14:00) are similar, as are the two nighttime ones (22:00 and 02:00). The std values suggest a moderate parameter variability across the years.

As pointed out previously, the parameters for NR daytime, KK summer daytime, and UR winter nighttime were not evaluated. As expected, the three smallest cities (UR, KK, and UT) exhibit smaller Ar61% footprints. The increased nighttime area of UR and KK reflects the behavior described in Section 4.1.

Figure 8 shows the mean magnitude a0 over the 14 years. All the cities demonstrate a cooling effect at nighttime. The SUHI magnitude during daytime is always higher than at night, again without seasonality differences. The nighttime mitigation of SUHI is evident for BG and CM. Indeed, as for UR and UT, the winter daytime a0 can be greater than the summer one. This is mainly ascribed to the winter cooling effect of the rural buffer, which exerts a consequent increasing of the SUHI (Equation (1)), even if the winter urban LST is lower than in summer.

The magnitude peak is reached at 14:00 in BG with 6.5 °C both in summer 2012 and in winter 2011. For all the cities, the a0 results at 14:00 confirm an interesting phenomenon discussed in [38]. Daytime SUHI magnitude generally increases as the background LST rises, showing an intensity peak at 14:00. When background LST increases from 10:00 to 14:00, the urban area absorbs more radiation compared to near rural areas, and urban LST experiences a faster increase than the rural surroundings. The a0 results also point out the different underlying mechanisms between day and night SUHI, as discussed later, considering that the background LST is more correlated to daytime SUHI compared to nighttime, as proved in [38].

Mean and std of orientation angle φ are reported in Figure 9. This parameter depicts the layout of the Gaussian footprint in each map. The results highlight how the orientation keeps the same behavior for both the seasons, but exhibits differences between daytime and nighttime for some cities.

For instance, the daytime footprint orientation of BG is about 65° (i.e., along the NE–SW quarters) and follows the urban area constituted by residential and commercial zones, distinctly deployed in the NE–SW quarters [28,39]. These urban features, combined with the traffic congestion in the road network due to intense human activities, generate for BG the significant daytime SUHI pattern, both in summer and winter. During nighttime, the absence of solar forcing and the reduced working activities contribute to mitigation of the heating intensity and to modification of the SUHI footprint and orientation of BG, as also visible in BG panels of Figure 5.

The same behavior is found for CM, where the daytime direction of the SUHI footprint along the NE–SW quarters follows the layout of residential and commercial zones [40]. During the nighttime, there is a significant cooling effect of the urban core and outskirts, causing a reduction of Ar61%, a0, and a φ modification, as also noticeable in CM panels of Figure 5.

Figure 10 shows the central location (x₀, y₀) of the Gaussian surface expressed as distance from to the center of the map (0, 0); the map center roughly corresponds to the city center. Except for BG, whose urban extension is considerably wider than the other cities, the mean shift of the center of the Gaussian footprints (with respect to the map center) is small both seasonally and daily. For BG, (x₀, y₀) moves always towards southeast, with a displacement of the order of 5–10 km in daytime. This central location is positioned in the core of the high-density residential and commercial zone, also including government institutes and public utility zones [28,39]. At nighttime, the BG footprint shift is pronounced in the south direction, indicating the permanence of an urban heating towards the coastal zone. This fact could suggest a nocturnal sea influence, as discussed in [41].

4.3. Yearly Trend of SUHI Footprint and Magnitude

A comparative analysis of the footprint area Ar61% and magnitude a0 as a function of the 14-year interval (2003–2016) was performed to provide a specific-city and intercity overview of the SUHI trend. These parameters were compared with the urban area evolution across the years, computed from LULC Landsat maps (Figure 2) to analyze how the urban growth affects the SUHI pattern for both daytime and nighttime. Since previous results showed a weak seasonality dependence, Ar61% and a0 were averaged over summer and winter maps for each year.

The upper panels of Figure 11 show the yearly trend of the urban area (black line) computed from LULC Landsat data inside the boundary map of Figure 3, compared with the daytime (red line) and nighttime (blue line) Ar61% for BG, CM, and NR. Lower panels report the daytime (green line) and nighttime (magenta line) a0 values. The curves were obtained as shape-preserving interpolants over the 14-year points. The same is reported in Figure 12 for UR, KK, and UT.

For each city, the urban growth during the years is either weakly correlated or not correlated with the SUHI footprint trend, whereas no correlation is found with the a0 trend. Particularly, daytime A61% values for BG slightly increases from 660 km² to around 800 km² between 2003 and 2016, exhibiting a correlation coefficient r = 0.9 with the urban area trend. Daytime A61% values for CM increase from 42 km² to 54 km², but with r = 0.6 with respect to the urban increase. Instead, nighttime A61% for BG, CM, and NR are unrelated to city growth. For the three smaller cities (i.e., UR, KK and UT), the urban expansion does not affect the daytime and nighttime SUHI footprint.

Figure 11 and Figure 12 show how the magnitude a0 is decorrelated with the city sprawl across the years. As expected, maximum SUHI intensities are mainly driven by specific yearly thermal conditions regardless of urban growth. These results confirm how the spatiotemporal urban heating patterns depend strongly on the background climate, such as precipitation, air temperature, solar radiation, and wind speed [18].

BG shows the highest annual a0 during daytime, with values ranging between 4 °C and 6 °C. During nighttime, a0 is rather similar for all the six cities, between 1 °C and 2 °C, even if CM shows mean values around 2.4 °C between 2003 and 2010.

The day–night magnitude difference may be ascribed to different factors, such as the different solar radiation absorption between urban and rural areas, since urban surfaces have a larger heat capacity and can absorb and store large amounts of thermal radiation with a higher rate than the rural surroundings [42]. Furthermore, anthropogenic heat released from energy consumption is an important contributor to urban warming and it is positively correlated with the urban heat island intensity [43,44]. Heat release from energy consumption (e.g., air conditioning systems, industrial heat generation, and energy use in transportation) usually increases during daytime due to intense human and economic activities [45]. The activities and transportation in larger urban areas, such as BG and CM, are expected to generate greater amounts of anthropogenic heat release, contributing to the rise of the daytime a0 values for BG and CM with respect to the smaller cities.

The lack of correlation for Ar61% and a0 with the urban sprawl over the years, except for the slight increase of BG and CM daytime footprints, can be ascribed to different mechanisms underlying the SUHI phenomenon [46]. A possible reason is linked to the outskirt sprawl. Considering the LULC Landsat maps and overlapping the different yearly scenes of Figure 2, it emerges that the new outskirt areas have a lower density of built-up pixels than the city core, especially for the small cities.

To better understand how the growth or modifications of a city affect the SUHI patterns, different factors at large and fine spatial scales should be considered. Some works found that urbanization warmed the LST on an annual mean scale, regardless of urban area size [47,48]. The size influence on SUHI depends on urban compactness: stretched or disperse cities are preferable for urban heat island mitigation [6]. Outskirts and suburban areas, if not densely urbanized, contain wide rural patches containing soil moisture, thus helping reduce the warming rate of the suburban zones [49]. Therefore, the results in Figure 11 and Figure 12 confirm the background climate influence on the heating trend over the years and suggest that the weak SUHI–urban growth relationship is ascribed to a lack of compactness in the urban sprawl, especially for the smaller cities investigated.

It is interesting to highlight other specific drivers of SUHI arising in the urban expansion and modification. Urbanization causes land cover changes, but the analysis of the urban heat island pattern should consider the heterogeneity of the new urban surfaces (materials and uses) and associated LST dynamics, as reported for Bangkok in [28]. Land cover changes from rural to urban affect the daytime and nighttime SUHI in different ways, depending on the rural category (grassland, cropland, shrub, forest) [50]. In the study cases, the rural surroundings exhibit a marked land cover variability with evident seasonal surface thermal changes.

Factors related to the urban texture at district level that affect the SUHI phenomenon—driving the heat trapping and absorption—rely on surface albedo, material, and color; on building layout, slope, and height; and on greenery cover [51,52,53,54]. At the same time, the urban texture can intensify the warming effect by reducing the wind circulation [55]. These factors suggest how the relationship SUHI–urban growth can show different trends and correlations, depending on the city specificity.

5. Conclusions

The SUHI pattern can be influenced by different factors, depending on external (location, climate) and intrinsic (city size and sprawl, land use, human activities) features. This work shows that the biggest cities, CM and BG, have higher daytime SUHI intensities, with pattern shape and layout related to the urban land use (e.g., high-density residential and commercial zones). At nighttime, SUHI intensities decrease for all cities and the heating pattern features generally change. This analysis highlights how the background climate and the intrinsic factors—such as human activities, transportation, and land use—contribute to the SUHI increase, especially at daytime and for the wider cities. Instead, the SUHI pattern is not correlated with the urban growth over the years: if the city sprawl is not compact and the new areas have low built-up density, the urban growth influence on SUHI is very weak.

The study of the SUHI temporal–spatial trend provides scientific support for urban land use policy to regulate and plan the intrinsic factors. For this purpose, the use of remote sensing data is unavoidable. Whereupon, suitable tools to process remote sensor data can be useful to perform quantitative analysis and comparisons of urban heating patterns, supplying indications that can be specific for a city but also can represent general signatures.

This study has also pointed out the usefulness of MODIS data for temporal sampling, but the low spatial resolution, suitable for Gaussian fitting, is a limit hindering SUHI investigation at district level. The possible options (Landsat and aircraft thermal images) have restricted or sporadic overpass times. These technical features motivate the future development of satellite constellations with the aim to improve, simultaneously, spatial and temporal resolutions.