Determining the Boundary and Probability of Surface Urban Heat Island Footprint Based on a Logistic Model

Abstract

Studies of the spatial extent of surface urban heat island (SUHI or UHI_Surf) effects require precise determination of the footprint (FP) boundary. Currently available methods overestimate or underestimate the SUHI FP boundary, and can even alter its morphology, due to theoretical limitations on the ability of their algorithms to accurately determine the impacts of the shape, topography, and landscape heterogeneity of the city. The key to determining the FP boundary is identifying background temperatures in reference rural regions. Due to the instability of remote sensing data, these background temperatures should be determined automatically rather than manually, to eliminate artificial bias. To address this need, we developed an algorithm that adequately represents the decay of land surface temperature (LST) from the urban center to surrounding rural regions, and automatically calculates thresholds for reference rural LSTs in all directions based on a logistic curve. In this study, we applied this algorithm with data from the Aqua Moderate Resolution Imaging Spectroradiometer (Aqua/MODIS) 8-day level 3 (L3) LST global grid product to delineate precise SUHI FPs for the Beijing metropolitan area during the summers of 2004–2018 and determine the interannual and diurnal variations in FP boundaries and their relationship with SUHI intensity.

1. Introduction

Urbanization results in land use change from non-urban to urban land, and is accompanied by increases in anthropogenic heat release [1,2,3]. Consequently, city energy and water balances are altered through the reduction of latent heat flux and the rise in sensible heat flux, which further affect regional and even global climates [4,5,6,7]. Among these impacts, the urban heat island (UHI), a phenomenon characterized by significantly higher air and land surface temperatures (LSTs) in urban areas than in suburban areas, negatively affects the health and comfort of urban dwellers and greatly increases energy consumption [8,9,10,11,12]. The UHI can be divided into four types, subsurface urban heat island (UHI_Sub), surface urban heat island (SUHI), canopy layer urban heat island (UHI_UCL), and boundary layer urban heat island (UHI_UBL), according to temperature differences in urban and rural cooling and warming rates at the surface, in the substrate and in the air [13]. In general, SUHI for cities or regions is easily estimated by means of satellite and aircraft sensors, which could be used to retrieve the differences between the interface of the outdoor atmosphere with the solid materials of the city and equivalent rural air to ground interface [13]. More importantly, the area affected by a SUHI is much larger than the area of the city [14,15,16,17]. Therefore, SUHI has gained considerable research interest during the past several decades [18,19,20,21,22,23,24].

UHI intensity (UHII), a well-known indicator of SUHI, is calculated as the temperature difference between urban land and a reference rural region [14,25,26,27,28,29,30,31,32]. Although UHII quantitative and spatiotemporal characteristics, as well as SUHI sources/sinks [8,33,34,35], have been characterized for many cities worldwide using multiple algorithms [14,25,26,27,28,29,30,31,32], in-situ station observations [36,37,38,39,40], and remote sensing techniques [41], the results of these studies lack stability and generalizability. The magnitude and extent of the SUHI can be biased by the selected reference rural regions, including a few specific observation stations or pixels, one or several types of land use, all areas within the buffer at a certain distance from the city boundary, etc. [42,43]. Therefore, accurate quantification of SUHI requires the selection of appropriate reference rural regions and relevant temperature thresholds for these regions.

The footprint (FP) of the UHI effect is a relatively new index that quantifies the area affected by the SUHI, by defining the extent of increased temperature with respect to the reference rural region [17,44,45,46,47]. This comprehensive index considers both the magnitude and spatial extent of the SUHI effect [17,46]. Gaussian surface models have been widely applied to calculate the FP due to its good performance in SUHI modeling [46,48,49,50]. A single exponential decay model has also been used to extract the FP by examining obvious urban/rural temperature “cliffs” [17]. In these previous studies, topographical factors and the SUHI sink (e.g., water or green land) were determined to eliminate the effects of abnormally low-temperature objects on the SUHI [35,51]. However, most studies have adopted the shape of the city as a default basis for the SUHI FP, with buffer zones being stretched from the border of the urban area [17]. These approaches appear not to have accounted for differences in the effects of land use types or landscape heterogeneity on LSTs outside of the urban area, and can impact the SUHI FP extent in multiple directions [52,53,54,55,56]. Therefore, the objective of this study was to develop a new method to determine LSTs threshold for reference rural regions in all directions, with the aim of determining a more precise boundary of the SUHI FP in a metropolitan region. The proposed new method will be more reasonable than previous methods in that it minimizes the effects of the city’s shape, topography, or landscape heterogeneity.

