This study selected the Wuhan Urban Agglomeration, an important strategic place in Central China, and constructed three major elements of the TSPurban based on the MSPA model and circuit theory "heat source, resistance surface, and corridor," and used "reverse thinking" to explore a variety of ways to optimize the TSPurban. The main conclusions were as follows: (1) regarding the heat source, heat island patches were extracted according to the MSPA model, and the patches in the top 30 important heat island source areas were evaluated using the connectivity analysis method. It was found that the spatial form of the heat island was consistent with the urban spatial distribution; (2) regarding the resistance surface, the resistance surface of heat island propagation and diffusion was constructed by combining the land cover resistance factors, and the resistance surface was evaluated by using the standard deviation ellipse and spatial autocorrelation analysis methods. It was found that the resistance surface was distributed in an "east–west" spatial pattern and that the gathering situation was characterized by "cold in the west and hot in the east"; (3) regarding the corridor, we identified the link (corridor) based on the circuit theory and used the Linkage Mapper software to connect the lowest-cost path of the SUHI "source," thus generating 56 heat island corridors with a total length of 1004.71 km. At the same time, Circuitscape software was used to identify key nodes to effectively control key points, and the ecological regulation function for alleviating the urban heat island effect was more prominent. In terms of optimizing the TSPurban, it is necessary to weaken the interconnection between heat island sources at the spatial level by "breaking" the links and key nodes in the TSPurban network. Ultimately, the focus is on increasing the fragmentation of heat island patches and their resistance to diffusion and reducing the connections between heat islands while weakening the heat source patches themselves. This study utilizes reverse thinking to effectively mitigate the urban thermal environment.

**Author Contributions:** C.H.: Writing—original draft, Data curation. H.L.: Data curation, Software, Supervision. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Natural Science Foundation of China under grant 41807319 and the Fundamental Research Funds for the Central Universities (YCJJ202204012).

**Data Availability Statement:** The original contributions presented in the study are included in the article, and further inquiries can be directed to the corresponding author.

**Conflicts of Interest:** The authors declare that the research was conducted without commercial or financial relationships that could be construed as a potential conflict of interest.

### **References**

### *Article* **Spatial Influence of Multifaceted Environmental States on Habitat Quality: A Case Study of the Three Largest Chinese Urban Agglomerations**

**Fei Liu 1,\*, Yuji Murayama 2 and Yoshifumi Masago 1**

**Abstract:** Habitat structure and quality in the urban agglomeration (UA) are subject to multiple threats and pressures due to ongoing anthropogenic activities and call for comprehensively effective solutions. Many approaches, including cartographic comparison, correlation analysis, the local entropy model, and GeoDetector, were jointly used to clarify the interplay between habitat quality and multiple environmental issues. In response to the overlapped risks of diverse environmental systems, this study presented an integrated research framework to evaluate the spatial influences of multifaceted environmental situations on habitat quality. We conducted the case study in the three largest Chinese UAs: Beijing–Tianjin–Hebei (BTH), Greater Bay Area (GBA), and Yangtze River Delta (YRD). The evaluation results show that the three UAs shared similarities and differences in relationship/impact types and their strengths. In 2015, most of the three UAs' landscapes delivered low–medium magnitudes of habitat quality (score <0.7) and emerged with unevenly severe consequences over space across different environmental aspects, highlighting the importance of maintaining habitat safety. Overall, habitat quality scores were synergistic with NDVI, but antagonistic to surface heat island intensity (SHII), PM2.5 concentrations, and residential support. However, locally structured relationships exhibited geographical complexity and heterogeneity between habitat quality and environmental systems. Regarding GeoDetector evaluation, PM2.5 concentrations in BTH, SHII in GBA, and NDVI in YRD played a dominant role in single-factor and interaction analysis. More importantly, the synergistic effect of various environmental issues on habitats was manifested as mutually enhanced rather than independent or weakened interactive effects, implying the aggravation of compound effects and the necessity of prioritization schemes. This study could provide beneficial insights into the interconnections between habitats' sustainability and multifaceted environmental situations in UAs.

**Keywords:** air quality; GeoDetector; habitat quality; thermal environment; urban agglomeration

### **1. Introduction**

Ongoing socioeconomic development has given rise to multifaceted problems and challenges for habitats. The Intergovernmental Panel on Climate Change (IPCC) highlighted the necessity of understanding the interrelationships and linkages between various environmental issues in the Sixth Assessment Report [1]. The natural ecosystem and human system are intertwined, and alterations in one component may stimulate others [1]. Thus, sustainable development strategies require a new understanding of compound environmental issues and overlapping consequences [2].

The United Nations declared that over half of the current world's dwellers reside in urban settlements in the World Urbanization Prospects [3]. This proportion will rise to two-thirds by the 2050s [3]. Accordingly, improving human well-being increasingly relies on upgrading the situations and services in human-dominated urban areas [4,5].

**Citation:** Liu, F.; Murayama, Y.; Masago, Y. Spatial Influence of Multifaceted Environmental States on Habitat Quality: A Case Study of the Three Largest Chinese Urban Agglomerations. *Remote Sens.* **2023**, *15*, 921. https://doi.org/10.3390/ rs15040921

Academic Editor: Riccardo Buccolieri

Received: 4 December 2022 Revised: 3 February 2023 Accepted: 4 February 2023 Published: 7 February 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Urban agglomerations (UAs) are a newly evolved manifestation of spatially clustered regions with a common interest and fate [6]. Regional integration development and drastic urbanization have propelled the formation of UAs [7]. The UA system is a complex ensemble comprising multiple cities and their extended outskirts with large-scale urban fabrics and artificial landscapes. These artificial landscapes provide adequate space to guarantee social-economic boom and human well-being [6]. In parallel, natural habitats provide a majority of essential ecosystem services to sustain the ecological security of the entire UA region. However, with drastic economic and population growth, UA regions have manifested vigorous trade-offs between ecosystem services provision and human social development, giving rise to enormous ecological deficits, especially in the core megacities and highly populated areas [4,8]. As such, the habitats of UAs are diverse but fragile, and might be subjected to a multitude of risk factors from natural and human interference [4].

Existing studies on UAs are typically confined to the individual environmental problem, formation mechanism, and mitigation strategy, including the thermal [9], air [10], housing [11], landscape [12], and resource [13] spheres. A vast quantity of studies have identified, for example, the heat island phenomenon [14,15], air pollution [16,17], and urban containment [18]. Abundant solutions have continuously contributed to a future sustainable world, such as nature-based solutions [19] and ecological intensification interventions [20]. Against the backdrops of UAs, implementing measures that safeguard environmental health and human well-being across multiple aspects is warranted [21,22]. However, few studies have been concerned with the fate of compound environmental risks and their interactive relationships. The environmental crisis in the UAs may be magnified due to the overlap of multiple issues, their influences on habitats may be aggravated, and comprehensive carrying capacity might have further deteriorated [1,23]. Considering the growing demands for comprehensive issues and compound risks in the UAs, it is imperative to understand the associations between diffident environmental problems and seek co-benefits.

Unraveling the relationships between habitats and environmental issues is the foundation of effective, comprehensive policies. Yet, one of the key challenges is obtaining insights into the complex interplay of multiple environmental situations [1,4]. The impact of human interference on different environment systems varies broadly from slight, modest, to severe, but may share similar driving factors [8,24–26]. Additionally, the measurement of the geographical association of multiple variables relies on their locations, distances, neighboring settings, and scales [24,27,28]. However, conventional methods, such as cartographic comparison, correlation, and regression analysis, generally only mirror simple relationships [29]. Exponential, quadratic, and interactive relationships are difficult to represent using these conventional methods. The lack of knowledge about complex relationships has obstructed the achievement of multifaceted environmental evaluation [29]. In light of this, this study presents a framework that integrates multi-method and multi-source data to evaluate the associations between multifaceted environmental states and habitat quality in different UA settings. Our proposed framework incorporates cartographic comparison, correlation analysis, the local entropy model, and GeoDetector to conduct a spatially multi-perspective relationship profiling. The local entropy model is a novel method for examining the locally varying relationships over space, which can capture linear and polynomial complex associations [29]. GeoDetector is a powerful tool for identifying spatial stratified heterogeneities and interactive effects [30,31].

We designated the UAs of Beijing–Tianjin–Hebei (BTH), Guangdong–Hong Kong– Macao Greater Bay Area (GBA), and Yangtze River Delta (YRD) as the study regions. In recent decades, the Chinese government has been committed to ecological civilization construction and sustainability management in these three UAs based on a series of relevant initiatives and policies [32–34], which have been in response to habitat degradation from dramatic socioeconomic growth [13,35–38]. Moreover, many scientific investigations have contributed to sustainable development schemes for the three UAs [37,39–44]. The thermal environment, air quality, living space, and green landscapes have been the focal topics

among environmental monitoring and assessments, which are also deeply discussed in this study. Thus, focusing on the three UAs of BTH, GBA, and YRD in 2015, this study consists of four goals, as follows: (1) Quantification of habitat quality; (2) Evaluation of diverse environmental situations from the thermal, air, living, and landscape spheres; (3) Identification of the interconnection and interplay between habitat quality and multidimensional environmental situations; and (4) Elucidation of a potential management pathway for multi-dimensional sustainability and security for the UAs.

### **2. Materials and Methods**

### *2.1. Study Areas*

This study focused on the three largest Chinese UAs: BTH, GBA, and YRD (Figure 1). The study regions were determined due to (1) their rapid socioeconomic development, (2) plentiful nature capitals, (3) critical geographical location, (4) urgent sustainability needs, and (5) superior policy support. BTH, GBA, and YRD are nationally and globally representative UA regions. These UA regions are not only primary engines and important pivots for China's socioeconomic development, but also the significant guarantee and forefront of China's ecological security. BTH UA (113◦27 E to 119◦50 E, 35◦03 N to 42◦40 N) is situated in Northeast China, encompassing 13 cities with 2.2% (2.2 × 107 ha) of the national territorial area. GBA UA (111◦27 E to 115◦42 E, 21◦57 N to 24◦40 N) has a scope of 9 cities and 2 special administrative regions in South China, accounting for 0.6% (5.6 × 106 ha) of the national territory. YRD UA (111◦27 E to 115◦42 E, 21◦57 N to 24◦40 N) is composed of 16 cities on the eastern coast of China, covering 2.2% of China's territory (2.1 × 107 ha).

**Figure 1.** Study areas: (**a**) geographic location of the study areas in China; (**b**) Beijing–Tianjin–Hebei (BTH); (**c**) Guangdong–Hong Kong–Macao Greater Bay Area (GBA); and (**d**) Yangtze River Delta (YRD).

At the end of 2015, approximately 8%, 5%, and 11% of the national population resided in BTH, GBA, and YRD, respectively [45]. BTH, GBA, and YRD had 10.2%, 14%, and 18.5% of the respective national totals in the gross domestic product (GDP) [45]. These three UAs emerge as compound habitat systems of multiple urban cores and their peripheral countrysides under different climate schemes characterized by abundant land details. BTH, GBA, and YRD are configured by a variety of ecosystems, landscapes, and landforms.  Significant alterations and transformations in the environmental landscapes and socioeconomic scales of these UAs in recent decades [38,42,46,47] have given rise to habitat degradation, ecological fragility, and climate crisis [48–51], such as intensified heat island effects [52] and deteriorated air pollution [16]. Thus, national development initiatives also highlight the urgency and necessity of sustainability-related, ecological-oriented strategies and projects in BTH, GBA, and YRD [32–34].

### *2.2. Data Preparation and Treatment*

Given the data availability, precision, and consistency of the three UAs, this study set 2015 as the year for the research. The year 2015 was an extraordinary time for economic growth and social boom in BTH, GBA, YRD, and even the entire territory of China due to a number of important development policies, including the 13th Five-Year Plan, Belt and Road initiatives, and regional collaboration plans [32–34,53].

