The Gradient Effect on the Relationship between the Underlying Factor and Land Surface Temperature in Large Urbanized Region

: Although research relating to the urban heat island (UHI) phenomenon has been signif-icantly increasing in recent years, there is still a lack of a continuous and clear recognition of the potential gradient effect on the UHI—landscape relationship within large urbanized regions. In this study, we chose the Beijing-Tianjin-Hebei (BTH) region, which is a large scaled urban agglomeration in China, as the case study area. We examined the causal relationship between the LST variation and underlying surface characteristics using multi-temporal land cover and summer average land surface temperature (LST) data as the analyzed variables. This study then further discussed the modeling performance when quantifying their relationship from a spatial gradient perspective (the grid size ranged from 6 to 24 km), by comparing the ordinary least squares (OLS) and geographically weighted regression (GWR) methods. The results indicate that: (1) both the OLS and GWR analysis conﬁrmed that the composition of built-up land contributes as an essential factor that is responsible for the UHI phenomenon in a large urban agglomeration region; (2) for the OLS, the modeled relationship between the LST and its drive factor showed a signiﬁcant spatial gradient effect, changing with different spatial analysis grids; and, (3) in contrast, using the GWR model revealed a considerably robust and better performance for accommodating the spatial non-stationarity with a lower scale dependence than that of the OLS model. This study highlights the signiﬁcant spatial heterogeneity that is related to the UHI effect in large-extent urban agglomeration areas, and it suggests that the potential gradient effect and uncertainty induced by different spatial scale and methodology usage should be considered when modeling the UHI effect with urbanization. This would supplement current UHI study and be beneﬁcial for deepening the cognition and enlightenment of landscape planning for UHI regulation.


Introduction
As one of the most evident features of humans' impact on the Earth's system, the urban heat island (UHI) phenomenon has attracted considerable attention in almost all major cities worldwide [1,2]. It is worth noting that the UHI effect has been exacerbated along with ongoing urban expansion and constant population influx in recent years, presenting substantial adverse threats to both residential and environmental health in and around cities [3][4][5]. When considering this, scholars across various disciplines have not stopped discussing the topic of UHI, in order to develop deeper understandings on the coupling relationship between UHI and the urbanization process [6,7].
As one of the many factors that influence the urban thermal environment, the landscape features of the underlying surface represent the most commonly used driving factor for examining the quantization relationship with UHI in most of the related studies [8]. This requires the conversion of these environmental factors into quantitative indicators in advance for statistical causal analysis. For UHI, previous studies have tended to use the difference in the surface air temperature (between the urban core and the suburbs) in order to characterize the UHI intensity [9]; with the development of Earth observation technology, the remotely sensed land surface temperature (LST) has been widely used as a more geo-convenient way to record the heterogeneity of the land surface thermal environment [1,10]. For investigating the landscape features of the underlying surface, researchers can easily access the detailed land cover and change information for most parts of the Earth's system while using modern remote sensing and geographic information technology. These land surface details can be further quantized into numerical landscape metrics following the landscape ecology theory, and they then fall into two categories, as landscape composition and configuration characteristics [11]. Following this, the spatial and statistical analysis can be employed in order to establish the statistical relationship between these numerical variables, i.e., the dependent (LST or UHI intensity) and independent (landscape indexes) variables.
A previous publication [6] has provided a systematic review of the statistical methodologies that were used in UHI related research, stating that the conventional statistics (e.g., ordinary least squares (OLS)) are the primary (>60%) tools that are employed by researchers in current studies investigating the impacts of underlying factors on UHI variation. However, most of these studies have been conducted within a single city, in order to better serve the planning task of a particular urban region [3,12]. However, multiple cities in economically developed regions gradually tend to exhibit centralization with urbanization, and then form urban agglomeration regions. This regional urbanization process will also have a profound impact on the land surface heat redistribution, so it has attracted the attention of many scholars [13,14]. At this point, there is still a lack of research on the broader metropolitan areas of multiple cities. In particular, in these regions, the environmental factors tend to be more complex and heterogeneous over time and space [15]. This presents a challenge to using conventional statistics in large-scale metropolitan areas because of the spatial non-stationarity, which refers to the more spatially varying relationships between the dependent and independent variables [12].
Moreover, UHI and spatial data both have multiscale features, which have been verified to change with different spatial analysis units [16,17]. Generally, the fine-scaled spatial unit can help to better reveal the spatial details of the study area. However, for research areas with various spatial sizes, the scale issue of the spatial analysis unit should be considered based on many external factors, e.g., spatial resolutions of both the land cover and LST data [18], or the modeling efficiency (a large study area with a finer unit means a higher sample size and computational load and vice versa) [17]. For this reason, the use of the spatial analysis unit in previous studies has not presented a unified methodology as a basis, and it has varied with different forms and sizes [19]. This context sensitivity for both UHI and land surface characteristics may further exacerbate the uncertainty that is posed by spatial-statistical analysis, especially for conventional (global) analysis methods that do not have good spatial suitability.
In general, there is still a lack of a continuous and clear recognition of the change pattern of the UHI-landscape relationship along different spatial gradients in current research. Thus, studying the potential dynamic between the heterogeneous LST and its contributory factors at different spatial scales is essential for providing insights into the dynamic mechanism of UHI [20]. Furthermore, few studies have employed comparative analysis in order to discuss the suitability of different statistical methods in multiscale UHI research [6]. When considering this, we chose a large-scaled urban metropolitan region in China as the case study area, in order to include more spatially complex and heterogeneous LST and land surface underlying characteristics. Subsequently, a comparative analysis that was associated with global-and local-based spatial-statistical models was applied in order to evaluate the validity of the results and robustly infer the association of UHI with its driving factors across different sizes of grid units. The results of our study can enhance the understanding of the potential gradient effect on the thermal influence by urbanization.