2. Materials and Methods

2.1. Study Area

The study area is the Beijing metropolitan region, China (Figure 1). This metropolitan region was selected for several reasons. The city has a sub-humid warm temperate continental monsoon climate and four distinct seasons, with a hot and humid summer. The phenomenon that there is an obvious summer SUHI intensity in Beijing has been widely observed through multi-source remote sensing data [35]. Meng et al. pointed out that the maximum value of the mean SUHI intensity can reach 5.88 °C and 2.46 °C in day and night, respectively, during 2003–2015 [57]. Mountains and plains are distributed to the northwest and southeast of Beijing, respectively, which is important for eliminating the effects of terrain on the SUHI FP and validating the proposed method.

2.2. Data Source

The moderate resolution imaging spectroradiometer (MODIS), which has fine spatiotemporal resolution and data accuracy, is widely used in the study and simulation of SUHI processes and mechanisms [11,12,35,58]. In this study, we used the Aqua MODIS 8-day level 3 (L3) LST global product (MYD11A2), which has a spatial resolution of 1 km, for the summers of 2004–2018 (June, July, and August). Its overpass times are around 13:30 (local solar time) in descending mode and 1:30 in ascending mode, representing daytime and nighttime LST observations, respectively. The dataset includes 12 data periods per year. Each period can be further divided into daytime and nighttime LST. The attached quality control (QC) files were used to exclude low-quality pixel values affected by clouds or other noise.