The land cover/use maps of the BTH, GBA, and YRD in 2015 were cited and pruned from Climate Change Initiative (CCI) global land cover (LC) products, European Space Agency (ESA) (https://www.esa-landcover-cci.org/ (accessed on 30 November 2022)). ESA CCI-LC product is a 300 m spatial resolution of a time-series global surface coverage dataset generated based on multi-source satellite archives [54]. The overall weighted-area accuracy of the ESA CCI-LC map in 2015 was validated to be 71.1% [54]. Thus, the ESA CCI-LC map can deliver abundant land details and flexibly serve diverse land-oriented applications. Referring to the official classification scheme [54] and IUCN habitat-related suggestions [55,56], we rearranged the land cover/use thematic legend into eight categories: agricultural habitat (AH), forest habitat (FH), grassland habitat (GH), wetland habitat (WH), shrubland habitat (SH), urban fabric (UF), vacant land (VL), and water (WA).

The annual average daytime and nighttime land surface temperatures (LSTs) of the BTH, GBA, and YRD in 2015 were composited and preprocessed based on the 1 km resolution of Moderate Resolution Imaging Spectroradiometer (MODIS) daily/night land surface temperature/emissivity products (MOD11A1 and MYD11A1, https://modis.gsfc.nasa. gov/data/dataprod/mod11.php (accessed on 30 November 2022)) using the Google Earth Engine (GEE). The 2015 yearly average normalized difference vegetation index (NDVI) distributions in the BTH, GBA, and YRD were also GEE-derived and composed using MODIS 1 km 16-Day vegetation indices (MOD13A2, https://modis.gsfc.nasa.gov/data/ dataprod/mod13.php (accessed on 30 November 2022)). The 1 km resolution spatial coverages of the 2015 average PM2.5 (particulate matter in the air with aerodynamic diameters smaller than 2.5 μm) concentrations in the BTH, GBA, and YRD were sourced from China High Air Pollutants (CHAP) datasets [57,58] (https://weijing-rs.github.io/product.html (accessed on 30 November 2022)). The China High PM2.5 datasets were developed by coupling in site PM2.5 data, MODIS aerosol optical depth (AOD) products, and various auxiliary data based on the space-time extra-trees model [57]. Various air-relevant studies have used the ChinaHighPM2.5 datasets due to their wide spatiotemporal range and high cross-validation accuracy (R2 = 0.86–0.90) [57,58]. Spatially gridded population density data (1 km of spatial resolution) of the BTH, GBA, and YRD in 2015 were extracted from the high-profile WorldPop population products. WorldPop (https://www.worldpop.org/ (accessed on 30 November 2022)), developed by the University of Southampton, has been dedicated to supporting the open utilization of various spatial demographic studies [59]. China-relevant spatial data (e.g., administrative boundary) were obtained from the Resource and Environmental Science and Data Center, Chinese Academy of Sciences (https://www.resdc.cn/Default.aspx (accessed on 30 November 2022)).

### *2.3. Research Framework*

This study aims to evaluate the interconnection and spatial heterogeneity between habitat quality and multiple environmental situations from multiple perspectives. Thus, the cascade diagram of our proposed framework is shown in Figure 2. The overall research framework is mainly composed of four modules: (I) data preparation and preprocessing; (II) correlation analysis; (III) local bivariate analysis; and (IV) GeoDetector analysis.

**Figure 2.** Overall research framework. Note: InVEST logo image sources from the website (https: //naturalcapitalproject.stanford.edu/software/invest (accessed on 30 November 2022)).

### *2.4. Quantification of Habitat Quality*

Habitat quality is geographically heterogeneous, which depends on habitat characteristics, suitability, and the presence and relative magnitudes of threat factors [60]. In light of this, the spatial patterns and magnitudes could be evaluated by the joint use of land cover/use and threats data [61,62]. The Integrated Valuation of Ecosystem Services and Tradeoffs–Habitat Quality (InVEST–HQ) model was developed by Natural Capital Projects in conjunction with Stanford University, the Nature Conservancy, and the World Wide Fund for Nature (https://naturalcapitalproject.stanford.edu/software/invest (accessed on 30 November 2022)). Given the low data demand, superior spatial visualization, and easy operation, the InVEST–HQ model is the dominant tool for habitat quality quantification. It has been widely applied to ecological studies across different regions and spatial scales [60]. Thus, we employed the InVEST–HQ model to estimate the habitat quality in the three UAs of BTH, GBA, and YRD. Habitat quality can be calculated using the following equations:

$$Q\_{xj} = H\_{\dot{f}} (1 - (\frac{D\_{x\dot{j}}^z}{D\_{x\dot{j}}^z + k^z})) \tag{1}$$

where *Qxj* represents the habitat quality of grid *x* in land cover/use type *j*; *Hj* refers to the habitat suitability of land cover/use type *j*; *Dxj* reflects the total threat level in grid cell *x* with land cover/use type *j*, which is judged by the types, intensities, and proximities of neighboring threats [62]; *k* is the half-saturation constant, which is equivalent to the half of the crest value of *Dxj*; and *z* is default at 2.5, which is a scaling parameter for mirroring spatial differentiation [60].

$$D\_{xj} = \sum\_{r=1}^{R} \sum\_{y=1}^{Y\_r} (\frac{w\_r}{\sum\_{r=1}^{R} w\_r}) r\_y i\_{rxy} \beta\_x S\_{jr} \tag{2}$$

where *R* denotes the number of threats; *y* indexes all grid cells on the raster layer of threat *r*; *Yr* indicates the set of grid cells occupied by the threat on the raster layer; *wr* and *ry* reflect the weight of the threat and the interference intensity of the threat for grid cell y, respectively; *irxy* is the habitat interference level of grid cell *x* from threat *r* on the grid cell y; *βx* refers to the reachability level (the anti-disturbance ability) of grid cell *x*; and *Sjr* represents the relative sensitivity of land cover/use type *j* to threat factor *r*. The habitat quality scores range from 0 to 1, which embodies serious to ideal situations for habitat quality [60]. In this paper, *irxy* was considered as follows:

$$d\_{rxy} = 1 - \left(\frac{d\_{xy}}{d\_{r\\_max}}\right) if\lim\_{} \text{linear attenuation} \tag{3}$$

$$d\_{\rm rxy} = \exp\left(-\left(\frac{2.99}{d\_{\rm r\,max}}\right)d\_{\rm xy}\right) \text{if } \exponentially \text{ decay} \tag{4}$$

where *dxy* is the linear distance between grid cells *x* and *y*, and *dr max* is the maximum effective influence distance of threat *r*. More specific information is provided in the InVEST user guide [60] and the relevant literature [35,49,63–66].

We primarily prepared three categories of documentation as input parameters of the InVEST–HQ model: (1) raster layers of land cover/use, (2) threat data, sources, and their accessibility, and (3) suitable habitat types and corresponding threat sensitivities [60]. Threat sources in this study included cultivated lands, bare lands, and built-up lands due to data availability and consistency. The raster layers of threat sources were derived in ArcGIS Pro Version 3.0.1 (https://www.esrij.com/products/arcgis-pro/ (accessed on 30 November 2022)). The setting of threat factors and sensitivity parameters were based on previous literature in the study regions, expert knowledge, and the official user guide manual [35,49,60,63–66], as shown in the Supplementary Materials.

### *2.5. Evaluation of Multifaceted Environmental States*

Land surface temperature (LST) is an essential key variable in surface energy balance, heat fluxes, and energy exchanges over different land cover types, playing a crucial role in geo-environmental and human health [28]. Thermal infrared (TIR) data from satellite platforms and retrieved LST are reliable barometers for diagnosing thermal circumstances in various fields [28,67]. Numerous studies have found different responses to daytime and nighttime LSTs and reported their different influences on environmental situations and human behaviors [28,67]. Thus, we utilized both daytime and nighttime LSTs. Instead of absolute LST observations, the intensity of the surface heat island (SHII) was estimated for evaluating the states of thermal environments in BTH, GBA, and YRD to minimize the effects and uncertainty of the climate background. The SHII was generated using the average daytime and nighttime surface heat island intensities, computed as Equations (5) and (6). The day and night SHII were, respectively, derived from the daytime and nighttime LSTs, which subtracts the corresponding mean LST value of non-UF category pixels from the LST value of each pixel in the study region [68]:

$$SHII = average \left( SHII\_{day}, SHII\_{night} \right) \tag{5}$$

$$SHII\_{day\text{ }or\text{ }night} = LST\_{turbnn} - LST\_{runl} = LST\_i - LST\_{non-lIF} \tag{6}$$

where *SHIIday* and *SHIInight* are the surface heat island intensities for day and night, respectively. *LSTurban* is the LST value of urban pixels, while *LSTrural* indicates the LST value of rural pixels. *LSTi* denotes the LST value of pixel *i* within the individual study regions of BTH, GBA, and YRD. The non-UF category pixels are substitutions for rural pixels, encompassing the pixels from various non-UF land cover/use categories. Thus, *LSTnon*−*UF* indicates the average LST values of non-UF pixels in the three UAs.

In this study, the yearly average concentration of PM2.5 was captured to reflect the air environmental situation and pollution in the three UAs. The values of the yearly average NDVI were used to mirror the environmental comfort and biological health degrees based on green spaces and vegetation situations.

Residential support (RS) refers to the spaces or services arranged for the necessities of human living and life, taking into consideration the land cover/use weights and population density [69]. BTH, GBA, and YRD are the most populated regions in China. Thus, residential support is an essential component of UA's environmental issues. The magnitude of RS in this study was evaluated based on the population distribution and built-up proportion, which is computed using the following equation [69]:

$$RS = Pop\_{std} \times P\_{BL} \tag{7}$$

where *Popstd* indicates standardized population density. *PBU* denotes the built-up proportion derived from the calculation of land cover/use.

### *2.6. Relationship Profiling*

### 2.6.1. Correlation and Local Bivariate Analysis

Given the data-processing capacity and spatial autocorrelation, we generated a set of 3000 m × 3000 m grid cells as the appropriate analytical scale for the spatial analysis. Thus, this study rescaled the habitat quality and all environmental indicators into 3000 m spatial fishnets for evaluating the relationships and interactions between habitat quality and multifaceted environmental situations. All the spatial operations and modeling in this study were completed in the ArcGIS Pro Version 3.0.1 platform. Initially, we utilized the corrplot R package (under R v4.1.0 statistical software, https://www.r-project.org/ (accessed on 30 November 2022)) to visualize the correlation matrix for disclosing the global association between various indicators across the entire UA region. Next, we adopted the local entropy model to explore the local connection between habitat quality and various environmental situations.

The local spatial associations between habitat quality and environmental situations were judged based on a local entropy model. The local entropy model is a nonparametric approach proposed by Guo [29] for multivariate data. It integrates local spatial analysis, Rényi entropy calculation, permutation-based distribution estimation, and a set of statistical tests, allowing the interactive detection of significant local multivariate associations over space [29]. Entropy is a gauge for quantifying uncertainties and randomness, such as Rényi entropy and Shannon entropy [29,70,71]. Rényi and Shannon entropies are computed as the following equations:

$$H\_{\lambda} = \frac{1}{1 - \lambda} \log \left( \int f(x)^{\lambda} dx \right) \quad \lambda \ge 0, \lambda \ne 1 \tag{8}$$

$$H = -\int f(\mathbf{x}) \log f(\mathbf{x}) d\mathbf{x} \tag{9}$$

where *Hλ* indicates Rényi entropy, *x* denotes a multi-dimensional vector in the data space, *f*(*x*) refers to the probability density function, and *λ* (≥ 0) represents the order of the *Hλ*. When *λ* is closing to 1, *Hλ* converges to the Shannon entropy (*H*) in the data space [29]. However, the unknown probability density function is the principal challenge when estimating the entropy for exploratory data analysis [29,70,71]. The power-weighted minimum spanning tree (MST) method is an effective alternative for entropy estimation in multivariate data analysis [29]. Accordingly, the local entropy model can be conducted by the following analytical procedures: (1) multivariate Rényi entropy estimation and observation for the local neighborhood, (2) construction for empirical distributions of local joint entropy using permutations based on the null hypothesis and MST method, (3) test for statistically significant local relationships using permutations, and (4) examination and classification of the local relationships [29,70,71]. Notably, using the appropriate parameter values for the size of the local neighborhood and the edge power of the MST is vital in implementing this local entropy model. Specifically, the scaling factor can modulate the sensitivity to subtle relationships. The proper parameter of the neighboring size can improve the likelihood and efficiency of detecting significant relationships and patterns. Additionally, optimizing the number of permutations can permit a neutralization between precision and processing time [29].

The local entropy model proposed by Guo [29] did not assume a prior relationship form, which can estimate multivariate entropy distribution without a probability density function and implement local spatial analysis without a regression model [29]. Thus, there is a growing application of local entropy maps in geography-related studies [29,70–72]. We designed the research steps of local bivariate analysis based on Guo's approach, as shown in Figure 3. The average habitat quality score of each analytical grid was extracted as the dependent variable. The respective average values of SHII, PM2.5 concentration, NDVI, and RS for each analytical grid were estimated as explanatory variable parameters. Further, we conducted a series of local bivariate relationship analyses to investigate the local association between habitat quality and each explanatory variable. We specified a 95% confidence level of the hypothesis test for significant relationships. Six categories of relationships can be identified based on the local entropy model: not statistically significant, positive linear, negative linear, concave, convex, and undefined complex.

**Figure 3.** Research flowchart of local bivariate analysis. Note: The flowchart of local bivariate analysis was designed and revised based on Guo's study [29], consisting of three steps: (I) local data

preparation; (II) permutation-based entropy distribution estimation and testing; and (III) local spatial

### 2.6.2. Geographical Detector Model

relationship mapping.

The Geographical Detector Model (GeoDetector) is a novel method for capturing stratified spatial heterogeneity and related driving forces, consisting of risk, factor, ecology, and interaction detectors modules [30,31]. GeoDetector can quantitatively analyze the association and similarity of the spatial distributions between a geostatistical variable and a series of stratified independent variables [31]. The merits of GeoDetector are identifying interaction influences between multiple factors and comparing their differences [31].

The GeoDetector model assumes that the study area is composed of multiple region blocks. Spatial heterogeneity is detected if the total variance of the regions is greater than the sum of the variances of the subregions [31]. In contrast, there is a statistical association between the two variables when the spatial patterns of the two variables tend to be identical [31]. Due to no presuppositions and constraints, the GeoDetector analysis is universal and transplantable [31]. Thus, the GeoDetector model has been used in various research fields, such as the geoenvironmental and sustainability sciences [14,73,74].

Given the consistency of the three UAs, the four numerical indicators—RS, SHII, PM2.5 concentrations, and NDVI—were stratified and categorized using equal intervals classification in conjunction with Jenks natural breaks optimization methods. These four stratified indicators served as the explanatory variables (X) in the GeoDetector analysis. The estimated score of habitat quality was still taken as the dependent variable (Y). Among the four modules of GeoDetector, the risk detector showed the responses of habitat quality on each stratified explanatory variable. The ecological detector was used to examine the significant impact difference of stratified explanatory variables on habitat quality. These two modules can determine the suitable range or type of explanatory variables. However, we were more concerned about the outputs of factor and interaction detectors in this study.

Based on the factor detector module, the magnitude of stratified spatial heterogeneity for habitat quality can be quantified using *q* statistical measures [31].

$$q = 1 - \frac{\sum\_{h=1}^{L} N\_h \sigma\_h^2}{N\sigma^2} = 1 - \frac{SSW}{SST} \tag{10}$$

where the dependent variable Y (habitat quality in this study) is divided into strata {*h*}, with *h* = 1, ... , *L*; *N* and *Nh* indicate the number of grids in the whole UA region and strata h, respectively; and *σ*2 and *σ*2 *h* denote the variances of the whole region and strata *h*, respectively. *SSW* is the sum of squares within, while *SST* is the sum of squares total. *q* stands for the explanatory power of a specific explanatory variable X (RS, SHII, PM2.5, and NDVI) for habitat quality, ranging from 0 to 1. The greater the *q* value, the stronger the influences of the explanatory variable X and the more apparent spatial differentiation of habitat quality [31]. Remarkably, habitat quality exhibits a random spatial distribution when the *q* value is 0, whereas habitat quality shows complete spatial differentiation when the *q* value equals 1 [31]. The *q* value satisfies the non-central *F* distribution, where *λ* is the non-central parameter and *Yh* is the mean of layer *h* [31]. -

$$F = \frac{N-L}{L-1} \begin{array}{c} q \\ \hline 1-q \end{array} \begin{array}{c} F \widehat{(L1\text{\textquotedbl{}} \text{\textquotedbl{}} NL \text{\textquotedbl{}})} \end{array} \tag{11}$$

$$\lambda = \frac{1}{\sigma^2} [\sum\_{h=1}^{L} \overline{Y\_h}^2 - \frac{1}{N} (\sum\_{h=1}^{L} \sqrt{N\_h Y\_h})^2] \tag{12}$$

According to the equations above, the *p*-value can be calculated as a statistical significance test for the *q* statistics. The interaction detector module can discern the interactive effect on the habitat quality between two environmental situations. Firstly, the impact (q values) of two stratified environmental situations on habitat quality was estimated as *q*(*Xi*) and *q Xj* , respectively. Next, we computed the *q* values of the interactive effect (*q Xi* ∩ *Xj* ). Ultimately, according to judging criteria (Table 1), the interaction type between two stratified environmental situations can be determined by comparing *q Xi* ∩ *Xj* with *q*(*Xi*) and *q Xj* . More detailed information on the GeoDetector model can be obtained from the official website (www.geodetector.cn (accessed on 30 November 2022)).

**Table 1.** Definition of the interaction types in the GeoDetector model [31].

### **3. Results**

### *3.1. Spatial Characteristics of Habitat Quality Based on Land Cover/Use Evaluation*

The thematic maps of land cover/use for three UAs were prerequisites for quantifying the habitat quality, as shown in Figure 4. The UF areas were highly urbanized areas dominated by large-scale built-up areas, highlighted in red. Natural habitat types are presented using different gradients of green and blue. AH referred to a semi-habitat system used for cultivation and farming activities, depicted in yellow.

Upon the completion of the InVEST–HQ analysis, the magnitudes of habitat quality for three UAs in 2015 were assessed, as shown in Figure 5. The habitat quality scores ranged from 0 to 1, representing the worst to perfect ecological situations. We assorted and mapped the habitat quality index into ten levels using an equal interval of 0.1, indicating the lowest (HL1: 0–0.1) to highest habitat quality (HL10: 0.9–1). The habitat qualities of the three UAs and their geographical variations differed by habitat types and suitability. Overall, the habitat situations of the three UAs in 2015 were not optimistic, and approximately one-third of the regions showed low magnitudes of habitat quality (HL1–HL3). The average habitat quality scores in BTH, GBA, and YRD were 0.3563, 0.3399, and 0.3336, respectively. The moderate-quality levels of habitats (HL4–HL6) made up the largest proportion of the three UA regions. However, the high habitat quality areas (HL7–HL10) added up to under 15% (Figure 5d). The spatial layouts of habitat quality in the three UAs were highly consistent with the configuration characteristics of land cover/use. The scores of non-habitat areas (UF and VL) were estimated to be zero. In contrast, FH and WA showed higher quality scores than other habitat types in the three UAs. However, the habitat quality score of the WA type in GBA was much lower than that of BTH and YRD. There was a large margin of habitat quality score for the GH type between BTH and the other two UAs. The average scores of AH, FH, GH, and WA in 2015 are shown in Figure 6, as these four habitat types dominated the three UAs (over 85% of area shares).

**Figure 4.** Thematic maps of land cover/use for three UAs: (**a**) BTH, (**b**) GBA, and (**c**) YRD.

**Figure 5.** Spatial patterns and distributions of habitat quality in the three UAs: (**a**) BTH, (**b**) GBA, (**c**) YRD, and (**d**) area summary.

**Figure 6.** The scores of habitat quality for the four dominant land cover/use.

Geographically, the areas with low habitat suitability in the three UAs in 2015 were chiefly located in human settlements, namely, urban and peri-urban fringes, indicating the adverse influences of anthropogenic activities on the regional habitats. The northwest part of BTH showed high-quality scores due to mountain ridges and massive forests, while the most serious habitat fragile areas were detected in Beijing and Tianjin. The severe habitat delicate areas in GBA were roughly configured as a triangle-shaped region, extending from Guangzhou–Foshan outwards to Macao, Zhuhai, Shenzhen, and Hong Kong. Only Zhaoqing and Huizhou exhibited favorable circumstances for habitats due to their abundant natural capital. In YRD, the southern part and the Taihu region were prominent due to their superior habitat suitability. Shanghai, the south of the Jiangsu Province, and the north of the Zhejiang Province suffered from a challenging habitat situation in 2015 due to the sparse vegetation availability.

### *3.2. Spatial Stratification of Multi-dimensional Environmental Situations*

Although the annual LST values showed a large difference in tendency and influence depending on the day and night, the spatial patterns of day and night heat island intensities in the three UAs were highly consistent. In light of this, this study utilized the mean value of yearly daytime and nighttime SHII to reflect the regional comprehensive thermal environment status in the three UAs. The spatial arrangements of daytime and nighttime thermal environments (both LST and SHII) in the three UAs are shown in the Supplementary Materials. This study divided the SHII index into five categories: ≤0, 0~1, 1~2, 2~3, >3 (unit: ◦C, Grades I–V) (Figure 7). The 2015 average values of daytime LST in BTH, GBA, and YRD were estimated to be 19.99 ◦C, 24.7 ◦C, and 20.83 ◦C, respectively, and those of respective nighttime LST were 21.84 ◦C, 26.79 ◦C, and 22.67 ◦C. The non-habitat areas and agricultural habitats dominated the areas with high SHII values. The area ratios of positive SHII regions were 56.36% in BTH, 57.12% in GBA, and 60.51% in YRD.

Significant spatial differentiations of yearly average PM2.5 concentrations were found in the three UAs. The spatial maps of the average PM2.5 concentrations in 2015 for the three UAs were stratified into five categories: ≤30, 30~40, 40~50, 50~60, >60 (unit: ug/m3, Grades I–V), as shown in Figure 8. The BTH regions showed the worst air condition among the three UAs, and about 43.78% of the regional territory was exposed to severe air pollution (Grade V). Most of the GBA regions (94.26%) in 2015 were clustered in Grade II and III PM2.5 concentrations. Only 13.8% of YRD regions did not suffer from severe air pollution.

This study weighed the UAs' living supporting capacity based on the product of built-up proportion and population density. Regions without living supporting capacity were assigned the value of zero. The RS value was equal to one, indicating perfect residential suitability. We defined five grades for the RS index: 0 (without RS capacity), 0~0.1, 0.1~0.2, 0.2~0.3, >0.3 (unitless, Grades I–V). There were multiple cores with extremely high living supporting capacity in the three UAs, as shown in Figure 9. Regions with residential suitability in BTH, GBA, and YRD were estimated to be 26.75%, 26.27%, and 30.1%, respectively.

**Figure 7.** Spatial stratification of the thermal environmental situations in the three UAs: (**a**) BTH, (**b**) GBA, and (**c**) YRD.

**Figure 8.** Spatial stratification of the air environmental situations in the three UAs: (**a**) BTH, (**b**) GBA, and (**c**) YRD.

**Figure 9.** Spatial stratification of the residential situations in the three UAs: (**a**) BTH, (**b**) GBA, and (**c**) YRD.

NDVI represents the aesthetic, recreational, and regulating functions of natural landscapes. We mapped the spatial patterns of the NDVI index with an interval of 0.2 for the three UAs (unitless, Grades I–V, Figure 10). In 2015, the average values of NDVI were projected to be 0.7263 in BTH, 0.7196 in GBA, and 0.6798 in YRD. Approximately 84.77% of BTH regions, 75.33% of GBA regions, and 72.30% of YRD regions had high values of NDVI, exceeding 0.6 in 2015.

**Figure 10.** Spatial stratification of the biological health situations in the three UAs: (**a**) BTH, (**b**) GBA, and (**c**) YRD.

### *3.3. Correlations of Environmental States*

We started with the correlation analysis between habitat quality and four environmental status indicators. In general, three UAs exhibited identical relation patterns in the entire regional scale, albeit with varying degrees. Habitat quality in the three UAs was negatively associated with the indicators of SHII (R2 ≥ 0.49, *p* < 0.001), PM2.5 concentrations (R2 ≥ 0.38, *p* < 0.001), and RS (R2 ≥ 0.26, *p* < 0.001), while positively related to NDVI (R2 ≥ 0.3, *p* < 0.001). The correlation coefficients between habitat quality and SHII were the highest in GBA (R2 = 0.73) and the lowest in BTH (R2 = 0.49). The association between habitat quality and NDVI degrees in GBA (R2 = 0.62) was the strongest among the three UAs. The most salient impact of PM2.5 concentrations (R2 = 0.44) on habitat was detected in BTH. YRD showed the highest coefficient between habitat quality and RS (R2 = 0.32). As for the four environmental status indicators, there were positive correlations among SHII, PM2.5 concentration, and RS, but NDVI negatively varied with other indicators (*p* < 0.001). In terms of UA, the strongest associations were SHII–PM2.5 concentration in BTH (R2 = 0.81) and SHII–NDVI in GBA (R2 = 0.74) and YRD (R2 = 0.64). In contrast, there were minor linkages between PM2.5 concentrations and RS in all the studied UAs (R2 ≤ 0.12). The relevant correlation coefficients are shown in Figure 11.

We identified six types of spatial structural relationships between habitat quality and four environmental status indicators based on the local entropy model. We assumed that there were local relationships between quality and individual environmental status indicators in the three UAs under the 95% confidence level. We set the 0.5 scaling factor, 30 neighborhoods, and 199 permutations to perform the hypothesis test for accommodating local spatial relationships.

**Figure 11.** Correlograms of multifaceted environmental situations in the three UAs: (**a**) BTH, (**b**) GBA, and (**c**) YRD. Note: The asterisk symbols ('\*\*\*') are representative of the statistical significance at the 0.1% level (*p* < 0.001).

The relationships were predicted and determined based on the adjusted R-square and the model's Modified Akaike Information Criterion (AICc). The classified local relationships in the three UAs are mapped in Figure 12, and the relevant categorical summary is shown in Table 2. We provided the entropy statistic and detailed information on the local bivariate analysis in the Supplementary Materials. The results of the local bivariate analysis displayed local structural relationships of various forms between habitat quality and environmental situations over space. Furthermore, the generated maps of local bivariate analysis differentiated the significance levels within a certain neighboring extent. Many regions showed statistically non-significant relationships in the local contexts of the three UAs, indicating a spatially random distribution pattern between variables presented in these regions. The spatially explicit local relationship maps provided an overview of the spatial diversities for significant local associations between habitat quality and environmental situations and navigated the analysis to focus on specific areas. The local bivariate analysis, thus, offered considerable interpretation clues. For example, about 4.86% of BTH regions showed a concave relationship between habitat quality and SHII, indicating that habitat quality in the corresponding BTH regions followed a concave curve as the SHII varied.

We found that the total area ratios of linear and quadratic models for the HQ–SHII relationship in BTH were 23.95% and 10.59%, respectively. This finding indicates that the polynomial and linear models can explain the 23.95% and 10.59% relationship between habitat quality and thermal environment in BTH, respectively. These explanations were also suitable for other estimated coefficients. Thus, the local bivariate analysis allowed us to capture the extent of local heterogeneity between habitat quality and different environmental situations: 27.37–43.14% in BTH, 23.7–65.37% in GBA, and 27.95–55.84% in YRD.

Comparing the performance of four environmental situation indicators, the spatial influences of SHII, NDVI, and RS on habitats in the neighboring local setting were more apt to be detected. However, understanding how the PM2.5 concentration locally influenced the UAs' habitat over space was not completed due to there being over 72% of non-significant areas. Looking at the statistically significant features, the negative linear relationship accounted for the largest share in the HQ–SHII, HQ–PM2.5, and HQ–RS results. In contrast, the positive linear relationship was the most remarkable in the HQ–NDVI results.

**Figure 12.** Local spatial relationship between habitat quality and multifaceted environmental states: (**a**–**c**) HQ–SHII for BTH, GBA, and YRD, respectively; (**d**–**f**) HQ–PM2.5 for BTH, GBA, and YRD, respectively; (**g**–**i**) HQ–RS for BTH, GBA, and YRD, respectively; and (**j**–**l**) HQ–NDVI for BTH, GBA, and YRD, respectively.

**Table 2.** Categorical summary of the local bivariate analysis.

### *3.4. GeoDetector-Based Interactive Effects Assessment*

Spatial stratified heterogeneity and interactive influence between habitat quality and multi-dimensional environmental status in the three UAs were captured using GeoDetector. The outputs show that the geographically stratified heterogeneities of habitat quality in the three UAs were, to a certain extent, affected by multifaceted environmental status. The habitat quality varied significantly over space for each graded environmental status indicator based on a series of *t*-tests at a significance level of 0.05. The related outputs of factor and ecological detector modules are shown in the Supplementary Materials. The stratified average values of habitat quality in each stratum of environmental situations for the three UAs are shown in Figure 13. The *q* statistic coefficients for four graded environmental status indicators in the three UAs under the factor detector module are shown in Figure 14. The higher *q* statistic coefficient indicated a greater explanatory power for the stratified heterogeneity of habitat quality. According to the coefficients of the *q* statistic, the most crucial influencing indicators on habitat quality in the three UAs were PM2.5 concentration (BTH, *q* = 0.17), SHII (GBA, *q* = 0.30), and NDVI (YRD, *q* = 0.21). Other indicators also showed a considerable impact on habitat quality. The outputs of the interaction detector module revealed the interactive influences from a combination of multiple environmental status indicators on habitat quality, as shown in Figure 15. With reference to the judging criteria of interaction (Table 1), we found that the *q* statistics for the bivariate interactions were greater than the *q* statistics for every single factor, indicating that the explanatory power of a single indicator could be mutually enhanced when interacting with another one. The cumulative effects of different environmental statuses could bring about more pressures and risks on habitat quality. Specifically, in BTH, the PM2.5 concentration played a dominant role in the single-factor influence analysis and strongly interacted with other indicators. The predominant interaction between PM2.5 concentration and NDVI showed the highest *q* statistic. A similar phenomenon was also found in the SHII indicator of GBA and the NDVI indicator of YRD.

**Figure 13.** The performances of the graded environmental situations under the risk detector: (**a**) SHII, (**b**) PM2.5, (**c**) RS, and (**d**) NDVI.

**Figure 14.** The *q*-statistic coefficients for the grading indicators in the three UAs.

### **4. Discussion**

### *4.1. The Impacts of Multifaceted Environmental States on Habitat Quality in the Three Urban Agglomerations*

Focusing on BTH, GBA, and YRD, this study spatially evaluated the relationships between habitat quality and multifaceted environmental states using various methods. The cartographic comparison allowed us to capture spatial differentiation and hotspots, the correlation analysis provided a global overview of the statistical relationships, the local entropy model inspected the local structural associations, and GeoDetector identified the spatial stratification heterogeneity and interactive influence. Through such incorporation of multi-perspective profiling, this study presents a comprehensive panorama of the interconnection of various environmental circumstances and their impacts on habitats for different UAs. As a series of biophysical and anthropogenic drivers can threaten the habitat quality and simultaneously trigger changes in different environmental mechanisms, there were similarities and discrepancies in the performances of the three UAs, reflected by the dominant impact and relationship strengths.

Based on the overlay of habitat quality and land/cover use maps in the three UAs, we found that more than four-fifths showed low–medium habitat suitability, spatially clumped in the human settlement and cultivation areas. Meanwhile, large-scale UF areas displayed prominent heat island phenomena, terrible air quality, over-intensive living space, and sparse green availability, albeit with varying degrees of three UAs. Habitat types with high greenness (e.g., forests) in the three UAs showed small gaps in habitat quality scores, ranging from 0.409 to 0.467. The scores of habitats with low greenness (e.g., grass, water) greatly varied across different UAs, depending on habitat types and extents (Figures 6 and 13d). Grounded on correlation analysis, SHII, PM2.5 concentrations, and RS showed an overall deteriorative impact, whereas NDVI had a beneficial influence on the habitat quality in the three UAs. This result was in line with previous environmental-related studies [10,14,15,44,49,74,75]. We also recognized that the most noticeable impacts on habitat were air pollution in BTH, thermal stress, and green capacity in GBA, and residential pressure in YRD. Thus, habitat quality tied in closely with natural assets' situations, indicating that habitat-related strategies depended on improving and maintaining the green and blue spaces.

Next, bivariate relationship profiling based on the local entropy model identified and delineated the local spatial influence of individual environmental status on habitat quality. The spatial mapping of local relationships showed the advantage of distinguishing diverse types of significant associations and highlighting their spatial patterns. Local spatial relationships were roughly in agreement with the global trend of the entire UA. Except for non-significant areas, the habitat quality scores of most regions in the three UAs were negatively related to SHII, PM2.5 concentrations, and RS, while positively associated with NDVI at the 95% confidence level. The local spatial influences of different environmental statuses on habitats were complicated with a rich mixture of relationship types, regardless of UAs. However, large shares of non-significant areas emerged in terms of HQ–PM2.5 (>72%) and HQ–RS (>51%). The PM2.5 concentrations used in this study were estimated based on in situ observations, AOD distribution, and auxiliary data (e.g., land cover/use and population) [57,58]. Our habitat quality scores were assessed using the land cover/use-based InVEST–HQ module [60]. It is difficult to simulate the habitat suitability and quality considering air-related threats due to model constraints, which interfere with the relationships captured between habitat and PM 2.5 concentrations. Regarding HQ– RS performance, the non-significant areas mainly belonged to natural habitats without residential function and service.

The upshots from GeoDetector analysis generally conformed with the foregoing paired relationship profiling and demonstrated the pile-up effect of any two different environmental issues. Relying on the GeoDetector investigation, this study emphasized spatial differences in the overlapping effects and interaction relationships of differently graded environmental strata on habitat quality. Using the interaction detector, we figured out that the impacts of the combination of different environmental issues on habitats in the three UAs were bivariate-enhanced rather than independent or weakened, indicating the aggravation and amplification of compound effect compared with the influence of a single environmental status. In terms of the BTH region, the cumulative effects of air pollution in conjunction with other environmental issues were more severe, coinciding with previous investigations [57,76,77]. BTH has struggled with the worst consequences of air pollution due to a chain of factors [16,78–80]. BTH is the national capital region of China, with over 110 million people [45]. The primary contributors were high energy consumption, vehicle emissions, biomass burning, and neighboring heavy industries [79,80]. As surrounded by mountains, BTH was a victim of its topography, where pollutants are trapped within the regional limits (mainly in the south and east of BTH).

In contrast, the overlapping influences of the heat island phenomenon and other issues in GBA have attracted more attention. We identified that the average annual daytime and nighttime LST values in GBA were the highest among the three UAs. On the other hand, a vast extent with low habitat quality scores in GBA was spatially congregated in the three contiguous megacity clusters (Dongguan–Guangzhou–Foshan, Zhuhai–Macao, and Shenzhen–Hong Kong). Meanwhile, these three contiguous megacity clusters exhibited forceful residential service capacities and low green gradients. However, GBA was superior to the other two UAs in the scale, extent, and gradient of air pollution. Most areas with over 50 ug/m3 annual PM2.5 concentrations in GBA were merely spatial distributed in Foshan. Accordingly, we observed the lowest impact of PM2.5 concentrations on GBA's habitats. In the YRD case, the interactions between NDVI and other indicators stood out. Additionally, there were minor gaps in impact magnitudes between indicators for YRD in comparison to those for BTH and GBA. We noticed that the gradient of NDVI and habitat quality in YRD were inferior to that of the other two UAs. These performances of YRD can be partly ascribed to development gradients and geomorphic features. Owing to the fast-growing urban built-up and settlements in recent decades, YRD is the most populated, urbanized, and prosperous region in China, with numerous megacities and more than 80 million urban dwellers [45]. Most areas of YRD render a plain terrain, whereas uplands and forests are spatially configured in the south region. Shanghai, Hangzhou Bay, and Suxichang Metropolitan Areas are the most remarkable regions with high urbanness. Such geographical configuration of YRD could respond to its layout and composition of habitat quality, revealing the lag and importance of YRD's greenness and habitat security.

### *4.2. Road Ahead and Implication for the Spatial Management of Multi-Dimensional Environmental Issues*

The uniqueness of this study was the multi-perspective profiling of numerous environmental issues and their interactions. Unlike other environmental studies, this study discerned the associations and interactive effects of different environmental situations instead of clarifying the mechanisms and driving forces of individual environmental problems. Initially, the circumstances of the designated onefold environment in the three UAs were separately evaluated using different indicators and methods. Subsequently, we constructed an integrated framework that combined multiple approaches to link the varying mixes and degrees of different environmental impacts on habitats. This framework propels the multi-objective sustainable management of natural-social capital and co-benefits-seeking policies.

The spatially explicit quantification, mapping, and profiling of various environmental impacts and their relationships serve as the chief references for development plans and policymaking. In recent years, numerous studies detected trade-offs and synergies among multiple environmental services, demonstrating how to move towards sustainability by optimizing the governance and service deliveries of different ecosystems [24,81]. Analogously, our understanding of multifaceted environmental situations and their influences on habitat quality can support work toward tailored sustainable pathways to specify (1) whether and where to adopt the countermeasures to safeguard nature and human habitats, (2) which sector and whose situation should be accounted for and take the corresponding precedence, and (3) how to reduce overlapped risks and advance the co-benefits, which we examined in the context of the three largest Chinese UAs.

Our evaluation results for these three UAs suggest the weightiness of embracing the spatial diversities and context heterogeneities of regional environmental situations for spatially targeted decision-making. Importantly, the landscapes delivering high habitat quality (score ≥0.7) in the three UAs occupied only under 15% of their respective areas. However, over half areas emerged with varying strengths of deteriorative environmental consequences across different aspects. This study focused on four environmental issues, namely, heat-related stress, air pollution, over-intensive living space, and limited urban afforestation, which have been the most common and urgent problems in UA regions and emphasized in national initiatives and local policies [23,25]. We identified that the compound impacts were amplified when one environmental issue overlapped with another. Since the interactive effects of different environmental issues were not independent, strategies to improve an onefold environment may contribute to ameliorating other environmental situations. It is imperative to take advantage of the interactive traits of environmental issues to safeguard and advance regional sustainability, preventing the formation and magnification of multifold risks from ecological, climate, and social aspects.

Furthermore, the motivation for UA formation is integrating regional resources, complementarity advantages, and beneficial cooperation [6]. The geographical heterogeneities and multi-dimensional gaps of UAs can offer opportunities for establishing rational resource allocation and cost-sharing mechanisms [6]. The core megacities and highly populated areas of the UA have been subjected to multiple environmental pressures. Due to the extremely constrained habitat availability and ecological carrying capacities, the supply of essential ecological services in these areas should be addressed by the enrichment and complementary of surrounding natural/semi-natural ecosystems rather than self-sufficiency [82]. Therefore, it is crucial to review the heterogeneities and relationships of the multifaceted environmental situations across the entire UA. By synthesizing different situations, the corresponding geographic areas should be identified to efficiently demarcate vulnerable environmental zones, implement the scopes of mitigation or restoration, and equitably assign environmental compensation [24,82].

Considering different hotspots and relationships between environmental issues, the spatial targeting of specific policies should be deliberated in potentially vulnerable areas to ensure effectiveness [38,83,84]. Spatially explicit mapping should visualize the performance of habitat quality and various environmental states, which can be more prone to identifying corresponding vulnerable locations. The most suitable strategies should be generated by integrating the different environmental profiles and minimizing overall negative feedback [23]. For example, since the most severe problem in BTH is air pollution, the interaction between air pollution and other issues can be aggravated. The mitigation of air pollution should be a considerable priority [19,79,80]. A series of strategies and interventions, including a clean energy plan, vast investments in urban greenery, the optimization of energy consumption structure and efficiency, and technological innovations [53], can take precedence to ameliorate the pollution condition of BTH. Simultaneously, thermal and biophysical environments in BTH can also be improved due to common management strategies. Similar arguments are suitable for GBA and YRD. A relief of the heat island phenomenon in GBA and the management intervention of urban afforestation in YRD should take precedence. Prioritization of decision-making in tandem with compound risk minimization sheds light on a promising amelioration pathway of both onefold and comprehensive UA environments, which can make the achievements of spatially comprehensive management and multi-dimensional sustainability possible and efficient from various aspects [19,23,82].

Although we are confident about the contributions of this study to sustainabilityrelated policymaking processes, there is still enormous potential for our framework. Firstly, we are considering the improvement of the habitat quality quantification in further studies by incorporating more threat source data and models, such as traffic data and heavy

industrial zone. Our habitat quality was mainly derived from land cover/use evaluation and, to some extent, affected the investigation of relationships between habitat quality and certain environmental situations (e.g., air quality). On the other hand, there can be a dynamic improvement when updating the indicators of environmental issues from a spatiotemporal perspective. In this study, restricted by data availability and consistency, we only evaluated the situations and impacts of four environmental issues on the spatial grid scale in 2015. The relationships between environmental issues and their influences may vary across different spatial scales and evolve over time. We desire to add more temporal dynamics and incorporate policy information to verify the proposed comprehensive research framework further. In particular, significant spatiotemporal disparities and trends may be reflected when comparing the situations before and after the COVID-19 pandemic. Moreover, we can endeavor to explore more environmental issues to enrich the multifaceted evaluation. Environmental profiles, such as energy emission, water purification, and food security, can be synthesized into our future research framework to produce more insights for sustainability support.

### **5. Conclusions**

In tailoring sustainable development schemes for the three UAs of BTH, GBA, and YRD, we profiled the relationships and interactions between habitat quality and multienvironmental status indicators from various perspectives by integrating cartographic comparison, correlation analysis, local entropy map, and GeoDetector. The cartographic comparison highlighted geographical hotspots, the correlation analysis weighed global heterogeneities, the local entropy map differentiated the local spatial associations, and GeoDetector recognized interactive relationships.

The contributions of this study stem from four aspects: (1) spatial valuation of habitat quality for the three UAs; (2) spatial overview of the multifaceted environmental situations for the three UAs; (3) multi-perspective profiling of the relationships and interactive influences of different environmental situations on habitat quality for the three UAs; and (4) implications for supporting habitat sustainability and co-benefit-sought decisions. The critical discoveries are summarized as follows:

In 2015, most of the regional landscapes in the three UAs delivered low/medium habitat availability levels overlapped with threats from the thermal, air, living, and biological environmental spheres. The worst compound environmental risks were geographically aggregated in highly urbanized areas, showing an execrable habitat quality. As expected, the habitat security of the UAs was diminished by the heat island phenomenon, air pollution, and residential support. In contrast, enhancing afforestation in the UAs can be a promising pathway to maintain habitat health. The delineation of the local structural relationships between habitat quality and different environmental statuses highlighted the spatial complexity of the influences and interplay among various environmental systems, indicating potential environmental injustice, spatially heterogeneous changes, and context dependency. Strategies targeting specific environmental issues may induce spatially varying consequences on habitats and produce various trade-offs in other environmental systems. The interaction between various environmental issues manifested bivariate-enhanced effects, implying that the overlay of different environmental issues can magnify the influence of each environmental risk. The seriousness of the worst environmental issue would deteriorate when overlapping with others. Thus, the prioritization of decision-making processes may be required in different geographical contexts, such as air pollution removal in BTH, heat island mitigation in GBA, and green space improvement in YRD.

In conclusion, the findings of this study support the multi-dimensional sustainability and co-benefits sought by taking advantage of the relationships and interactions between habitats and diverse environmental issues in the UAs.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/rs15040921/s1, Table S1. Threat parameter settings for estimating habitat quality; Table S2. The sensitivity parameter settings of habitat types to threat factors; Table S3. Entropy results summary for local bivariate analysis; Figure S1. Spatial stratification of the land surface temperatures (LSTs) in the three UAs; Figure S2. Spatial stratification of the surface heat islands in the three UAs; Figure S3. Results for ecological detector module; Figure S4. Results for risk detector module.

**Author Contributions:** Conceptualization, F.L. and Y.M. (Yuji Murayama); Data curation, F.L.; Formal analysis, F.L.; Funding acquisition, Y.M. (Yuji Murayama); Investigation, F.L.; Methodology, F.L.; Project administration, Y.M. (Yuji Murayama); Resources, F.L., Y.M. (Yuji Murayama) and Y.M. (Yoshifumi Masago); Software, F.L.; Supervision, Y.M. (Yuji Murayama); Validation, F.L.; Visualization, F.L.; Writing—original draft, F.L.; Writing—review and editing, Y.M. (Yuji Murayama) and Y.M. (Yoshifumi Masago). All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partly supported by the Japan Society for the Promotion of Science (JSPS) grants of 21K01027 and 18H00763.

**Data Availability Statement:** The data that support the findings of this study are available upon reasonable request from the corresponding author.

**Acknowledgments:** The authors highly appreciate the editors and anonymous reviewers for their valuable and constructive comments and suggestions to improve the manuscript.

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

### **References**

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

### *Article* **Driving Factors of Land Surface Temperature in Urban Agglomerations: A Case Study in the Pearl River Delta, China**

**Wenxiu Liu 1,2, Qingyan Meng 1,2,3,\*, Mona Allam 4, Linlin Zhang 1, Die Hu 1,2 and Massimo Menenti 5**

	- Stevinweg, 12628 CN Delft, The Netherlands<EMAIL_ADDRESS>
**Abstract:** Land surface temperature (LST) in urban agglomerations plays an important role for policymakers in urban planning. The Pearl River Delta (PRD) is one of the regions with the highest urban densities in the world. This study aims to explore the spatial patterns and the dominant drivers of LST in the PRD. MODIS LST (MYD11A2) data from 2005 and 2015 were used in this study. First, spatial analysis methods were applied in order to determine the spatial patterns of LST and to identity the hotspot areas (HSAs). Second, the hotspot ratio index (HRI), as a metric of thermal heterogeneity, was developed in order to identify the features of thermal environment across the nine cities in the PRD. Finally, the geo-detector (GD) metric was employed to explore the dominant drivers of LST, which included elevation, land use/land cover (LUCC), the normalized difference vegetation index (NDVI), impervious surface distribution density (ISDD), gross domestic product (GDP), population density (POP), and nighttime light index (NLI). The GD metric has the advantages of detecting the dominant drivers without assuming linear relationships and measuring the combined effects of the drivers. The results of Moran's Index showed that the daytime and nighttime LST were close to the cluster pattern. Therefore, this process led to the identification of HSAs. The HSAs were concentrated in the central PRD and were distributed around the Pearl River estuary. The results of the HRI indicated that the spatial distribution of the HSAs was highly heterogeneous among the cities for both daytime and nighttime. The highest HRI values were recorded in the cities of Dongguan and Shenzhen during the daytime. The HRI values in the cities of Zhaoqing, Jiangmen, and Huizhou were relatively lower in both daytime and nighttime. The dominant drivers of LST varied from city to city. The influence of land cover and socio-economic factors on daytime LST was higher in the highly urbanized cities than in the cities with low urbanization rates. For the cities of Zhaoqing, Huizhou, and Jiangmen, elevation was the dominant driver of daytime LST during the study period, and for the other cities in the PRD, the main driver changed from land cover in 2005 to NLI in 2015. This study is expected to provide useful guidance for planning of the thermal environment in urban agglomerations.

**Keywords:** land surface temperature; spatial analysis; urban agglomeration; driving factors; geodetector metric

### **1. Introduction**

The urban heat island (UHI) is a global phenomenon caused by urbanization [1]. UHI affects air quality [2], threatens the health of urban residents [3,4], influences building energy consumption, and leads to the risk of overheating in outdoor thermal environments [5]. Generally, UHI can be assessed by using the air temperature or land surface temperature

**Citation:** Liu, W.; Meng, Q.; Allam, M.; Zhang, L.; Hu, D.; Menenti, M. Driving Factors of Land Surface Temperature in Urban Agglomerations: A Case Study in the Pearl River Delta, China. *Remote Sens.* **2021**, *13*, 2858. https://doi.org/ 10.3390/rs13152858

Academic Editors: Yuji Murayama and Ruci Wang

Received: 16 June 2021 Accepted: 18 July 2021 Published: 21 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

(LST). Urban agglomerations, which represent groups of cities that have a compact spatial organization and close economic connections, have become the most prominent feature of global urbanization in recent decades [6,7]. Especially in China, urban agglomerations have become the major form of national urbanization [8]. Therefore, it is important to understand the spatial distribution patterns of LST and the main driving factors affecting LST in order to formulate informed urban policies for urban agglomerations in the future [9].

Most previous studies focused on the spatial distribution characteristics and driving factors of LST [10–18]. The spatial distribution of high LSTs is significant for understanding the landscape patterns of LST [19,20]. According to previous studies, the potential factors of LST can be divided into three types: (1) land cover factors, which include land use/land cover (LUCC) [21], the normalized difference vegetation index (NDVI) [22–24], and the normalized difference built-up index (NDBI) [25,26]. This type also includes the impervious surface area (ISA) and the impervious surface distribution density (ISDD), which indicate the urbanization intensity of an area from the perspective of impervious surfaces [27,28]. (2) Socio-economic factors, which include the gross domestic product (GDP) [29], population density (POP) [30,31], and industrial production activities [32]. In addition, the nighttime light index (NLI) can indicate the human density [33], energy consumption [34], built-up area [35], and socioeconomic dynamics [36]. However, few studies have explored the relationship between NLI and LST. (3) Physiographical factors, which include elevation and slope, have also been studied [37,38].

The correlations between LST and various factors should be clarified in order to determine the main drivers of LST. Several methods, including the use of Pearson correlation coefficient [39], linear regression analysis [8,30], and stepwise multiple-linear regression [40], have been widely used in previous studies. LST is usually affected by multiple factors, and their interactions are quite complex. However, the above-mentioned methods cannot adequately measure the non-linear correlations between LST and the various drivers [9], or they cannot deal with categorical data [41]. In addition, few studies have taken the combined effect of different factors into consideration when determining the dominant factors of LST. Therefore, this study applied the geo-detector (GD) metric [42] to explore the correlations between LST and three types of drivers in order to avoid the previous limitations. The GD metric is applicable to the non-linear relations between the drivers and LST, and can also measure the combined impact of different drivers on LST. Due to its good applicability, the GD metric has been widely used in studies of drivers, including air pollution [43–45], public health [41,46,47], land use [48,49], regional economies [50,51], and urban environment [9,52].

This study was conducted in the Pearl River Delta (PRD) urban agglomeration, China. The PRD is one of the most important economic centers in China, and is also one of the most highly urbanized and populous regions in the world [53]. The cities included in the PRD were merged into an urban agglomeration in 2005 [54]. The PRD contributed about 9.12% of the national gross domestic product in 2015 [55], and is known as an "economic miracle" and the "world factory". Its rapid urbanization and industrialization led to significant land cover changes and brought a series of environmental problems, including urban flood risk, air pollution, and the UHI effect in the PRD [56–58]. The PRD has a subtropical climate and is located along the coast of the South China Sea. The dual effects of rapid urbanization and climate make the PRD highly susceptible to extreme heatwave events [20]. However, the current studies on the thermal environment in the PRD are still not adequate. Most previous studies focused solely on individual factors of LST, such as LUCC, ISA, or greenspace in the PRD [19,20,37,59–62]. In addition, more attention has been paid to the core cities, such as Guangzhou, Dongguan, Shenzhen, and Zhongshan cities [19,63–66], while the cities in the east and west of the PRD have received less attention.

To address the above problems, this study aimed to explore the spatial distribution characteristics of LST and to detect its main driving factors by using the GD metric in the PRD. First, we studied the spatial distribution of LST and used Moran's Index to indicate the presence of hotspot areas (HSAs). Then, the HSAs were determined by using hotspot analysis in a GIS environment. Second, the Hotspot Ratio Index (HRI) was developed to evaluate and rank the thermal gradient of each city in the urban agglomeration. Finally, the GD metric was used to measure the influence of each factor and the combined influence of the different factors on LST. Seven factors were selected based on the literature review and available data. The effects of the seven factors on the LST in the nine cities were compared to reveal the complex mechanism of LST in the PRD. Our study contributes to a better understanding of the spatial distribution characteristics of LST and its main driving factors, which is helpful in providing insights into the optimization of the urban thermal environment.

### **2. Study Area and Data Sources**

### *2.1. Study Area*

The PRD urban agglomeration is one of the regions with the highest urban densities in the world and is one of the most developed regions in China. It is located in the southcentral coastal region of Guangdong Province, and it includes nine cities: Guangzhou, Foshan, Dongguan, Shenzhen, Zhongshan, Zhuhai, Huizhou, Zhaoqing, and Jiangmen (Figure 1). The PRD is surrounded by hills and mountains to the north, west, and east. The total area of the PRD is 55,000 km2 [61]. The PRD experienced rapid population growth and economic growth from 2005 to 2015, with its population increasing by 13.27 million, and its GDP increased by about 3.8 trillion yuan (http://stats.gd.gov.cn, accessed on 15 September 2020). The PRD became the largest mega-region in the world according to the World Bank Group (2015). The average annual temperature ranges from 21 to 23 ◦C, and the average annual precipitation is over 1500 mm [60,67]. The PRD was selected for this study because its thermal environment became a serious problem due to the acceleration of its urbanization. Table 1 shows some key attributes of the nine cities in the PRD. The air temperature data were obtained from the National Meteorological Data Center (http://data.cma.cn, accessed on 10 April 2021). The permanent population and electricity consumption data were derived from the Statistical Yearbook of Guangdong Province.

**Figure 1.** Maps of the Pearl River Delta (PRD) urban agglomeration: (**a**) map of Guangdong province, China; (**b**) PRD urban agglomerations, including nine cities; (**c**) elevation map of the PRD.

Note: "+" indicates the positive changes from 2005 to 2015, and "−" indicates the negative changes from 2005 to 2015; the rate of change is calculated with the formula: (T2015- T2005)/T2005\*100.

### *2.2. Data Sources*

In this study, the LST data for 2005 and 2015 were obtained from the 8-days (version 6) MODIS Aqua LST composite products (MYD11A2) (https://search.earthdata.nasa. gov, accessed on 6 June 2019). The MYD11A2 LST product was validated with in situ measurements to yield a bias of <0.5 K [68]. We used the annual mean LST observed in the daytime/nighttime at 13:30/01:30 local solar time [37]. The potential driving factors were selected according to previous research and available data, including the elevation, LUCC, NDVI, ISDD, GDP, POP, and NLI (Table 2). The elevation, LUCC, NDVI, POP, and GDP data were obtained from the Resource and Environment Science and Data Center (www.resdc.cn, accessed on 15 June 2019). The LUCC data included six types: cultivated land, woodland, grassland, water area, construction, and unused land. The NDVI product was computed from the continuous time series of the SPOT/VEGETATION NDVI remote sensing dataset. The GDP and POP were grid data with a 1 km resolution. In addition, the impervious surface data were obtained from published global urban datasets (https://doi.org/10.6084/m9.figshare.11513178.v1, accessed on 18 June 2020). The calculation of the ISDD is detailed in [27]. The NLI was obtained from nighttime light imageries, which were downloaded from the Resource and Environment Science and Data Center (www.resdc.cn, accessed on 15 June 2019) and the National Geophysical Data Center (http://ngdc.noaa.gov, accessed on 15 June 2019). All data were resampled to a 1 km resolution to keep the consistence of the spatial resolution of the data analyzed. The sources of the driver data are summarized in Table 2.

### **Table 2.** Influencing drivers and data sources.

### **3. Methods**

### *3.1. Spatial Variability of LST*

The MYD11A2 LST products were used to get the annual average daytime LST and annual average nighttime LST for 2005 and 2015. The MODIS Reprojection tool was used to preprocess the MYD11A2 LST data, which included format conversion, re-projection, and clipping. The mask tool was used to remove cloud pixels; only pixels with a high quality (average LST error < = 1 K) were used to calculate the annual average LST for 2005 and 2015.

### 3.1.1. Global Moran's Index

Spatial distribution analysis is often used to explore geographical phenomena. The Global Moran's Index was used to analyze the spatial distribution of the LST. The value of the Global Moran's Index is between −1.0 and 1.0, where 1 indicates perfect positive spatial autocorrelation, −1 indicates perfect negative spatial autocorrelation, and 0 indicates a perfect spatial randomness [69]. The Global Moran's Index is calculated as follows:

$$I = \frac{n\sum\_{i=i}^{n}\sum\_{j=1}^{n}\mathcal{W}\_{ij}(\mathbf{x}\_{i}-\overline{\mathbf{x}})\left(\mathbf{x}\_{j}-\overline{\mathbf{x}}\right)}{\sum\_{i=1}^{n}\sum\_{j=1}^{n}\mathcal{W}\_{ij}\sum\_{i=1}^{n}\left(\mathbf{x}\_{i}-\overline{\mathbf{x}}\right)^{2}}\tag{1}$$

where *n* is the number of pixels; *xi* and *xj* are the LSTs at pixels *i* and *j*; *x* is the mean LST; *wij* is the spatial weight determined with the spatial correlation of the LSTs at pixels *i* and *j*.

### 3.1.2. Hotspot Analysis

Many previous methods for delineating high-temperature regions have applied numerical classification [70], Gaussian surface fitting [71–73], radial sampling [74,75], and hotspot analysis [76,77]. Hotspot analysis can be used to identify the spatial clusters of high LSTs (hotspots) and low LSTs (cold spots) by using only LST data, thus avoiding the subjective influence that may be caused by having too many parameter settings. Therefore, hotspot analysis was applied to identify the areas with high (hotspot) or low (cold spot) LSTs in this study. Based on the clustering patterns of the spatial distribution of the LST, the Getis–Ord Gi (*Gi\**) index was applied to measure the degree of clustering for each pixel [77]. For each pixel, the hotspot analysis returns a z-score. A higher positive z- score indicates a higher degree of clustering of high LSTs (hot spots) [78]. In this study, the areas with z-scores ≥ 2.58 (corresponding to the 99% confidence level) were defined as HSAs [76]. The *Gi\** is defined as:

$$G\_i^\* = \frac{\frac{\sum\_{j=1}^n \omega\_{ij} x\_j - \overline{X} \sum\_{j=1}^n \omega\_{ij}}{S \sqrt{\frac{\left[n \sum\_{j=1}^n \omega\_{ij}^2 - \left(\sum\_{j=1}^n \omega\_{ij}\right)^2\right]}{n-1}}} \tag{2}$$

$$\overline{X} = \frac{\sum\_{j=1}^{n} x\_j}{n} \tag{3}$$

$$S = \sqrt{\frac{\sum\_{j=1}^{n} \mathbf{x}\_{j}^{2}}{n} - \left(\overline{X}\right)^{2}}\tag{4}$$

where, *xj* is the LST value of pixel *j*, *ωij* represents the spatial weight between pixels *i* and *j*, and *n* is the total number of pixels; the spatial weight is determined according to the Queen's adjacency connectivity matrix [76].

### *3.2. Hotspot Ratio Index*

The urban-heat-island ratio index was applied in previous studies to evaluate and rank thermal environments [79,80]. It is effective in measuring the characteristics of a thermal environmental, but its calculation requires urban or rural boundary data [80,81]. In this study, the hotspot ratio index (HRI) was used to evaluate and rank the thermal environments of the nine cities. The HRI gives scores and ranks the areas according to the value of the LST. The HRI is calculated as the weighted sum of the percentage of HSAs with different classes in the selected region. The HRI reflects the thermal gradient in the selected region. A higher HRI indicates that a city has a higher thermal gradient, which implies that the city is more likely to experience greater thermal stress. The HRI is calculated as follows:

$$HRI = \sum\_{i=1}^{n} i \ast P\_i \tag{5}$$

where *n* is the total number of classes. To avoid skewed data, the quantile classification method was applied to classify the HSAs into 5 classes according to the LSTs. *Pi* is the ratio of the area of class *i* to the total study area. *i* is the class index, i.e., *i* = 1, 2, 3, 4, and 5. A higher class index refers to a higher LST. To calculate the HRI, the areas of 1st, 2nd, 3rd, 4th, and 5th classes of HSAs were divided by the total study area, then multiplied by the class index.

### *3.3. Geo-Detector Metric*

The geo-detector (GD) metric is an effective metric for revealing the drivers of geographical phenomena [42]. The main idea assumes that if an independent variable affects the dependent variable, then their spatial distributions should be similar. The GD metric includes four detectors: factor, interaction, risk, and ecological. We used the factor detector to measure the influence of each driver on the LST. We also used the interaction detector to measure the combined influence of the interactions between drivers on the LST. The GD metric is freely available from www.geodetector.cn, accessed on 18 July 2021.

### 3.3.1. Factor Detector

The factor detector signifies the effects of factors on the LST. The greater the probability distribution (q) value is, the greater the influence of the factors on the LST will be. The value of q is between 0 and 1; in extreme cases, the q value equals 1, indicating that factor X completely affects the spatial distribution of Y, and a q value of 0 indicates that factor X has nothing to do with Y. The value of q is calculated as follows [42]:

$$q = 1 - \frac{1}{N\sigma^2} \sum\_{h=1}^{L} N\_h \sigma\_h^2 \tag{6}$$

where the LST (Y) and the drivers (X) are composed of *L* classes (*h* = 1, 2, ... *L*); *N* and *Nh* represent the numbers of cells in the entire area and the *h* class, and *σ*2 and *σ*2 *h* represent the variance in the LST.

### 3.3.2. Interaction Detector

The interaction detector was used to measure the influence of the interactions between different factors on the LST. It identifies whether the interactions between drivers (X1 and X2) have an effective influence on the LST or not. First, we separately calculate the influence of each driver (q(X1) and q(X2)) on the LST; then, the influence of the interaction between X1 and X2 (q(X1∩X2)) on the LST is calculated based on an overlay of the two drivers (X1 and X2). According to the values of q(X1), q(X2), and q(X1∩X2), the effect of the interaction between the two drivers on the LST can be determined (Table 3) [9].

**Table 3.** The effects of interactions between two drivers.

### 3.3.3. Data Preparation

The GD metric captures the spatial heterogeneity of the attributes, and it requires the discretization of input data. Many methods have been applied for discretization [82]. The natural breaks method depends on the principle of the "maximum and minimum distance between classes" that can allow the original characteristics of the data to be kept [83]. This method has provided good results [43,52]. Accordingly, we applied natural breaks classification in order to classify the drivers. Sample points were generated in ArcGIS according to the 1 km grids over the entire study area. Therefore, the numbers of sample points in each city depended on the area of the city. The values of each driver and mean LST

were extracted for the sample points. The natural breaks method was applied to classify all drivers (except LUCC, as it was already classified into cultivated land, woodland, grassland, water areas, construction, and unused land). For the GD metric, LST was the dependent variable, while the classified drivers were independent variables (Table 2).

### **4. Results**

### *4.1. Spatial Distribution of LST*

As seen in Figure 2, in 2005, higher daytime LSTs were mostly concentrated in the central area of the PRD, including the cities of Dongguan, Shenzhen, Guangzhou, Foshan, and Huizhou. By contrast, in 2015, higher daytime LSTs were concentrated in Dongguan, Shenzhen, Foshan, Zhongshan, Jiangmen, and Guangzhou. The daytime LSTs decreased in the cities located in the northern part of the PRD, including the cities of Zhaoqing, Guangzhou, and Huizhou. For the nighttime LST, the higher LSTs were also concentrated in the central area of the PRD in 2005. In 2015, the highest nighttime LSTs increased in the cities of Dongguan, Shenzhen, Foshan, Zhongshan, Jiangmen, and Guangzhou.

**Figure 2.** Daytime and nighttime land surface temperature (LST) maps of the PRD:(**a**) daytime LST in 2005; (**b**) daytime LST in 2015; (**c**) daytime LST difference between 2005 and 2015; (**d**) nighttime LST in 2005; (**e**) nighttime LST in 2015; (**f**) nighttime LST difference between 2005 and 2015.

4.1.1. General Spatial Pattern of LST

The spatial pattern of the LST was estimated by using the Global Moran's Index. The Global Moran's Indexes for the daytime LST were 0.903 and 0.915 at the 99% significance level in 2005 and 2015, respectively. For nighttime LST, the Global Moran's Indexes were 0.940 and 0.941 at the 99% significance level in 2005 and 2015, respectively. The results for the Moran's Index indicate that the LST distribution was a clustered pattern—areas had aggregations of high and low temperatures—and there was less than a 1% likelihood that the results were random, which indicates the presence of HSAs in the study area.

### 4.1.2. Spatial Distribution of Hotspot Areas (HSAs)

Hotspot analysis was applied to identify the HSAs in the PRD. For each LST pixel, the hotspot analysis returned a z-score. We defined areas with z-scores ≥ 2.58 as HSAs according to [78] (Table 4). For the daytime LST, the HSAs pattern was consistent with the pattern of the high-temperature regions (Figures 2 and 3).

**Figure 3.** Hot-/cold-spot distributions in the PRD: (**a**) hot-/cold-spots for daytime LST in 2005; (**b**) hot-/cold-spots for daytime LST in 2015; (**c**) hot-/cold-spots for nighttime LST in 2005; (**d**) hot-/cold-spots for nighttime LST in 2015.

Three patterns of HSAs in the PRD were observed—stable, reduced, and expanded by comparing the HSAs maps for 2015 with those for 2005 (Figure 4). For the daytime LST, the stable HSAs distribution indicates that the HSAs tended to be clustered in the center of the PRD. The HSAs expanded towards the central and southwestern parts of the PRD. A considerable growth of HSAs was observed in Jiangmen compared with Huizhou and Zhaoqing. In contrast, the HSAs were significantly reduced in the northwest and northeast of the PRD.

**Table 4.** Average LSTs and the percentage of hotspot areas (HSAs) in each city during the study periods.

Note: HSAs were delineated at the 99% significance level.

#### *Remote Sens.* **2021** , *13*, 2858

**Figure 4.** Maps of changes in the spatial distributions of hotspot areas (HSAs) at the 99% significance level in the PRD: (**a**) the changes in daytime HSAs in the PRD; (**b**) the changes in nighttime HSAs in the PRD (the "blank" areas did not include HSAs).

For the nighttime LST, the stable HSAs were bigger and more continuously distributed than those in daytime. Most of the stable HSAs were located near the coast and the Pearl River estuary, including in Zhuhai, Zhongshan, Foshan, Guangzhou, Shenzhen, and Dongguan.

### *4.2. Hotspot Ratio Index in the Nine Cities*

The HSAs were classified according to [79] (Table 5 and Figure 5). The values of the HRI varied across cities in the daytime (Figure 6). Low HRI values were observed in Zhaoqing, Jiangmen, Zhuhai, and Huizhou (Figure 6). The HRI values of Dongguan, Shenzhen, Foshan, and Zhongshan were higher than those of the other cities.

**Table 5.** Criteria for HSAs classification.

The HRI also varied across cities in the nighttime. The highest HRI values were recorded in Zhongshan and Zhuhai. The HRI values of Zhaoqing and Huizhou were minimal. The higher HRI values were recorded in Dongguan, Shenzhen, Zhongshan, Zhuhai, and Foshan cities. These results indicate that the spatial distribution of HSAs was highly heterogeneous among the cities in both daytime and nighttime.

### *4.3. Dominant Drivers of LST in the PRD*

The factor detector quantifies the influence of each driver on the LST by calculating the q values (Table 6). The q value shows that driver X explains (100\*q) % of the LST. The q

values at the 95% significance level are shown in Table 6. The results show the differences in the influence of each driver on the LST in 2005 and 2015.

**Figure 5.** HSAs classification in each city in 2005 and 2015: (**a**) the percentage of different levels of HSAs in the total area of each city in the daytime; (**b**) the percentage of different levels of HSAs in the total area of each city in the nighttime.

**Figure 6.** Hotspot Ratio Index (HRI) in the nine cities.

**Table 6.** The q values of the LST drivers in 2005 and 2015.

For the daytime LST, in 2005, elevation had the greatest influence on LST (50%), followed by the NDVI (40%). In 2015, the NLI was the dominant driver of LST, which indicated the increase in the influence of socio-economic drivers on daytime LST due to the continuous development of the urban agglomeration [84]. For the nighttime LST, in 2005, the NLI (45%) and POP (45%) had the greatest influence, followed by the NDVI (44%) (Table 6). In 2015, the NLI (50%) also had the greatest influence, followed by the NDVI (45%) (Table 6). The results indicate that the NLI and NDVI were important factors for the nighttime LST.

The explanation rate of the NLI for the daytime LST increased from 36% in 2005 to 58% in 2015, and the explanation of the NLI for the nighttime LST increased from 45% in 2005 to 50% in 2015, which indicated the important influence of the NLI factor on LST in the PRD.

### *4.4. Dominant Drivers of LST in the Nine Cities*

### 4.4.1. Factor Detector Analysis

The results for the factor detector for each city are presented in Figures 7 and 8, where a color-coding scheme was applied to facilitate the interpretation. Lower q values are shown in a lighter red color, while darker red colors indicate higher q values. There are two ways to explain these results. When read from left to right, the differences in the influences of all drivers in one city are clear. When reading from top to bottom, the differences across cities can be explained for each driver. The results show that the influence of the selected drivers on LST varied from city to city (Figures 7 and 8).

For Huizhou, Zhaoqing, and Jiangmen, the elevation had the greatest influence on daytime LST, as there are many mountains and forests in these three cities. Similarly, in Guangzhou, the elevation also exhibited the highest influence rate (61%) in 2005. Socioeconomic factors also showed great influences on nighttime LST in the above four cities (Figure 8).

For Dongguan and Shenzhen, the NDVI had the greatest influence on daytime LST in 2005, while in 2015 the NLI was the dominant driver. The NLI was the most influential factor for nighttime LST in Dongguan and Shenzhen in 2015, which indicated the importance of socio-economic development in these areas and its effect on LST. The NDVI also had the high explanation rate for nighttime LST in these two cities during the study period.

For Foshan and Zhongshan, the NLI was the dominant driver of daytime LST during the study period (Figure 7). The NDVI showed the greatest influence on nighttime LST in these two cities (Figure 8). The influences of the drivers on the nighttime LST in Guangzhou and Zhuhai are shown in Figure 8. None of the seven drivers had a significant influence on nighttime LST in Zhuhai, but all seven drivers had a strong influence on nighttime LST in Guangzhou.

The socio-economic drivers (GDP, POP, and NLI) were highly correlated with daytime LST in the central cities in the PRD, including Dongguan, Shenzhen, Guangzhou, Foshan, and Zhongshan (Figure 7). The socio-economic drivers had a low level of correlation with the daytime LST in Huizhou, Zhaoqing, and Jiangmen, where elevation was the dominant driver of daytime LST (Figure 7). Similarly, the land cover drivers, i.e., the LUCC, ISDD, and NDVI, were highly correlated with daytime LST in Dongguan, Shenzhen, Guangzhou, Foshan, and Zhongshan, while the correlation was weaker in Huizhou, Zhaoqing, and Jiangmen. Due to topographical issues, urbanization is concentrated in the central region

of the PRD [37]. For the nighttime LST, the NDVI showed a strong influence in most of the cities during the study period (Figure 8).

**Figure 7.** The q values of drivers for each city in the daytime (the blank values indicate that the results of corresponding factors are not significant at a *p*-value of 0.05).

**Figure 8.** The q values of drivers for each city in the nighttime (the blank values indicate that the results of corresponding factors are not significant at a *p*-value of 0.05).

### 4.4.2. Interaction Detector Analysis

In this study, the interaction detector was calculated for certain drivers: elevation, LUCC, NDVI, ISDD, GDP, POP, and NLI. Strong interactions were observed among all seven drivers (Figures 9–12). The interactions between elevation and the other drivers were significantly strong in Huizhou, Jiangmen, Zhaoqing, and Guangzhou (Figures 9 and 10). The greatest effects of the drivers' interactions on daytime LST were Elevation ∩ NDVI (68%), Elevation ∩ LUCC (46%), Elevation ∩ GDP (68%), and Elevation ∩ NDVI (71%) in Huizhou, Jiangmen, Zhaoqing, and Guangzhou, respectively, in 2005. In 2015, the greatest effects of drivers' interactions on daytime LST were Elevation ∩ NLI in Huizhou (69%), Jiangmen (57%), and Guangzhou (78%). The explanation rate of Elevation ∩ LUCC, Elevation ∩ ISDD, and Elevation ∩ POP for daytime LST in Zhaoqing was 63%, which illustrates that elevation had the greatest influence on daytime LST in this city.

**Figure 9.** The interactions among drivers with respect to daytime LST in each city (2005).

The interactions between elevation and the socio-economic factors showed a high explanation rate for nighttime LST in Zhaoqing, Jiangmen, Huizhou, and Guangzhou (Figures 11 and 12). The explanation rates of the interactions between elevation and the three socio-economic factors exceeded 30% (in Zhaoqing), 40% (in Jiangmen), 42% (in Huizhou), and 72% (in Guangzhou) in 2005. In 2015, the highest combined explanations of elevation and the socio-economic factors reached 39% (in Zhaoqing), 43% (in Jiangmen), 57% (in Huizhou), and 72% (in Guangzhou).

In Dongguan and Shenzhen, in 2005, the interactions between the NDVI and the three socio-economic factors exceeded 63% and 74%, respectively (Figure 9). In 2015, the greatest influence of the interactions on daytime LST in Dongguan and Shenzhen was from NLI ∩ LUCC (72%) and NLI ∩ NDVI (79%), respectively (Figure 10). The results of the interaction analysis for daytime LST in Dongguan and Shenzhen were consistent with the results of the factor detector analysis for daytime LST in these two cities. For the nighttime LST in Dongguan and Shenzhen, the predominant interactions between drivers were NDVI ∩ Elevation, NDVI ∩ NLI, NDVI ∩ POP, and NDVI ∩ GDP (Figures 11 and 12). This result indicates the important influence of the NDVI and socio-economic factors on the nighttime LST in these two cities.

**Figure 10.** The interactions among drivers with respect to daytime LST in each city (2015).

**Figure 11.** The interactions among drivers with respect to nighttime LST in each city (2005).

**Figure 12.** The interactions among drivers with respect to nighttime LST in each city (2015).

Similarly, the interactions between land cover and the socio-economic drivers were most predominant in Foshan, Zhongshan, and Zhuhai. NLI ∩ LUCC had the greatest influence on daytime LST in Foshan with 64% in 2005 and 63% in 2015. NLI ∩ NDVI had the greatest influence on daytime LST in Zhuhai with 43% in 2005 and 56% in 2015. The strongest interaction in Zhongshan was between the drivers NLI ∩ ISDD (65%) in 2005, and the interactions between NLI and the three land cover drivers explained 70% in 2015 (Figures 9 and 10). For the nighttime LST in Foshan, Zhongshan, and Zhuhai, the interactions between the NDVI and socio-economic factors were predominant during the study period (Figures 11 and 12).

### **5. Discussion**

### *5.1. Spatial Patterns and Drivers of LST*

This paper aimed to study the spatial distribution characteristics and main drivers of LST in the PRD urban agglomerations. The HSAs were mainly located on both sides of the Pearl River estuary. Our findings agree with those of [37]. The HSAs were mainly distributed in Shenzhen, Dongguan, Guangzhou, Foshan, Zhongshan, and Zhuhai, while there were fewer HSAs in Huizhou, Jiangmen, and Zhaoqing (Figure 4). The urbanization rates of the nine cities are as follows: Shenzhen > Guangzhou > Zhuhai > Dongguan > Foshan > Zhongshan > Huizhou > Jiangmen > Zhaoqing [85]. The HSAs results indicate that high temperatures are more likely to be observed in areas with high urbanization rates. There was a remarkable decrease in the daytime HSAs in the northern PRD (Figure 4). There are many mountains and forests in the northern PRD, where the environment has been improved significantly since 2006. A series of ecological projects have continuously promoted the increase in vegetation cover in these areas (http://gz.gov.cn, accessed on 15 September 2020), which decreased the LST in the northern PRD. This finding is consistent with the findings of [20]. There was a difference between the spatial distributions of daytime and nighttime HSAs. This could be explained by the impacts of impervious

surfaces and human activities on the LST in areas with higher urbanization rates during the day, while the impact of water became significant at night due to its high thermal capacity and properties [10]. In other words, the differences between daytime and nighttime HSAs may be related to the different driving forces of daytime and nighttime LST [86].

The HRI was developed to quantitatively evaluate and rank the thermal environment of each city. The results for the HRI indicated that Dongguan and Shenzhen, which were characterized by high urbanization rates, were the most thermally cities in day and night (Figure 6). This can be explained by the fact that Dongguan is an industrial city, while Shenzhen has a big port with intensive traffic, as it is a coastal city, which increased the LSTs in both cities [66]. The daytime HRI index also increased more in Foshan and Zhongshan than in other cities (Figure 6), which can be due to the impervious surfaces, which obviously increased in these cities from 2005 to 2015 [37]. However, Zhaoqing, Jiangmen, and Huizhou, which were characterized by the lowest urbanization rates, showed the lowest HRI values. These findings can provide guidelines for the improvement of the thermal environment in urban agglomerations.

The drivers of LST were analyzed by applying the GD metric. All of the drivers showed significant influences on daytime and nighttime LST in the PRD. However, the main driving factors varied from city to city. Vegetation has a cooling effect through evapotranspiration and by shading other urban facets. Vegetation can absorb a large fraction of solar radiation and use a large fraction of it for evaporation and transpiration, thereby reducing sensible heat transfer to the boundary layer, i.e., reducing the temperatures of both the vegetation and urban air [19,38]. The shade from vegetation protects other urban surfaces from solar radiation and prevents increases in air and surface temperatures [87]. Our results showed that the NDVI had a strong correlation with daytime LST in Shenzhen and Guangzhou. Our results agree with those of [19,30]. The NDVI also showed a strong influence on nighttime LST in Dongguan, Shenzhen, Foshan, Zhongshan, and Guangzhou. This result is in agreement with the findings of [88].

Variations in LST due to elevation differences should be considered in LST studies in large areas where the terrain is not flat [89]. A negative correlation was observed between LST and elevation [90]. In our study, the elevation factor had a stronger influence on daytime LST in Huizhou, Zhaoqing, Jiangmen, Guangzhou, and Shenzhen. The elevation values of these five cities are higher than those of the other cities in the PRD (Figure 1 and Table 1). The high density of vegetation in the mountainous areas led to a great impact on temperature [91]. In addition, the low urbanization rates in Zhaoqing, Jiangmen, and Huizhou indicated the weak influences of the ISDD, GDP, POP, and NLI factors on daytime LST, which was in line with our results (Figure 7). LST is usually affected by multiple factors, and the higher correlation between the elevation factor and LST in these three cities may be the result of the weak influences of the other factors that were related to the urbanization rate on daytime LST. Similarly, the differences in correlation intensity between elevation and LST in 2005 and 2015 can also be partially explained by the variations in the influences of other factors on LST. The influence of elevation on nighttime LST in the PRD was lower than that on daytime LST in the PRD (Table 6, Figures 7 and 8); this could be due to the location of the PRD. The influence of the ocean at night makes the driving mechanism of nighttime LST more complex. This finding is in agreement with those of [37].

In this study, POP had a significant influence on daytime LST in Dongguan, Shenzhen, Foshan, Zhongshan, and Guangzhou, which were characterized by high urbanization rates (Figure 7). However, a case study in the Yangtze River Delta urban agglomeration found that there was no significant correlation between the surface UHI intensity and population density [92]. The statistical population census data aggregated by the district were used in the study of the Yangtze River Delta, while gridded population data with a 1 km spatial resolution were used in our study, and these were consistent with the spatial resolution of the LST data used in our study. Population indirectly reflects the development of an area and the complexity of the urban surface. An increasing population leads to rapid transformation of natural land cover into impervious surfaces, such as buildings, streets, and other human-made features, which can reduce thermal admittance and increase the surface temperature through a modification of the energy balance [31]. In addition, a high population density can directly (by metabolic heating) or indirectly (by anthropogenic heat emissions) affect urban surface temperature.

In 2005, the socio-economic drivers explained most of the observed variations in daytime LST in Dongguan, Shenzhen, Guangzhou, Foshan, and Zhongshan, which were characterized by high urbanization rate. The daytime LST is affected by landscape changes and anthropogenic activities [9]. Socio-economic activities are usually associated with energy consumption, which increases anthropogenic heat and local surface temperature [93]. The change rates of the population and electricity consumption were positive in these cities (Table 1), which indicated the possible increase in the intensity of the population and socio-economic activities. The NLI represents the intensity of nighttime lighting, and it is an indicator of human density [86]. Anthropogenic heat emissions caused by human activities contribute to daytime LST increases. A low intensity of nighttime lighting is the result of agricultural activities, especially in underdeveloped areas with low urbanization rates. In addition, the brightest areas in nighttime light images could represent factories or industries, which also consume energy and generate heat emissions in the daytime, causing daytime surface temperature changes [94]. Consequently, it can be inferred that the greater influence of NLI on the daytime LST in Dongguan and Shenzhen cities was related to the high urbanization rates and industrialization characteristics of these two cities. Our results agree with those of [95]. The results for the drivers in these nine cities further revealed the mechanisms that drive the LST.

### *5.2. Management of Urban Agglomerations*

The results of the HSAs and HRI showed that the spatial distribution of LST in the PRD urban agglomeration was highly heterogeneous among the nine cities in both daytime and nighttime. Therefore, it is recommended that more attention and resources be committed to improving the thermal environments of cities that are characterized by higher thermal stress in daytime, such as Dongguan, Shenzhen, Zhongshan, and Foshan. Human-made surfaces absorb solar radiation and lead to heat storage in urban areas [39]. As the urbanization is greater in Dongguan and Shenzhen compared to Huizhou, Jiangmen, and Zhaoqing, the thermal stress in daytime is higher in Dongguan and Shenzhen than in Huizhou, Jiangmen, and Zhaoqing. To reduce the thermal stress, the spatial expansion of urban areas should be controlled in these cities. On the other hand, increasing green vegetation areas is an effective way to mitigate the increasing LSTs in regions with high thermal stress. The significant reduction of HSAs in the northern part of the PRD is due to the implementation of greening measures by local governments [20]. Therefore, increasing the green areas may also be an effective way to alleviate the thermal stress in Shenzhen, Dongguan, Guangzhou, and Huizhou, as the NDVI explained over 40% of the variations in the daytime LST in these cities (Figure 7). In addition, the NDVI also showed a strong influence on nighttime LST in Dongguan, Shenzhen, Foshan, Zhongshan, and Guangzhou. The dominant drivers of daytime LST in the central cities of the PRD were the socio-economic ones in 2015, including in Dongguan, Shenzhen, Foshan, Zhongshan, Guangzhou, and Zhuhai. The HRI greatly increased in Dongguan, Shenzhen, Foshan, and Zhongshan, which could be due to the increased intensity of socio-economic activities. Therefore, it is suggested that the heat emissions caused by socio-economic activities should be considered in the urban planning in order to reduce the thermal stress in these cities. The results of this study can provide useful guidance for planners towards better management of the thermal environment.

### *5.3. Limitations and Future Work*

This study provided a comprehensive framework for identifying the spatial distribution characteristics and main driving factors of the LST at different spatial scales in the PRD. There are still some limitations. First, this study analyzed the spatial patterns and the influencing factors of the LST for just two years: 2005 and 2015. Therefore, further studies

are needed in order to address this aspect by applying the most recent multi-temporal remote sensing data in the future. Second, the annual scale was selected in this study, as it provided good results in many previous studies [31,37,96]. Seasonal changes also need to be studied in the future. Third, we considered five levels of LST classification, but more classification criteria need to be considered in future studies. Fourth, the combined effects of two drivers on the LST were measured in this study. However, the combined effects of multiple variables on the LST are important and need to be studied in future works. Finally, the spatial resolution of the remote sensing data must be improved by using higher-resolution remote sensing data in the future.

### **6. Conclusions**

This study took the PRD urban agglomeration, a typical urban agglomeration in Southeast China, as a case study and focused on the spatial distribution and main driving factors of the LST in 2005 and 2015. A spatial autocorrelation analysis, hotspot analysis, and the HRI were used to quantitatively analyze the spatial distribution characteristics of the LST. Seven drivers of LST were studied, including physiographical, land cover, and socio-economic drivers. The GD metric was applied to explore the explanation rate of the driving factors for the LST. The main conclusions from this study can be summarized as follows:

This study is expected to provide useful guidance for the optimization of the thermal environment in PRD urban agglomerations. In addition, the research framework used in this study can also be applied in other urban agglomerations. The results of this study revealed the complex mechanisms of the variability of the LST among the nine cities in the PRD urban agglomeration, and corresponding improvement measures have been recommended to help urban planners in better articulating urban agglomeration management strategies.

**Author Contributions:** Conceptualization, W.L. and Q.M.; methodology, W.L. and M.A.; software, data curation, writing-original draft, W.L.; writing—review, Q.M., L.Z., D.H. and M.M.; writingreview and editing, W.L, M.A. and M.M.; visualization, M.A. and D.H.; investigation, L.Z.; supervision, Q.M.; funding acquisition, Q.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Tsien Hsue-shen Urbanology Award of Hangzhou International Urbanology Research Center & Zhejiang Urban Governance Studies Center, grant number no.20QXS007; the Major Projects of High-resolution Earth Observation Systems of National Science and Technology, grant number 05-Y30B01-9001-19/20-1.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding website.

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

### **Abbreviations**

### **References**

### *Article* **A Remote Sensing Approach for Surface Urban Heat Island Modeling in a Tropical Colombian City Using Regression Analysis and Machine Learning Algorithms**

**Julián Garzón 1,2,\*, Iñigo Molina 1, Jesús Velasco 1 and Andrés Calabia 3**

**Abstract:** The Surface Urban Heat Islands (SUHI) phenomenon has adverse environmental consequences on human activities, biophysical and ecological systems. In this study, Land Surface Temperature (LST) from Landsat and Sentinel-2 satellites is used to investigate the contribution of potential factors that generate the SUHI phenomenon. We employ Principal Component Analysis (PCA) and Multiple Linear Regression (MLR) techniques to model the main temporal and spatial SUHI patterns of Cartago, Colombia, for the period 2001–2020. We test and evaluate the performance of three different emissivity models to retrieve LST. The fractional vegetation cover model using Sentinel-2 data provides the best results with R2 = 0.78, while the ASTER Global Emissivity Dataset v3 and the land surface emissivity model provide R2 = 0.27 and R2 = 0.26, respectively. Our SUHI model reveals that the factors with the highest impact are the Normalized Difference Water Index (NDWI) and the Normalized Difference Build-up Index (NDBI). Furthermore, we incorporate a weighted Naïve Bayes Machine Learning (NBML) algorithm to identify areas prone to extreme temperatures that can be used to define and apply normative actions to mitigate the negative consequences of SUHI. Our NBML approach demonstrates the suitability of the new SUHI model with uncertainty within 95%, against the 88% given by the Support Vector Machine (SVM) approach.

**Keywords:** Surface Urban Heat Island (SUHI); Land Surface Temperature (LST); Principal Component Analysis (PCA); Multiple Linear Regression (MLR); Machine Learning; Naïve Bayes

### **1. Introduction**

Urban expansion transforms natural areas into surfaces covered with concrete, asphalt, and buildings (highly impervious materials), reducing evapotranspiration and decreasing the cooling capacity of the air, which in turn helps to reduce the impacts of high urban surface temperature on the urban surface. Due to the existing urban growth, the climate in these areas becomes warmer than the regional areas of the suburban and rural regions, resulting in the phenomenon of Urban Heat Islands (UHI) [1]. The UHI refers to a phenomenon in which urban areas tend to have higher air or surface temperatures than their surroundings [2]. Traditionally, terrestrial observation methods, such as ground meteorological stations that record specific values of air temperature, have been used to model UHI [3]. The difference between air temperature measurements recovered from urban and rural meteorological stations is a direct method used to model UHI [4]. However, the high heterogeneity in urban areas makes temperature spatially diverse, making it difficult for a small number of stations to realistically represent the real variability [5]. When the UHI phenomenon is monitored by remote sensing, it is referred to as Surface Urban Heat Island (SUHI). The reason is that the parameter considered here is the Land

**Citation:** Garzón, J.; Molina, I.; Velasco, J.; Calabia, A. A Remote Sensing Approach for Surface Urban Heat Island Modeling in a Tropical Colombian City Using Regression Analysis and Machine Learning Algorithms. *Remote Sens.* **2021**, *13*, 4256. https://doi.org/10.3390/ rs13214256

Academic Editors: Yuji Murayama and Ruci Wang

Received: 16 September 2021 Accepted: 19 October 2021 Published: 22 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Surface Temperature (LST), which differs from studies of air temperature [6]. Therefore, LST is an essential variable to characterize SUHI, which has been listed as an essential climate variable of the World Meteorological Organization program. This variable is an important indicator of the energy balance between the atmosphere and the surface of the Earth [7]. Zhou et al. [8] presented a broad review of the SUHI phenomenon and suggested several methods of analysis through the integration of remote sensing data, thermal trends, field observations, and numerical modeling. Sekerteking and Zadbagher [9] suggest that to model and simulate LST, it would be important to investigate the performance of various Machine Learning methods associated with statistical and numerical models. Li et al. [10] and Mei et al. [11] pointed out that one of the problems to be solved in inference models of LST or geographic variables is to evaluate the influence of contributing factors.

Land-use land-cover changes (LULC) and spectral indices from satellite data, such as, e.g., Normalized Difference Vegetation Index (NDVI) or Normalized Difference Water Index (NDWI), have been used extensively to investigate the relationships between urban and biophysical systems, as well as their impact on surface temperature [12–15]. LULCs are due to complex interactions between the urban system and the biophysical environment that produce significant changes in local temperatures. In a recent work published by Shi et al. [16], eight parameters were referred to as urban design factors, where thermal properties of building materials, vegetation, vegetation cover ratio, and ground emissivity were taken into account. The authors suggested the need to detect more potential factors affecting this phenomenon.

Understanding and quantifying urban temperatures in space and time are significantly relevant for city planners in defining policies that generate adaptation strategies to mitigate the SUHI effects. A very useful tool is Principal Component Analysis (PCA), a multivariate statistical technique that aims to preserve the total variance and reduce the dimensionality of the data set, while eliminating redundancy in the data [17]. Several authors have used this technique to detect spatial patterns of biophysical factors by synthesizing information from a set of images [18–20]. Multiple Linear Regression (MLR) analysis is an approach used to evaluate the relationship between independent and dependent factors [21]. This method has also been widely used to determine the relationship of various environmental factors [22–24].

Advanced nonlinear analysis techniques, such as Machine Learning, have been applied in numerous studies that require analyzing variables related to urban thermal changes. Some examples are population density, land cover, and urbanization [25–29]. Voelkel and Shandas [30] implemented a UHI model to detect a daily distribution of temperatures. Their results revealed that a random forest (RF) model performed better in predicting temperature. Furthermore, Zumwald et al. [31] developed a model to create high-resolution air temperature maps. This model makes predictions by integrating an RF algorithm with low-cost weather stations. It is important to note that the behavior of SUHI varies over time and is associated with factors such as human development and changes in land use. In this sense, Kafy et al. [32] formulated a seasonal thermal prediction influenced by LULC through Cellular Automata and Artificial Neural Network algorithms. Their findings indicated that by 2039, the urban growth of Cumilla, and Bangladesh, plus the decrease in land cover, will cause 30% of cities to experience temperatures above 33 ◦C. Shi et al. [33] noted that the use of time series and Machine Learning techniques is a growing trend in SUHI research.

Several studies have widely documented the influence of spectral indices such as the NDBI, NDVI and Normalized Difference Water Index (NDWI) on the SUHI phenomenon [34,35]. However, none of these authors considered the weighted contribution of the factors to temperature changes, while it is an interesting analysis that can identify which factors generate the greatest influence on SUHI. On the basis of these factors, specific adaptation measures to thermal change can be defined and applied. The main novelty of the approach proposed in this work is the application of a weighted Naïve Bayes

Machine Learning (NBML) algorithm to segment the geographical space into regions of different thermal intensity, not explored in previous literature.

Understanding and quantifying urban temperatures in space and time is significantly relevant for city planners in defining policies that generate adaptation strategies in the face of adverse effects of SUHI. Here we study the application and assessment of modeling procedures that allow evaluating the contribution of various factors to SUHI. For this purpose, a combination of PCA and MLR techniques applied along with Machine Learning Algorithms is used to detect high thermal intensity patterns in the tropical Colombian Andean city of Cartago. Although SUHI is a derived quantity, expressed as the difference between urban and rural LST, the delimitation of thermal zones using LST ranges allows establishing comparisons with other zones, e.g., rural areas, and classifies the space into zones with greater or lesser thermal activity. In this study, the LST ranges are taken from Wang et al. [36], as they are based on statistical criteria, and they appear to conveniently reflect the LST differences of urban areas with their surroundings. The spatial patterns of the SUHI phenomenon can be represented through LST ranges, which, combined with the weights of the involved variables, are further classified using Machine Learning algorithms.

The methodology suggested in this article establishes an effective method for assessing SUHI patterns, locally, and attempts to draw several recommendations for planning sustainable urban development and for the regeneration of areas with thermal excesses.

### **2. Materials and Methods**

### *2.1. Study Area*

The city of Cartago is located in the south-west area of Colombia in the Andean region at an altitude of 917 m above sea level. It has an extension of 279 km2 with moderate topographic relief. The geodetic coordinates of the city center are 4.75◦ N and 75.9◦ W. This area belongs to Valle del Cauca Department and is surrounded by the Cauca and La Vieja rivers.