Study Area
The Beijing-Tianjin-Hebei (BTH) region was selected as the study area, and it is located in the North China plain (36 • , covering an area of nearly 2.16 ×10 5 km 2 . This region is dominated by a typical temperate monsoon climate, with an annual mean temperature of 12 • C and mean precipitation of approximately 500 mm [21]. Most of the precipitation and hottest temperatures occur in the summer season (June, July, and August), and the winter season is usually cold and dry [14]. The study area includes two municipalities (i.e., Beijing and Tianjin) and 11 prefectural cities ( Figure 1). As one of the largest socio-economic urban agglomerations in China, the BTH region has experienced a rapid urbanization process in the past few decades. Until the year of 2015, the population urbanization rate of this region reached >60%, and 69.7 million people lived in the urban areas (7509.63 km 2 ) [22]. However, dramatic changes of land cover by urbanization have caused serious environmental degradation in this area and they have posed a substantial health threat to the densely populated residents [23].  An increasingly prominent UHI phenomenon is one of these environmental threats, which has been frequently discussed in the BTH region [13,24]. This provided us with the initial impetus to choose this area as an ideal case. Moreover, the BTH region presents additional research advantages for this study. One is that, although the BTH region has experienced significant changes in urban forms, the 13 cities that are within this region are experiencing different levels of urban development [22]. For example, there was a nearly 20-fold GDP gap between the city of Beijing (22,968 RMB) and Hengshui (1220 RMB) in 2015 [21]. Farmland and forest are the dominant land cover types in this region (mainly distributed in the mountainous area), and the urban areas are mainly surrounded by farmland [13]. Various levels of urbanization have formed the underlying landscape with clear spatial heterogeneity at the regional scale, thus providing us with diverse landscape characteristics for discussing our study topic.

Land Cover and Land Surface Temperature Data
There are two types of data source used in this article: land cover and LST data. The land cover details of the BTH region were identified based on remote sensing technology. The remote sensing images (cloud-free, <10%) were acquired from the Geospatial Data Cloud (http://www.gscloud.cn/), including both Landsat-7 (ETM+, with a spatial resolution of 30 m) and -8 (OLI, 30 m) data. The acquisition time was during the summer (June to August) in 2000, 2005, 2010, and 2015, respectively. Subsequently, the land cover information was obtained using the backdating interpretation approach (taking the classification result of 2015 as the benchmark reference), for which the detailed procedure can be referred to in [25,26]. A total of six types of land cover were eventually divided, as forest, grassland, wetland, farmland, built-up land, and bareland, respectively ( Figure 2 and Table 1). The accuracy of all the interpretation results was verified >95% [25], in order to fulfill the usage purpose of our study.
The LST data of the BTH region were obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) daytime datasets. First, the Terra MODIS LST 8-day composite products (observed at 10:30 a.m. local time, with a spatial resolution of 1 km, MOD11A2) that were recorded during the summer season between 2000 and 2015 were downloaded from NASA (http://earthdata.nasa.gov/), and they were retrieved with a generalized split-window algorithm that generally has an absolute bias of <1 K [27,28]. Second, the LST values ( • C) were converted by the surface heat radiation intensity (i.e., the pixel value) of the MOD11A2 datasets by the following conversion formula: where DN is the digital number of MOD11A2 datasets, 0.02 is the radiation scaling ratio (from metadata of the datasets), and 273.15 is the conversion gap between the Kelvin (K) and Celsius temperature ( • C). Afterwards, we further composited these LST data (eightday composite) from June to August into single data in order to represent the average LST in the summer seasons of 2000, 2005, 2010, and 2015, respectively ( Figure 3). Thus, we adopted these data to conduct the following analysis, by referring to the conclusion of previous studies that the strongest UHI intensity is usually found during the summer season in northern China [29][30][31].

Grid Analysis Design
The theory of landscape ecology suggests that multiscale information is required to obtain a better understanding of the relationship between the regional landscape structure and ecological process in the context of spatial heterogeneity [11,32]. This is because both of the elements may have different performances at different spatial and temporal scales, thus leading to complexity in geological/ecological research [33,34]. Therefore, discussing the gradient effect will be beneficial in revealing the potential rules of specific ecological processes.
In research that is related to the urban thermal environment, issues of the gradient/scale effect are often raised, as we summarized in the introduction section. However, the proper analytical unit and its size are still being debated in the literature [35][36][37]. In this study, we thus applied a multiscale polygon grid-based analysis in order to examine the strength of the relationship between the land cover composition and LST. This approach looks into the statistical relationship between land cover density (LCD) and mean LST, but no proper grid sizes have been suggested by previous studies [38]. Therefore, we simply set up a gradient of different sizes of grid for subsequent comparative analysis, taking 6, 9, 12, 15, 18, 21, and 24 km as the analysis scales, respectively. All of the grid sizes were designed to be an integer multiple of the pixel size of the thematic data (1 km for LST and 30 m for land cover), in order to ensure that the pixel values of these thematic data within each grid could be effectively summarized.
Subsequently, using spatial overlay analysis, we summarized the two types of thematic information for each grid, i.e., land cover ( Figure 2) and LST ( Figure 3).  Specifically, for the land cover data, the density of each type of land cover (%) was calculated, as follows: where D is the LCD, i = [forest, grassland, wetland, farmland, built-up land, and bareland represents the type of land cover, j = [6,9,12,15,18,21, and 24 km] represents the type of grid unit, n is the pixel number for land cover i, 9 × 10 −4 is the fixed area for each land cover pixel, and A j is the total area for the j type of grid. For the LST data, the zonal statistics tool in the ArcGIS platform was used to aggregate the mean LST value ( • C) for each grid, representing the thermal environment characteristics. Additionally, all of the (intact) grids were constrained within the study area, in order to obtain complete and correct land cover and LST information. Grids that were located at the edge of the study area were deleted, because of their LCD (the total land cover area < the area of the grid) or LST value (null value). Figure 4 shows a sample graph of the thematic characteristics for different scales of grid.

Statistical Analysis
Previous studies have documented that land use and cover change (LUCC) have been treated as the major driving determinants for the alteration of the surface thermal environment, especially referring to the conversion between natural and artificial construction land [5,14,39,40]. In addition, other types of variables, such as landscape pattern characteristics, surface biophysical parameters, socio-economic factors, etc., have also been demonstrated to have significant effects on LST [36,37,41]. As stated, the main purpose of this study is to discuss the potential gradient effect in the relationship between the landscape features of the underlying surface and LST in a large region. Under the a priori expectation that LUCC influences LST and, for the sake of simplicity, only LCD was used as the independent variable, while the summer daytime LST was treated as the dependent variable. Both types of variable were processed, as described in Section 2.3 ( Figure 4).
The Spearman correlation matrix was first developed to reveal the bivariate association strength between LST and the density of land cover across the grid sizes and time periods [19,42]. This can help to preliminarily determine which type of land cover has the most effective contribution to LST.
Subsequently, based on the results of correlation analysis, we conducted further regression analysis in order to elaborate on the deeper quantitative relationship between the key land cover element and summer daytime LST. The review by [6] concluded that most relevant studies (>68%) tend to use global-based analysis to quantify the relationship between the spatial variables and UHI variables, of which ordinary least squares (OLS) regression is the dominant method among the reviewed studies. However, when considering the topic of UHI at a regional scale, the UHI effect may show context sensitivity and vary significantly over time and space [20,43]. Therefore, it is suggested that the UHI-related phenomenon should be locally modeled for an entire study area. In this study, accordingly, both global-based (OLS) and local-based (geographically weighted regression, GWR) analysis were conducted in order to explore the relationship between LST and the land use density in the BTH region.
OLS regression is a widely used model for quantifying the dependent variable in terms of its relationship with explanatory variable [44,45], with the following equation: where Y denotes the dependent variable (i.e., the summer LST), X is the explanatory variable (i.e., the built-up density), β is the slope coefficient, β 0 is the intercept value, and ε is a random error term.
The GWR method is one of several spatial regression approaches increasingly used in geography and other disciplines [12,46]. GWR provides a local model for incorporating the spatial characteristics of both dependent and explanatory variables and creating conditions for analyzing the spatial characteristics of the regression relationship, expressed as where k indicates the kth grid unit, (x k , y k ) represents the spatial coordinates for the grid unit k, and the meanings of Y, X, β 0 , β and ε are the same as those of OLS regression. In contrast to the OLS model, the estimation for Y k in the GWR model is mainly determined using limited independent variables X near the location k. Each of the independent variables has a specific contribution weight for Y k , while that in the OLS model has a weight of unity. A Gaussian distance-decay-based kernel function is commonly used in order to calculate the weight (W) for GWR model while using the following function: where d indicates the Euclidean distance between grid unit k and its neighboring grid unit, and b is referred to as the bandwidth parameter. Based on this, the GWR tool will automatically find the optimal distance (for a fixed kernel) by using the Akaike Information Criterion (AICc)-based bandwidth method. More details regarding the two regression methods can be found in [44][45][46]. Finally, the performances of the OLS and GWR regression models were compared across different grid sizes and time periods. The adjusted coefficient of determination (R 2 ), the global Moran's I of the standard residuals, and AICc were used in order to evaluate the performances of the built models with respect to the goodness-of-fit and residual spatial autocorrelation. According to [46], the (adjusted) R 2 can be interpreted as the proportion of the variance of the dependent variable that is covered by the regression model, which can be used to test the fitness of the regression model (ranging from 0.0 to 1.0, where a larger value indicates a better model fitness). The AICc can be used to test the model performance and compare the regression models. When considering the complexity of the model, the model with a lower AICc value will better fit the observed data. Comparing the AICc value of OLS and GWR models is one way to assess their advantages. The I index is used for testing the spatial distribution characteristics of the model standard residuals. If the I index loses its statistical significance (p > 0.05), then this indicates that the standard residuals are spatially randomly distributed. These diagnostic factors can also be spatialized (by local Moran's I analysis, etc.) in order to reveal potential findings in depth.
In this study, the statistical analysis of attribute data, e.g., the Spearman correlation analysis, was conducted while using the statistical software of OriginPro. The spatial analysis procedures, such as OLS and GWR regression analysis, were completed using the platform of ArcGIS Desktop [46].

Results
The Spearman correlation analysis showed that all of the LCDs were significantly related to LST across different grid sizes and periods, except for bareland (in some cases), as shown in Table 2. Among these, the density of built-up land displayed the most significant and positive correlation with the LST variation in all cases, while that of the forest showed the strongest negative correlation. This indicates that the built-up land and forest mainly affected the LST variations in the BTH region. Therefore, we will focus on the performance of this positive driver for the UHI phenomenon in the subsequent analysis, i.e., the built-up land density. Table 2. The Spearman correlation coefficients for the land cover density and land surface temperature (LST) across different grid sizes and periods (** p < 0.01 level and * p < 0.05 level). The global regression models estimated statistically significant (p < 0.01) relationships for the built-up density and LTS for all cases, as shown in Table 3. Generally, from the year 2000 to 2015, the global R 2 for these models exhibited an increasing trend. Moreover, in each year, the higher R 2 showed that the OLS models have a better measurement of goodness of fit as the grid scale increases. However, the standard residuals for all of the OLS models displayed significant spatial clustering (i.e., spatial correlation) characteristics (I, p < 0.01), which indicated that the OLS model may not be suitable for quantifying the relationship between the built-up density and LTS in a relatively large region. Table 3. Comparison of ordinary least squares (OLS) regression and geographically weighted regression (GWR) for the built-up density and LST across different grid sizes and periods: Summary of the coefficients and diagnostics. R 2 is the adjust coefficient of determination, AICc is the corrected Akaike Information Criterion, and I is the global Moran's I of the standard residuals (** p < 0.01 level). In contrast, the GWR models seem to be better suited to the investigation of the region thermal environment and the underlying influencing factors (Table 3). For each case, the GWR model exhibited an apparently higher R 2 value than that of the OLS model. This is further confirmed by the AICc of each pair of models. If the difference between the AICc values of the two models is >3, then the model with a lower AICc value can be considered to be the better model. When considering the complexity of the GWR regression algorithm, these results demonstrate that the GWR models function better than the OLS models. Additionally, the GWR models present a more stable performance than the OLS models across different grid sizes and periods. In particular, the significance of these I values for most of the cases (p > 0.05) further demonstrated that the prediction results that are produced by GWR models are more applicable to the situation of this study.

The Relationships between LST and the Built-Up Density Revealed by the OLS and GWR Models
This study discussed the potential relationship between LST and underlying physical factors across various spatial scales. The OLS and GWR analyses both confirmed that the underlying physical factors exert a significant impact on the LST variation, as shown in Table 3.
In terms of the spatial-temporal scale, as the most frequently used regression method, the increasing R 2 of OLS models from 2000 to 2015 indicated that the urban sprawl of the study region further strengthens its influence on the UHI phenomenon. Moreover, notably, various R 2 values for different spatial analysis units implied that a spatial scale issue does exist in regional UHI studies. Similar findings have been discussed by [17], which recorded a lower strength of linear regressions (R 2 ) between urban ecological land and LST with an increasing grid size (spatial analysis unit) in the city of Beijing. This gradient effect can be mainly attributed to the changes in sample statistic responses to different grain sizes. In other words, most of the studies that focused on the relationship between LST and land surface coverage were inclined to use the average (or other types of statistical indicators) values as the statistical variables that are based on the specific spatial analysis unit [8,18,36,38]. The sample size and range (max-min) both narrowed with larger grid sizes, leading to a reduction in the complexity and heterogeneity of the variables for regression analysis, as shown in Figure 5. This generalization with a larger spatial analysis unit is very much in line with the global regression models, i.e., the OLS model [12].    However, the performance of OLS models was much lower than that of the GWR models. The main reason for this is that, as a classic regression approach, OLS analysis tends to provide an estimate of the general situation for an entire region, regardless of the spatial pattern [47]. In this study, the underlying physical factors and LST in the region of BTH both exhibited significant spatial heterogeneity, as shown in Figures 2 and 3. This will weaken the quantitative expression of the relationship between the two variables in this study, i.e., built-up density and LST, by classic regression methods. Evidence for this can be found in relevant studies. For example, while using global analysis methods, [19,48] found an R 2 value of >0.7 for quantifying the relationship between urban built-up land and LST in a highly developed urban region. In contrast, [13,37] found an R 2 of nearly 0.5 in order to quantify the OLS regression models for urban built-up land and LST in a region with urban and suburban areas. Subsequently, a further cut in R 2 was found in this study with a larger region including urban, rural, and natural areas. Moreover, although the OLS regression results are statistically significant (p < 0.01, Table 3), the I index with the same statistical significance (p < 0.01, Table 3) further indicates that the OLS model has certain limitations in this study. Because of that, ideal spatial regression prediction results require both the high and low predicted (standard) residuals that are to be randomly distributed [46,49]. Furthermore, the local Moran s I analysis results highlight the obvious spatial clustering characteristics of the standard residual by the OLS models, as shown in Figure 6.
In view of the above, the UHI phenomenon exhibits high spatial heterogeneity for characterizing the thermal content of the underlying surface, which is difficult to characterize with global regression methods. Instead, the GWR model exhibits a better performance than the OSL model for depicting the spatial non-stationarity relationship between LST and underlying variables, with a higher R 2 (>0.79) and lower AICc (Table 3). Moreover, for most of the cases, the insignificant Moran index (I) indicated that the residuals of the GWR models were randomly distributed, in order to meet the verification of the regression results. Additionally, the local Moran s I analysis showed that the spatial cluster characteristics of the predicting residuals for the GWR model are weaker than those of the OLS model in terms of both breadth and intensity ( Figure 6). Specifically, for the OLS model, the cluster areas with high predicted deviation have no obvious spatial regularity, being distributed in all kinds of regions, including urban, rural, and natural regions. On the contrary, the cluster regions for the GWR model are mainly distributed in the mountainous areas of the study area. This is understandable, because the global and local regression models have both been argued to provide unexpected results for elevation [12,29,38], and our results also confirmed that the prediction error of the regression model will increase with a higher elevation (Figure 7).

Implications and Limitations of this Study
This study suggests a significant relationship between the LST and built-up density in the BTH region. However, most previous studies have tended to adopt global-based methods to derive the spatial relationship between the spatial variables and UHI phenomenon, as we mentioned above [20]. This may lead relevant studies to overlook the locally detailed differentiation of the underlying mechanisms in large-extent urbanized regions. In addition, the OLS model will further face prediction uncertainties that are caused by various spatial scales, as shown in Table 3. However, for this, there is still no final conclusion on the so-called optimal spatial scale or unit for UHI study, which is still within the limits of the specific case study [8,18]. However, what is certain is that the LST and its landscape drivers both present a significant spatial non-stationarity, which indicates that the UHI effect is context-sensitive over space. Therefore, the GWR model seems to outperform an aggregated model in order to depict this non-stationarity with multiscale applicability, by inferring a spatial dynamic approach to parameter estimations. When considering the heterogeneity background that was presented in this study (Figures 2 and 3), to some extent, the GWR analysis can provide more robust information for understanding the interrelation between UHI spatiotemporal variation and landscape dynamics across multiple scales.
However, it is worth noting that this does not mean that the GWR model can be treated as a "perfect" model for subsequent UHI-related studies. Our results also showed that, in the cases with smaller grid sizes, the regression performances of the GWR model are not ideal (i.e., statistically significant I in Table 3). This indicates that GWR also has a certain limitation that is related to the spatial scale, although this limitation is more lenient than that of the OLS method. For studies that focus on a smaller urbanized region, such as an individual city or the core built-up region of a metropolitan area, the heterogeneity and complexity used to describe the land surface thermal environment and landscape characteristics will both be further enhanced, due to the usage of more sophisticated spatial and climate data [18,50]. This study adopted seasonal average LST data to avoid the issue of missing data, due to cloudy weather conditions in summer [51,52]. However, this generalization means that our results ignored many research details that are required in studies that are based on air temperature (with a finer temporal resolution) and medium/high resolution remote sensing data (with a finer spatial resolution) [18,53]. In addition, it did not involve the usage of various spatial analysis units, such as a regular grid [47,54], concentric ring [38,55], spatial buffer [13,56], or irregular urban block [37,48], etc. All of these factors will further bring uncertainties for quantifying the relationship between UHI and underlying variables [36]. Beyond that, this study only used the built-up density as the key independent variable for explaining the variation of LST. It has been frequently reported that there are many other types of key factors that have significant impacts on the urban land surface thermal environment, e.g., the landscape pattern index and land surface biophysical parameters [1,6]. Many of these factors have been shown to exhibit various responses to changing spatial analysis scales, which may lead to scaledependent uncertainty in subsequent regression analysis [16]. For addressing these issues, further in-depth spatial-statistical discussions are needed to cover the potential changes in thermal landscape patterns at different grain sizes and in different regional environments, as well as the multi-factor responses, in order to obtain a better understanding of the UHI phenomenon and its associated driving factors across the metropolitan area.

Conclusions
In this study, we chose a large-extent urban agglomeration-the BTH region in Chinaas the study area. While using multi-temporal land cover and LST data based on remote sensing technology, we examined how the underlying surface characteristics affect the UHI effect in the study area. On this basis, this study further discussed the modeling performance in quantifying the relationship between LST and built-up density across various spatial analysis gradients, by comparing the OLS (global) and GWR (local) methods. The main findings can be summarized, as follows.
The OLS and GWR analysis both confirmed that the composition of built-up land contributes as an essential factor that is responsible for the UHI phenomenon in a large urban agglomeration region. The difference is that the OLS model results exhibited a significant (spatial) gradient effect, namely its performance in quantifying the relationship between LST and its driving factors, which changed with different grid sizes (that ranged from 6 to 24 km). This also indicated that using a global-based model may generalize and ignore the heterogeneous characteristics of LST and the landscape in a large urbanized region. In contrast, the GWR model revealed considerably stronger relationships in interpreting these driving forces for LST. Moreover, the results showed a robust and lower scale dependence performance for accommodating the spatial non-stationarity in a UHI study than that of the OLS model. In general, our study acts as a supplement to the field that is related to the scale and methodology issues in UHI research. The main findings provide insights into the dynamic causal mechanism between the UHI effect and landscape factors in large extent urban agglomeration areas and, thus, provide useful implications for future UHI studies and site-specific landscape planning.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the privacy agreement among co-authors.