Land use datasets were used to remove pixel values affected by cold sources (e.g., ponds, lakes, rivers, parks, or golf courses) in urban areas. The 30-m-resolution dataset of the period of 2005 and 2015 was manually interpreted via Landsat Thematic Mapper (TM)/Enhanced Thematic Mapper (ETM) images at the Chinese Academy of Sciences and was found to have high classification accuracy based on massive field sampling [59]. The Beijing metropolitan region comprises six types of land use, including cropland, forest, grassland, water, urban land, and rural settlement in the dataset. The land use of 2017 with 10-m-resolution comes from Finer Resolution Observation and Monitoring of Global Land Cover (http://0-data-ess-tsinghua-edu-cn.brum.beds.ac.uk/fromglc2017v1.html) [60]. The dataset includes cropland, forest, grassland, shrubland, wetland, water, impervious surface, and bare land in the Beijing metropolitan region.

A digital elevation model (DEM) with a spatial resolution of approximately 90 m was downloaded from the Space Shuttle Radar and Topography Mission (SRTM) and used to take into account terrain effects [12].

2.3. Method

To quantify UHII, we first defined the reference rural LST background field. The LST background field must be relatively stable and unaffected by the urban area. However, depending on the topography and landscape heterogeneity, the selection of different reference points can result in different LST background fields and SUHI intensity values [37]. We assumed different reference points for different directions, setting the urban center (or center of gravity) as the polar coordinate origin. Thus, the LST of each pixel and the reference rural LST in a certain direction can be expressed as T_S(θ,r) and T_ref(θ). The UHII for each pixel, noted as UHII(θ,r), can thus be calculated as follows.

$$\mathrm{UHII}\left(\mathsf{\theta }, \mathrm{r}\right)={\mathrm{T}}_S\left(\mathsf{\theta },\mathrm{r}\right)-{\mathrm{T}}_{ref}\left(\mathsf{\theta }\right),$$

(1)

We first determined the center (or center of gravity) of the study area. Buffer rings were then drawn from the urban center to the surrounding rural region (Figure 2).

To ensure equal sampling, we used buffer rings of equal area. In repeated experiments, the initial radius r₁ was set to 5 km, with each subsequent ring area S_i equal to the initial ring area S₁. Thus, the radius r_i of each subsequent ring was determined as $5\sqrt{\mathrm{i}}$ km, where i is the number of the buffer ring. The average LST within each buffer ring was then calculated as follows.

$$\overline{{\mathrm{T}}_S}\left(\mathrm{r}\right)=\frac{1}{2\pi }{{\displaystyle \int }}_0^{2\pi }{\mathrm{T}}_S\left(\mathsf{\theta },\mathrm{r}\right)\mathrm{d}\mathsf{\theta },$$

(2)

Theoretically, the air temperature should follow a relatively well-defined spatial pattern radiating from the urban center to the suburbs [17,47]. Although the LST would show a certain fluctuation because of the surface properties, including geometric, radiative, thermal, moisture and aerodynamic, during the radiation process (see the blue points in Figure 2), the overall performance of the fitted logistic curve (see the red line in Figure 2) could reflect the similar decline trend compared to air temperature [61,62,63]. When these regions are far enough away from the urban center, the LST should remain relatively stable as it is no longer affected by the SUHI effect. During the process of LST decline, abnormal values can occur when LST is affected by topography, small cold sources (e.g., ponds, lakes, rivers, parks, and golf courses), or clouds. An algorithm must be developed to eliminate these abnormal values prior to the fitting process.

$$\overline{{\mathrm{T}\prime }_S}\left(\mathrm{r}\right)=\frac{1}{2\pi }{{\displaystyle \int }}_0^{2\pi }{\mathrm{T}}_S\left(\mathsf{\theta },\mathrm{r}\right)·\mathrm{g}\left(\mathrm{DEM},\mathrm{water},\mathrm{greenland},\mathrm{cloud}\right)\mathrm{d}\mathsf{\theta },$$

(3)

After the removal of outliers, the logistic model was finally applied to fit the T_ref(θ) curve.

$${\mathrm{T}}_{ref}\left(\mathsf{\theta }\right)=f\left(\overline{{\mathrm{T}\prime }_S}\left(\mathrm{r}\right)\right),$$

(4)

The blue points in Figure 2 indicate average LST values the buffer rings, extracted after the outlier removal. These LST points were fitted to the smoothed LST series curve (red solid line) using the logistic algorithm. This algorithm was applied to characterize the LST series using the most important transition point derived from the LST spatial series data, i.e., T_ref(θ), which is the location where the LST becomes relatively stable. The peripheral area at this location is no longer the SUHI FP.

The change in the LST data for a single urban structure can be modeled as follows:

$$f\left(\overline{{\mathrm{T}\prime }_S}\left(\mathrm{r}\right)\right)=\frac{c}{1+e^{a+br}}+\mathrm{d},$$

(5)

where r is the buffer ring number, $f\left(\overline{{\mathrm{T}\prime }_S}\left(\mathrm{r}\right)\right)$ is the LST value at buffer ring r, a and b are the parameters of the logistic curve that require fitting. The sum of c and d represents the highest LST value among all buffer rings, and d is the lowest LST value among all buffer rings.

Mathematically, these transition points can be calculated using the rate of change of the curvature (K) curve (Figure 2, black dotted line) of the fitted LST series curve. Values of K for Equation (5) at different buffer rings can be calculated as follows:

$$K=\frac{d_{\beta }}{d_s}=-\frac{b^{2}cz\left(1-z\right){\left(1+z\right)}^3}{{\left[{\left(1+z\right)}^4+{\left(bcz\right)}^2\right]}^{3/2}},$$

(6)

where z = e^a+br, β is the angle of the unit tangent vector at buffer ring r along a differentiable curve, and s is the unit length of the curve. The K’ value, which represents the rate of change of the curvature curve, can be calculated as follows.

$$K^{\prime }=b^3\mathrm{cz}\left\{\frac{3z\left(1-z\right){\left(1+z\right)}^3\left[2{\left(1+z\right)}^3+b^2{c}^{2}z\right]}{{\left[{\left(1+z\right)}^4+{\left(bcz\right)}^2\right]}^{5/2}}-\frac{{\left(1+z\right)}^2\left(1+2z-5z^2\right)}{{\left[{\left(1+z\right)}^4+{\left(bcz\right)}^2\right]}^{3/2}}\right\},$$

(7)

As the LST decreases from the urban area to the surrounding area, three extreme values of K’ in the LST fitted curve can be inferred from Equation (7). The last minimum point (black) corresponds to the location of the reference rural region. T_ref(θ), the average LST in the reference ring, was identified as the reference rural LST, and the difference between the maximum LST and T_ref(θ) is the UHII at an angle of θ. Thus, each pixel can be considered part of the SUHI FP based on whether the LST in the pixel is greater than the threshold T_ref(θ) (Figure 3).

To avoid the impact of remote sensing data instability, we further explored the probability of seasonal SUHI FP for each pixel via Equation (8), such that higher probability indicated a greater likelihood that the pixel was part of the SUHI FP during that season.

$$Pro{b}_{FP}\left(\mathsf{\theta },\mathrm{r}\right)=\frac{Nu{m}_{FP}\left(\mathsf{\theta },\mathrm{r}\right)}{Nu{m}_{logistic}\left(\mathsf{\theta },\mathrm{r}\right)}\times 100\%,$$

(8)

where $Pro{b}_{FP}\left(\mathsf{\theta },\mathrm{r}\right)$ means that the probability of seasonal SUHI FP for each pixel. $Nu{m}_{FP}\left(\mathsf{\theta },\mathrm{r}\right)$ means that the number of the pixel being judged as the SUHI FP in summer. $Nu{m}_{logistic}\left(\mathsf{\theta },\mathrm{r}\right)$ means that the number of the pixel passed the significance test through the logistic fitting. For the 8-day MODIS LST product, the value of $Nu{m}_{logistic}\left(\mathsf{\theta },\mathrm{r}\right)$ should range between 0 and 12. This probability was calculated only if the logistic fit was found to be significant within the 8-day MODIS LST product.

UHI capacity was introduced to characterize the degree of UHI effect in the UHI FP.

$$C_{UHI}=\underset{D}{\overset{}{{\displaystyle \iint }}}UHII\left(\theta ,r\right)d\sigma ,$$

(9)

where $C_{UHI}$ is UHI capacity, i.e. the cumulative amount of UHII in the UHI FP (shown in the shaded part in Figure 2). The unit is °C·km². D is the UHI FP area. $UHII\left(\theta ,r\right)$ is the UHII in each pixel.

3. Results

3.1. Spatiotemporal Variation in the FP

3.1.1. Daytime Spatiotemporal Variation in the FP

Spatially, high-probability (>70%) summer SUHI FPs were mainly concentrated in urban centers, including Dongcheng, Xicheng, Haidian, Chaoyang, Shijingshan, Fengtai Districts, the northern edge of Daxing District, and the Capital International Airport in Shunyi District (Figure 4). The summer SUHI FPs showed significant expansion between 2004 and 2018, especially from the northern edge to the center in Daxing District, from the southern edge to central Changping District, and outward from the Capital Airport. Even in the southern parts of Miyun and Huairou Districts, there were large patches of high-probability summer SUHI FPs in the last few years of the study period.

Medium-probability (30–70%) summer SUHI FPs were distributed at the periphery of the high-probability summer SUHI FPs, mainly concentrated in the eastern part of Fangshan District, the southwestern part of Daxing District, the central parts of Tongzhou and Changping Districts, the northern part of Shunyi District, and the southern parts of Miyun and Huairou Districts. Distinct expansion of medium-probability FPs was observed in the eastern part of Fangshan District, the southwestern part of Daxing District, and the southern part of Tongzhou District, and there were clear conversions from medium-probability to high-probability FPs in the central part of Changping District, the northern part of Shunyi District, and the southern parts of Miyun and Huairou Districts.

Low-probability (<30%) summer SUHI FPs were mainly distributed in Yanqing, Huairou, Miyun, Pinggu, Mentougou, and Fangshan Districts at the periphery of the Beijing metropolitan area. There was no significant spatiotemporal evolution between 2004 and 2018.

Low-probability summer SUHI FPs had the largest area, followed by medium-probability and high-probability FPs (Figure 5). A 5.5:1:1 ratio was observed among the three types of FPs. Overall, low-probability SUHI FPs remained generally stable, accounting for approximately 72% of the total Beijing metropolitan area. The area of medium-probability FPs was basically the same as that of high-probability FPs, but with different trends during the period 2004–2018. The proportion of medium-probability SUHI FPs declined from a maximum of 18.96% in 2007 to a minimum of 9.85% in 2016. The proportion of high-probability SUHI FPs increased from a minimum of 9.64% in 2005 to a maximum of 18.77% in 2017. These increases mainly resulted from the highest-probability (>90%) SUHI FPs. These areas accounted for most of the high-probability SUHI FPs.

3.1.2. Nighttime Spatiotemporal Variation in FPs

Spatially, high-probability (>70%) summer SUHI FPs were also concentrated in urban centers, including Dongcheng, Xicheng, Chaoyang, Shijingshan, and Fengtai Districts, the eastern part of Haidian District, the northeastern part of Fangshan District, and the northern part of Tongzhou District (Figure 6). The area of high-probability summer SUHI FPs was significantly smaller than that in the daytime. Notably, there were no high-probability summer SUHI FPs in the northern part of Daxing District or the northern and central parts of Haidian District. However, a large portion of the summer SUHI FPs appeared around the Miyun Reservoir in the nighttime. From 2004–2018, the nighttime high-probability summer SUHI FPs showed a weaker expansion trend than that in the daytime. There were four significant expansions of nighttime high-probability FPs. The largest of these was located in the northern part of Daxing District, with additional contiguous expansion between Fengtai and Fangshan Districts, and between Chaoyang District and the Capital International Airport in Shunyi District. Notably, high-probability summer SUHI FPs expanded around the Miyun Reservoir during the nighttime.

Medium-probability (30–70%) summer SUHI FPs were also distributed at the outer periphery of the city, especially in the northern parts of Tongzhou and Daxing Districts near the urban area. There was also a large medium-probability summer SUHI FP in the central to the southern part of Fangshan District. From 2004–2018, medium-probability summer SUHI FPs showed the most significant expansion trend, even compared with the overall daytime trend. Medium-probability summer SUHI FPs sprawled rapidly to the south in Tongzhou, Daxing, and Fangshan Districts. There was also a contiguous northwestern expansion trend in Haidian and Changping Districts.

Low-probability (<30%) summer SUHI FPs were the largest areas distributed far from the city center, with no significant spatiotemporal evolution except for an area encroached by medium-probability summer SUHI FPs between 2004 and 2018.

Nighttime low-probability SUHI FPs were larger than those in the daytime, accounting for 76.69% of the total Beijing metropolitan area during the study period, reaching 84.13% in 2005 (Figure 7). This result was mainly due to the decreased proportion of high-probability SUHI FPs compared to that in the daytime. The proportion of high-probability SUHI FPs was 9.80% between 2004 and 2018, 4.10% less than that in the daytime. The proportion of medium-probability SUHI FPs was similar to that in the daytime. The area of medium-probability and high-probability SUHI FPs increased, whereas the area of low-probability SUHI FPs decreased from 2004–2018. The proportion of low-probability SUHI FPs declined from a maximum of 84.13% in 2005 to a minimum of 66.33% in 2017. This loss of area was mainly due to encroachment by medium-probability SUHI FPs, which increased in area from a minimum of 7.85% in 2005 to a maximum of 21.81% in 2017. High-probability SUHI FPs showed a relatively stable growth trend. The highest-probability (>90%) SUHI FPs accounted for most of the high-probability SUHI FPs.

3.2. Comparative Analysis of SUHI FPs Obtained by Gaussian Surface Model and Logistic Model

To validate the results from the proposed logistic method, we compared the derived UHI FPs with those calculated from the Gaussian surface model, which is the most commonly used method in the UHI FP literature [46]. Specifically, we used the 8-day SUHI FP in the summer of 2017 as an example of the comparative analysis.

Overall, the spatial patterns of SUHI FPs based on the two models are different, as all the SUHI FPs estimated by the Gaussian surface model were ovals that are continuous in space (Figure 8). For the Gaussian surface model, although the centers of the ovals appear to be near the city center, their morphological distributions greatly vary, as affected by the unstable LSTs retrieved from the remote sensing images. However, the SUHI FPs by the logistic method does not show such oval patterns.

Table 1. The areas of the probability levels of SUHI FPs estimated by the Gaussian surface model and the logistic model (Unit: km²).

The Probability Levels of SUHI FPs	Gaussian Surface Model		Logistic Model
	Daytime	Nighttime	Daytime	Nighttime
0.0–0.1	9823	6723	8128	5180
0.1–0.2	712	1678	934	2557
0.2–0.3	567	1083	672	1983
0.3–0.4	329	1014	614	1237
0.4–0.5	436	683	523	772
0.5–0.6	259	357	525	658
0.6–0.7	202	307	508	529
0.7–0.8	156	172	563	610
0.8–0.9	166	256	612	808
0.9–1.0	2004	2381	1575	320

provides the areal percentages by the probability level of SUHI FP between the two methods. The areas of high-probability SUHI FPs estimated by the two models are almost identical in the daytime, accounting for 15.87% and 18.77% of the total area, respectively. In the nighttime, however, the proportion of high-probability SUHI FP by the Gaussian surface model (19.17%) was significantly higher than that by the logistic model (11.86%). Moving onto the medium probability, the areas of medium-probability SUHI FPs by the Gaussian surface model, for both daytime and nighttime, were substantially smaller than that by the logistic algorithm. In particular, the SUHI FP estimated by the Gaussian surface model in the daytime only accounted for 8.37% of the total study area, which was the smallest type of SUHI FP and was only half of the SUHI FP by the logistic algorithm. Last, for low probability, the nighttime SUHI FP areas by the two models are almost identical, accounting for 66.43% and 66.33% of the total area, respectively. However, in the daytime, the proportion of low-probability SUHI FP by the Gaussian surface model was 9.34% higher than the low-probability SUHI FP calculated by the logistic algorithm.

To further analyze the reasons for the differences in SUHI FPs between the Gaussian surface model and the logistic model, we calculated the area differences in different probability levels of SUHI FPs by cross-classification tables (Tables S1 and S2 in Supplementary Information).

In the daytime, among all the probability ranges, the SUHI FPs with the highest (0.9–1.0) and the lowest (0–0.1) have the highest agreement of recognition between the two models, with their agreement percentages reaching up to 80.51% and 64.60%, respectively. However, for the other high-probability ranges, namely 0.8–0.9 and 0.7–0.8, only 6.21% and 3.20% of the SUHI FPs recognized by the logistic model were identified as the same level by the Gaussian surface model, with most being classified to the highest level (i.e., 0.9–1.0) by the latter due to the continuity of its shape.

In the nighttime, the distribution of SUHI FP probabilities was generally consistent with the that in the daytime. The agreement extent for the high-probability SUHI FPs between the two models is 74.91%, while 14.67% of the logistic model were classified to low-probability by the Gaussian surface model. The agreement extent for the low-probability SUHI FPs between the two models is 80.06%, while 6.94% of the logistic model were classified to low-probability by the Gaussian surface model.

We also analyzed the spatial differences of SUHI FPs between the two methods by overlaying different probability levels of SUHI FPs with land use maps (Tables S3–S6 in Supplementary Information). The major differences exist in cropland and impervious surface in the daytime, while in cropland, water bodies, and impervious surface at night. During the day, 11.46%, 24.55%, and 63.99% of the impervious surface were identified as low-probability, medium-probability, and high-probability SUHI FPs estimated by the logistic algorithm, respectively. However, the proportion of the impervious surface with high-probability SUHI FPs estimated by the Gaussian surface model is 51.09%, which is lower than that by the logistic model. During the day, 39.57%, 33.26%, and 27.17% of cropland were identified as low-probability, medium-probability, and high-probability SUHI FPs estimated by the logistic model. However, the proportion of the cropland being low-probability SUHI FPs estimated by the Gaussian surface model is 59.96%, which is also higher than that by the logistic model. In fact, many cropland pixels have been proven to be a heat source for the city [35]. In the nighttime, the proportions of impervious surface with high-probability SUHI FPs estimated by the Gaussian surface model and the logistic algorithm are 55.34% and 40.67%, respectively. Thus, the proportion of high-probability SUHI FPs estimated by the logistic algorithm in cropland was smaller than that estimated by the Gaussian surface model. However, the proportion of water bodies with high-probability SUHI FPs estimated by the logistic model reaches 59.62% at night, which is far higher than that by the Gaussian surface model. This is also consistent with the argument that the role of water could be a heat source for cities at night [35].

3.3. Temporal Changes in UHII

We first calculated the threshold LST value of the reference rural region for each 8-day MODIS LST product using Equation (4). Then, the UHII of each 8-day MODIS LST product was calculated using Equation (1), and the summer average maximum UHII value was calculated for all 8-day MODIS LST products for each year.

The average daytime maximum UHII was 4.27 ± 0.52 °C in summer, and the average nighttime maximum UHII was 3.04 ± 0.23 °C (Figure 9). We further used the Mann–Kendall Test and homogeneity tests to determine whether the day and nighttime series of interannual maximum UHII between 2004 and 2018 has a monotonic upward or downward trend. Results show that both the day and night time series of interannual maximum UHII have an upward trend (Sen’s slopes are 0.5048 and 0.0286 for the daytime and the nighttime, respectively). However, the upward trend for the daytime did not pass the 90% significance test, but the slight upward trend for the nighttime passed the 99% significance test. Radiation from the sun is the most important driver of climates near the ground [13]. In cities, urbanization, human-made structures and impervious surfaces replace vegetation, greenery, and water, altering radiation exchanges. The multi-faceted urban ‘surface’ that emits and reflects to itself can generate myriad distinct radiation budgets. These complex urban surfaces also absorb and store large amounts of solar shortwave radiation, while the reduction of sky-view factors decreases the loss of longwave radiation. In addition, cloud affects the transmission and emission of incoming radiation and increases the interception of outgoing radiation. Finally, air pollution may also affect the process of radiation exchanges [13].

To further understand the temporal changes of UHII, we calculated the summer average LSTs of high-probability SUHI FPs and the reference rural region between 2004 and 2018 (Figure 10).

There were increasing trends for the average LSTs of high-probability SUHI FPs and the reference rural region both in the daytime and in the nighttime. The extent of increase in the LST of the reference rural region was greater than that of the average temperature of SUHI FPs with the probability ranges of 0.8–0.9 and 0.9–1.0. The difference between day and night exists in terms of the level of significance of the increasing trend. The significant level for the daytime did not pass 90% significance test, but the significant level for the nighttime passed 90% significance test. These results explain why the UHII did not change significantly with the time between 2004 and 2018.

3.4. Temporal Changes in UHI Capacity

Compared to UHII, UHI capacities revealed more distinct temporal trends for different levels of UHI FPs. Generally, there were increasing trends for UHI capacities both in the daytime and in the nighttime between 2004 and 2018 (Figure 11). In particular, the UHI capacity of the SUHI FP with the highest probability range (0.9–1.0) has a more significantly increasing trend than the other two probability ranges. For example, in the daytime, the slope of the increasing trend for highest probability range is more than 15–27 times of those for the other two probability ranges, while in the nighttime, the growth rates are similar among the three high-probability UHI FPs.

4. Discussion

4.1. Elimination of Artificial Bias in Background Temperatures

With rapid economic development and population flow [64,65], many places in the world are experiencing rapid urbanization [66], catalyzing the formation of the urban heat island (UHI) [3]. To capture the UHI FP is of critical importance in understanding the process of urbanization and promoting sustainable urban development. In this study, we developed a new method for determining the SUHI FP. SUHI parameters have been calculated in previous studies using a variety of mathematical approaches, including the exponential decay and Gaussian function methods [17,44,45,46,47]. These approaches are generally limited by the method used of identifying the reference rural region, which uses the urban built-up area as a core to develop buffers outward toward rural regions [17]. These methods tend to use the shape of the city as the default shape of the SUHI FP. However, there is no consistent spatial relationship between these two shapes due to the impacts of topography and landscape heterogeneity on the LST. Other methods typically choose specific points or land use types as reference rural regions to calculate the UHII [42]. These methods ignore variation in LST attenuation in different directions. The UHII can change due to randomness and contingency among the reference objects. Furthermore, the reference rural regions can be subject to land use and land cover change through time [67], giving rise to further uncertainties. Therefore, determining the LST threshold in the reference rural region in each direction is key to determining the SUHI FP for a metropolitan region. In this study, we applied the polar coordinate system to determine different LST thresholds for each direction, which no longer needs to consider the shape of the city in the model.

4.2. Enhanced Reliability of Remote Sensing Data and Surface Attribute

Previous methods have had drawbacks in determining the spatiotemporal characteristics of SUHI effects using remote sensing data. Differences among the SUHI characteristics determined in this study, including maximum and reference LST values and SUHI intensity, on different dates within the same season also demonstrate the instability of remote sensing data. Therefore, if only single-phase remote sensing images or seasonal average LST data are used to study the UHII, it would lead to the uncertainty of quantifying SUHI FP and UHII in the study [35]. Besides the application of remote sensing data, we usually ignored much other information from the surface itself including walls, facets, and windows. Since the surface temperature is extremely sensitive to different properties (geometric, radiative, thermal, moisture and aerodynamic, etc.) that can change over time (day and night, season, and year) or under certain weather conditions, the SUHI is not static in terms of time or space. Considering these uncertainties resulting from remote sensing data and surface properties, we do not propose an absolute SUHI FP, but rather the application of pixel statistical probabilities to determine SUHI FPs.

4.3. Relationship Between SUHI FP and UHII

The average summer maximum UHII was 4.27 ± 0.52 °C and 3.04 ± 0.23 °C in the daytime and nighttime between 2004 and 2018, respectively. The Mann–Kendall Test confirmed that the upward trend (Sen’s slopes is 0.5048) did not pass the 90% significance test during the daytime. Although the upward trend did pass the 99% significance test during the nighttime, the upward trend with the Sen’s slopes value of 0.0286 was very slight. However, although there was no significant increase in maximum UHII, why was the SUHI FP getting larger between 2004 and 2018? The maximum UHII results were determined from the temperature difference between the maximum temperature in urban area and the threshold temperature in reference rural regions since these are essentially areas unaffected by the urban area. The reference rural region calculated by the logistic model means that the region should remain a long-term stable temperature background field for the city. This means that the location of the area affected by the city may change, but the temperature threshold for determining whether it is affected by the city does not significantly change. According to the algorithm of deriving the UHI FP, the UHI FP and the reference rural regions are all dynamically changing. Theoretically, the UHII of each pixel in the UHI FP can be calculated, but such dynamic references make the calculated UHII less meaningful than the UHI FP itself as well as the UHI capacity. Therefore, it is more important to use the logistic model to quantitatively calculate the temperature threshold of the reference rural region to accurately delineate the boundaries of the UHI FP. Both the UHI FP and the UHI capacity are increasing resulting from these areas transforming from non-SUHI FP areas to SUHI FP that are mainly derived from urbanization, i.e., natural surfaces are being replaced by urban structures (buildings, impervious surfaces).

5. Conclusions

We developed a logistic method to determine the spatiotemporal evolution of SUHI FPs and calculate the UHII in the Beijing metropolitan region during the summers of 2004–2018 using MODIS 8-day L3 LST global products. The logistic model well represented the trend of LST attenuation from the urban center to the surrounding rural regions. The key transition point, calculated using the rate of change in the curvature of the fitted logistic model, represented the LST threshold in reference rural regions, and was used to further determine the SUHI FP in the study area. To avoid the influence of remote sensing data instability, we calculated the probability of each pixel as part of the SUHI FP. The proposed method is relatively advantageous in that it minimizes the effects of the city’s shape, topography, or landscape heterogeneity. Our results also show that low-probability SUHI FPs had the largest area, followed by medium-probability and high-probability FPs. However, high-probability SUHI FPs increased significantly during the study period. The area of daytime high-probability SUHI FPs was larger than that in the nighttime. The daytime SUHI maximum intensity was higher than that during the night. However, the trend fluctuated, with no regular increase or decrease during either period.

If the rays can be used to fully represent the buffer rings from the urban center to the surrounding rural region in all directions in future studies, then the accuracy of SUHI FP identification may be further improved. Improving the accuracy of remote sensing data remains another issue for future study.

This research can be extended to different cities to automatically identify their SUHI FP and calculate UHII. More importantly, this algorithm can accurately draw the boundaries of SUHI FP for different UHII (1.5 ℃ or 2.0 ℃, etc.). It is of great significance to study urban climate and its mitigation measures in the context of global change. In addition, as the urban heat island footprint boundary is confirmed, the results can be applied in the field of urban planning. Urban planning aims at optimizing the allocation of land use in the urban region by scientific designing for new urban development land and ecological protection zone [68,69], which would directly affect the urban climate. Therefore, the scope of the SUHI FP provides a spatial decision reference for urban planning.