Spatio−Temporal Changes and Key Driving Factors of Urban Green Space Configuration on Land Surface Temperature

Abstract

Changes in land cover by rapid urbanization have diminished the cooling effect of urban green spaces (UGS), exacerbating the upward trend of land surface temperature (LST). A thorough and precise understanding of the spatio-temporal characteristics of UGS and LST is essential for mitigating localized high temperatures in cities. This study identified the spatio-temporal changes in UGS configuration and LST in Shanghai from 2003 to 2022. The correlation between UGS configuration and LST was explored using spatial autocorrelation analysis and causal inference. The results show that (1) the high-temperature space had grown from 721 km² in 2003 to 3059 km² in 2022; (2) in suburbs, the largest area of UGS tended to decrease, while the number of patches tended to increase, indicating a distinct feature of suburbanization; (3) changes in the largest area of UGS had more significant spatial correlation, indicating that urban sprawl primarily impacts large UGSs; and (4) compared to the number and shape of UGS, changes in the largest area are the key factor influencing regional LST. These findings enrich the knowledge of the spatio−temporal relationship between the UGS configuration and its cooling effect in urbanization, offering valuable insights for building cooler cities.

1. Introduction

In recent decades, the profound impacts of urbanization on sustainable development issues have emerged as global concerns [1,2]. These impacts show more significant hazards in cities with high population densities [3]. The urban heat island (UHI) effect, resulting from extensive conversion of natural surfaces to impervious ones, poses a significant threat to the resilience of urban ecosystems [4], as well as to human health [5,6]. Therefore, how to cool down the city and mitigate the UHI effect has become a hot research topic.

Mitigating the summer UHI effect has become increasingly challenging, as a significant number of inhabitants are likely to endure both intense solar radiation and heat build-ups [7]. To build a comfortable living environment, adjusting the pattern of urban green space (UGS) within the city to mitigate the UHI effect has received focused attention from scholars around the world [8]. There is unanimous consensus among scholars that land cover types such as parks, water bodies, and green roofs play a significant role in reducing the land surface temperature (LST) [9]. Compared to other types, UGS do not only remove from the atmosphere large amounts of carbon dioxide released by human activities to mitigate urban warming [10]. Trees also influence the urban climate through shading, transpiration, and facilitating airflow and heat exchange [11,12]. These influences are closely tied to the spatial pattern of UGS within the city.

Composition and configuration are two important dimensions that impact the cooling effect of UGS [13]. Composition refers to the percentage of UGS area, with differences mainly in the number of trees and the shading area [14]. Numerous studies have examined differences in the cooling effect of UGS composition during the diurnal hours. For instance, research indicates that trees exhibit a stronger cooling effect during pre-dawn and evening, attributed to heightened shading [15]. The magnitude and intensity of the diurnal UHI effect are governed by variances in plant transpiration rates and leaf area [16]. However, fewer studies have focused on the relationship between UGS configuration and LST. Furthermore, there is some controversy in the research on the cooling effect of UGS configuration [17]. One study has posited that the spatial distribution of UGS holds greater significance than the shape of individual patches [18]. Peng, J. et al. found that the morphology of blue-green spaces is more influential than spatial distribution in the cooling effect [19]. Thus, it is important to study the influences of UGS configuration on LST variations.

The UGS configuration tends to change strongly during rapid urbanization [20]. UGS witness a predominant increase in number during urban sprawl growth, while UGS patches tend to expand in size during the compact development phase [21,22]. These UGS configuration changes may have different effects on the LST [23]. Meanwhile, there has been a limited number of studies investigating the change of UGS configuration impact on LST from an inter-annual perspective. This leads to the fact that the strategy of urban planning to improve the cooling effect through the GI configuration is still inadequate [24]. Thus, with the urban sprawl, it becomes imperative to explore the spatial correlation between inter-annual variations in the UGS configuration and LST, identifying UGS features with significant cooling effects.

This study aims to reveal the changing characteristics of the UGS configuration and LST in a rapidly urbanizing area, and the spatial correlations between them. Based on the results of existing studies, this study makes the hypothesis that changes in UGS configuration will produce a significant effect on land surface temperature. The key issues to be clarified are (1) revealing the trends in the UGS configuration and LST; (2) clarifying the spatial differences in the UGS configuration on nearby LST; and (3) seeking the UGS configuration factor that has a significant influence with LST changes. Answering these questions helps to understand the impact of UGS configuration on the LST from a spatio−temporal change perspective. It will help decision-makers to develop measures for sustainable development and UHI mitigation in Shanghai.

2. Materials and Methods

2.1. Study Area

Shanghai (120°52′ E–122°12′ E, 30°40′ N–31°53′ N) experiences a subtropical monsoon climate characterized by high temperatures and rainy summers [25]. As of the end of 2022, the resident population of Shanghai had reached 24.76 million people (https://tjj.sh.gov.cn/, accessed on 10 March 2024). As one of the most rapidly urbanizing and economically developing cities in China, it has experienced massive urbanization since 2003. The urbanization rate has surged from an initial 75.52% to 89.30% by the end of 2022. The study area encompasses the urban extent of the mainland region, spanning an area of 5110.7 km². Human activities in this region are highly intensive and concentrated, contributing to a highly heterogeneous urban green space (UGS) pattern [26] (Figure 1).

Due to the concentration of population and industry in the short term, the UHI effect in the study area becoming more significant, with frequent and long-lasting heat waves [27]. Since 1982, global mean temperature has risen about 0.72 °C and the average rising rate was 0.072 °C/10 years [28], while the difference between urban and suburban meteorological station monitoring data in the study area has shown an upward trend, with a sudden increase of 0.1 °C every 10–13 years [29]. It indicates that warming in the study area has been occurring at a faster rate than globally. Additionally, in the summer of 2022, Shanghai experienced its most severe heatwave event since 1961 [30].

2.2. Data Source and Preprocessing

To describe the significant changes that have occurred at UGS configuration, the rapid development phase of the study area was selected: from 2003 to 2022. In this study, all data were from open-source websites. Landsat images were used for land surface temperature (LST) inversion. Radiometric calibration and atmospheric correction were pre-processed using “Radiometric Calibration” and “FLAASG Atmospheric Correction” of ENVI 5.6 software, respectively. As the study area involves 2 Landsat images, image stitching was performed using “Seamless Mosaic” based on preprocessing. The road data were obtained from OpenStreetMap (www.openstreetmap.org, accessed on 10 March 2024). The land cover data were collected from the global 30 m land cover dataset. According to the Chinese Academy of Sciences land use/cover classification system and the actual situation in Shanghai, the data were reclassified into six types: farmland, forest, grassland, water bodies (including seawater, freshwater, and wetland), construction land, and bare land, using ArcGIS Pro 3.0. Detailed information is listed in

Table 1. Detailed information on the data used.

Data	Time	Detailed Information
Landsat 5	2 August 2003 02:01:44 Coordinated Universal Time	Path: 118, Row: 38 Path: 118, Row: 39. Resolution = 30 m.
Landsat 8	22 August 2022 02:25:32 Coordinated Universal Time	Path: 118, Row: 38 Path: 118, Row: 39. Resolution = 30 m.
Land cover type	2003 and 2022	Six land cover types: construction land, farmland, forests, grasslands, water bodies, and bare land.

.

The UGS within the study area predominantly consists of parks, residential green spaces, and roadside green spaces, which collectively form a regional network. Beyond the outer ring, there are primarily medium-sized patches interspersed with regular farmland. The average area of UGS in Shanghai in 2022 is about 500 m². Therefore, the size of the research unit should not be too small to ensure the integrity of the UGS configuration. Based on the existing studies on the range of UGS cooling, the maximum distance was estimated to be 300 m [31]. Considering the statistical analysis, a square research unit size of 800 m was utilized in this study. Following the exclusion of units near the sea and lake, the study area was divided into 8768 units.

2.3. Land Surface Temperature Retrieval

Given that the United States Geological Survey no longer provides the Atmospheric Correction Parameter Calculator, it is not feasible to conduct LST inversion using the heat radiation equation. Thus, the LST required for this study transitions to using the Landsat Collection-2 Surface Temperature products (https://www.usgs.gov/landsat-missions/landsat-collection-2-level-2-science-products, accessed on 10 March 2024). The Scaling Factor provided on the web page was used to convert the pixel values in the thermal infrared band for Collection 2 Level 2 to the LST values. The formula is as follows:

$$LST=0.00341802\ast {Value}_{pixel}+149.0,$$

(1)

where LST represents the surface temperature (°C); ${Value}_{grid}$ represents the value of each pixel.

2.4. Indicators of UGS Configuration

The landscape pattern index is a quantitative description of the spatial configuration and composition of the landscape [32]. It has been extensively utilized to evaluate the characteristics of land cover change. Compared to composition, landscape configuration refers to the expansion pattern, spatial characteristics, arrangement, and orientation of landscape elements [33,34]. It is more dramatically stressed during the urbanization process. Therefore, this study considered forests, grasslands, and farmlands as UGS, focusing on the impact of landscape configuration changes of UGS on UHI. Based on the previous studies and considering the practical applicability in guiding urban planning, the Large Patch Index (LPI), number of patches (NP), and Landscape Shape Index (LSI) were selected to characterize the UGS configuration (

Table 2. Detailed information on UGS configuration indicators.

UGS Configuration Indicator	Formula	Explanation
Large Patch Index (LPI)	$\mathrm{L}\mathrm{P}\mathrm{I}=\frac{max\left(a_{ij}\right)}{A}\ast 100$	$max\left(a_{ij}\right)$ is the area of the patch in square meters and A is the total landscape area in square meters. 0 < LPI ≦ 100. Units = Percent.
Number of Patches (NP)	$\mathrm{N}\mathrm{P}=n_i$	$n_i$= number of patches in the landscape of patch type (class) i. NP ≥ 1, without limit. Units = None.
Landscape Shape Index (LSI)	$\mathrm{L}\mathrm{S}\mathrm{I}=\frac{0.25E^{\mathrm{*}}}{\sqrt{A}}$	E* = total length (m) of edge in the landscape; includes the entire landscape boundary and some or all background edge segments. A = total landscape area (m2). LSI ≥ 1, without limit. Units = None.

). LPI represents the single UGS patch with the largest area. NP represents the total number of UGS patches. LSI represents the complexity of the UGS patch’s shape within the region. While other indicators like Patch Density and Fractal Dimension Index may also indicate the level of configuration, they do not offer easily interpretable spatial information for decision-makers [35].

LPI signifies the percentage of area covered by a specific UGS type. It was selected because changes in land cover type within patches may be influenced by the abundance of patches with the same type in the surrounding landscape. UGS patches with a higher LPI may pose greater difficulty in changing to another type [36]. NP provides a visual representation of the quantity of UGS patches within the research unit. By examining changes in LPI values within the research units, the increasing or decreasing trends of NP indirectly illustrate the spatial characteristics influenced by urban sprawl [37]. For instance, a decrease in LPI coupled with an increase in NP within the same unit indicates heightened UGS fragmentation and the destruction of important patches. LSI quantifies the proportion of patch boundaries exposed to edge effects [38]. Higher LSI values indicate that a larger proportion of the patch boundary is susceptible to the effects, particularly in UGS patches with smaller areas and more irregular perimeters [39]. In this study, it suggests that UGS are more influenced by heat from adjacent construction land. Analyzing the results of LSI calculations in conjunction with NP can provide a dependable basis for the fragmentation analysis of UGS.

The indicators were computed using FRAGSTATS 4.2 software, which requires specifying the object of analysis. Therefore, patch data for UGS in the study area was extracted, and the coordinate system was converted to WGS_1984_UTM_Zone_50N. The moving window size was set to 800 m.

2.5. Bivariate Moran Analysis

Moran’s I, as the most common indicator for spatial autocorrelation analysis, considering the stochastic phenomenon of the spatial distribution in two or more dimensions, characterizes the clustering and disaggregation of the research units [40]. It is more useful in characterizing spatial patterns than standard landscape indices, especially in Landsat data [41]. Bivariate Moran can reveal the spatial correlation between the independent variable and the dependent variable of neighboring research units. The results visualize how changes in the independent variable affect the neighboring units in terms of spatial distribution [42].

Thus, this study investigates the characteristics of spatio−temporal changes between UGS configuration and LST through global and local Moran statistics. Initially, to characterize the interannual variation of UGS configuration indicators, the difference between 2003 and 2022 was calculated using raster precision:

$$∆I=\overline{I_{2003}}-\overline{I_{2022}},$$

(2)

where $∆I$ is the value of the interannual change in indicator I (in this study I represents LPI, NP, LSI, and LST) from 2003 to 2022, $\overline{I_{2003}}$ is the mean value of the indicator for each research unit in 2003, and $\overline{I_{2022}}$ is the mean value of the indicator for each research unit in 2022.

Secondly, this study utilized global Moran’s I to assess the global spatial correlation between UGS configuration and LST in the study area. When the value is positive, it represents that the two variables are positively correlated with each other as well as they are clustered. The larger the positive value, the greater the positive correlation. When the value is negative, it represents that the two variables are negatively correlated with each other and they are discrete. If the value tends to 0, it shows that the two variables are not correlated with each other [43]. Bivariate local Moran’s I (BLMI) was computed to unveil the spatial correlation between UGS configuration and adjacent LST at the scale of the research unit. As an extension of the univariate local Moran’s I, BLMI assumes that closely related variables tend to exhibit similar spatial patterns [44]. The formula is as follows:

$$I_{kl}=Z_k^i{\sum }_{j=1}^n{W}_{ij}{Z}_l^j$$

(3)

where $Z_k^i=\frac{\left[X_k^i-\overline{X_k}\right]}{{\sigma }_k}$, $Z_k^i=\frac{\left[X_l^j-\overline{X_l}\right]}{{\sigma }_l}$.$X_k^i$ is the mean value of ∆LST for research unit i, $X_l^j$ is the mean value of ∆LPI, ∆NP, and ∆LSI for the neighboring research unit j. $\overline{X_k}$ and $\overline{X_l}$ represent the mean values of the variations of ∆LST and UGS configuration indicators, respectively. ${\sigma }_k$ and ${\sigma }_l$ represent the variance in X for variables ∆LST and UGS configuration indicators. $W_{ij}$ is the spatial weight, which is calculated based on the neighbor matrix.

The results indicate the relationship between a particular research unit and its neighboring units: (1) a High–High (H-H) cluster indicates that the dependent variable of the analysis unit i and the neighboring unit j (∆LST) are both high; (2) a Low–Low (L-L) cluster indicates the unit i and the neighboring unit ∆LST are both low; and (3) the units of Low–High (L-H) and High–Low (L-H) indicate that the independent variable of unit i is low/high, while the dependent variable of the neighboring unit ∆LST is high/low [44]. This study focused on the H-H and L-L clusters as it unveiled spatial aggregation of research units with high UGS configuration indicators and ∆LSI. The above steps were performed by GeoDa 1.8 (http://geodacenter.github.io/, accessed on 10 March 2024) software.

2.6. Geographical Convergent Cross Mapping

A full and accurate understanding of the causal relationship between independent and dependent variables is an important element in the effective use of natural resources and the achievement of sustainable development. Due to spatial heterogeneity and spatial spillover, this means that the potential outcomes of a spatial unit affect its spatial neighbors [45]. Traditional regression models are often used to analyze the relationships between variables within the same research unit, ignoring the impact on the surrounding environment. It only detects the magnitude and direction of the effects of factors on the response variable, but correlation alone cannot be used to infer causality between variables [46]. Thus, the causal model was introduced in this study to try to reveal the causal relationship between UGS configuration and LST.

Current causal inference approaches still have limitations in handling spatial cross-sectional data. Therefore, a novel spatial causal inference model called Geographical Convergent Cross Mapping (GCCM) was proposed [47]. This model effectively addresses the limitations of existing spatio−temporal causal models by employing dynamic systems theory and general embedding theory to identify causal relationships between spatial cross-sectional variables. In this study, ∆LST was considered as the dependent variable, and the spatial causal associations between ∆LST and ∆LPI, ∆NP, ∆LSI were investigated. The modeling process was implemented using the “stats” and “pcalg” packages in R Studio 3.1 [48]. Moreover, to verify the results of causality, this study also used the Pearson correlation coefficient to quantify the strength of the association between variables.

3. Results

3.1. Spatio−Temporal Changes in UGS Configuration and LST

Figure 2 shows the spatial distribution of LST in the study area in 2003 and 2022, as well as the interannual variation characteristics of LST. The high-temperature area of Shanghai expanded from 721 km² in 2003 to 3059 km² in 2022, using a threshold of more than 35 °C to define high-temperature space [49]. The UHI center shifted from its original location in the northern area to the southwest. New areas experiencing high temperatures have emerged in towns outside the outer ring. Surprisingly, there was a decreasing trend in LST observed in the central city within the outer ring road.

For the UGS configuration indicators, LPI exhibited a ring around the central city in 2003, which was due to the presence of continuous large areas of farmland within the suburbs. However, with urban sprawl, there were significantly fewer areas with larger UGSs in 2022, and they showed a decreasing radial pattern from the central city to the suburbs. Regions presenting a reduction in the largest UGS seemed to be related to the road’s form. Meanwhile, regions with significant inter-annual changes in the number of UGS were distributed near the outer ring road and across the suburbs. Compared to 2003, the maximum number of UGS in the study area in 2022 increased from 27 to 77. There was a notable increase in the number of areas with more UGS patches. Changes in the shape of UGS showed characteristics of regional differences, with a decrease in the central city and an increase in the suburbs. There was a general tendency for UGS morphology to be simpler in the research units located in the central city, particularly along the Huangpu River. The regions where the shape of UGS became complex were observed near the suburban ring road in the southwest, directly influenced by urban sprawl. Additionally, spatial overlap was observed in these regions where the number and shape of UGS were simultaneously increasing, indicating a trend of fragmentation and complexity in UGS.

3.2. Spatial Correlation between UGS Configuration and LST

Table 3. Bivariate global Moran’s I among ΔLPI, ΔNP, ΔLSI, and ΔLST.

	Independent Variables			Dependent Variable
	∆LPI	∆NP	∆LSI	∆LST
∆LPI	0.856	−0.580	−0.316	−0.270
∆NP	−0.580	0.794	0.581	0.067
∆LSI	−0.316	0.582	0.729	0.015
∆LST	−0.272	0.070	0.017	0.773

Note: p < 0.05.

shows the global Moran statistics between ∆LPI, ∆NP, ∆LSI, and ∆LST (p < 0.01). The global Moran’s I of ∆LPI was 0.856, which was 7.8% and 17.4% higher than ∆LSI and ∆NP, respectively. Thus, it can be argued that there was a stronger spatial correlation between the largest area of UGS and LST, as opposed to its shape and number. This result indicated that urban sprawl initially impacts the largest area of UGS in the region.

The global Moran’s I for both ∆NP and ∆LSI were negative, implying that the research units with changes in the number and shape of UGS showed a discrete distribution. It was evident that most of the bivariate global Moran’s I was positive, except for the difference between ΔLST and ΔLPI. The negative correlation was strongest between ∆LST and ∆LPI (Moran’s I = −0.270). It indicated a spatial clustering of research units with changes in the largest UGS. Both ΔNP and ΔLSI showed a significant positive correlation with ΔLST, but the correlation coefficients were smaller than those of ΔLPI. This indicated that a change in the largest UGS had a stronger effect on LST in the surrounding area than shape and number. Combining the regions located with significant changes in NP (Figure 2), spatial correlation between ΔNP and ΔLST indicated that the variation in the number of UGS was related to the specific local development pattern.

The BLMI revealed five clusters of correlations between the UGS configuration indicators and the LST (Figure 3), with the white units being non-significant spatial correlations. There was substantial spatial differentiation between the five clusters, as confirmed by the global Moran’s I. For all UGS configuration indicators, the spatial distribution of the H-H and L-L clusters was largely aligned with the regions where changes in ∆LST produced significant alterations (Figure 2). The H-H clusters were labeled as concentrated in the Shong-jiang-Qingpu-Jinshan districts, with a large number of agglomerations of different sizes. In contrast, the L-L clusters were mainly in the Baoshan-Jing’an-Yangpu districts, showing large continuous gathering regions.

Figure 3a,b depict the difference in the distribution of H-L and L-L clusters based on the BLMI between ∆LPI, ∆NP, and ∆LST. Compared to changes in ∆NP, clusters with high ∆LPI and neighboring low ∆LST occupied a larger area in the Baoshan–Jing’an–Yangpu districts and Pudong Coastal, with the spatial distribution becoming more aggregated closer to the central city. H-H cluster characterized by high ∆NP and neighboring high ∆LST was prevalent in the Songjiang–Jinshan districts, with the total number of research units in this cluster being more than ∆LPI but similar to ∆LSI. There was spatial overlap between the L-L clusters of ∆NP, ∆LSI, and ∆LST, suggesting a significant spatial correlation between the reduction in the number of UGS, morphological regularization, and the decrease in LST.

3.3. Driving Factors of UGS Configuration Influencing LST

Table 4. GCMM results for ∆LST and influencing indicators.

Method	Relationship	∆LPI		∆NP		∆LSI
		$\mathit{r}/\mathit{\rho }$	p	$\mathit{r}/\mathit{\rho }$	p	$\mathit{r}/\mathit{\rho }$	p
Pearson Correlation		−0.352 **	0.000	0.011	0.000	0.063 *	0.000
GCCM	x xmap y	0.14	0.00 **	0.17	0.00 **	0.15	0.00 **
	y xmap x	0.19	0.00 **	0.11	0.00 **	0.12	0.00 **

** correlation is significant at the 0.01 level, * correlation is significant at the 0.05 level; r/$\rho$ r is for Pearson correlation, $\rho$ is for GCCM; p is the significance p-value for Pearson correlation and GCCM; xmap means cross-mapping prediction; x means ∆LPI, ∆NP, and ∆LSI; y means ∆LST.

shows the spatial correlation coefficients between UGS configuration indicators and ∆LST, including the results of the Pearson correlation and GCMM. The correlation between ∆LPI and ∆LST was significantly negative, while the correlation between ∆LSI and ∆LST was significantly positive, albeit with a small correlation coefficient. ∆NP showed a non-significant correlation with ∆LST. The GCCM produced results that are numerically close to the Pearson correlation coefficient. x xmap($\to$) y means UGS configuration indicators are influenced by ∆LST, and y xmap($\to$) x means ∆LST is influenced by UGS configuration indicators. The result of the causal association between ∆LPI and ∆LST revealed that the cross-mapping prediction skill ($\rho$) of ∆LST $\to$ ∆LPI was about 0.19, larger than ∆LPI $\to$ ∆LST (Figure 4a). It suggests that the largest patch size of UGS has a strong influence on LST changes. The $\rho$ of ∆LST $\to$ ∆NP and ∆LST $\to$ ∆LSI had smaller values of 0.11 and 1.12, respectively (Figure 4b,c). This result represented a weak causal association between changes in the number and shape of UGS and LST.

It is worth noting the positional change of the blue and red lines in Figure 4b,c, compared to Figure 4a. The $\rho$ of ∆NP $\to$ ∆LST and ∆LST $\to$ ∆NP were 0.11 and 0.17, while the $\rho$ of ∆LSI $\to$ ∆LST and ∆LST $\to$ ∆LSI were 0.15 and 0.12, respectively. These values were very close. Combined with the Pearson correlation coefficient, this indicated that changes in the shape and number of UGS did not significantly influence the LST of the research unit where it was located and its neighboring units. Meanwhile, the coefficient between ∆NP and ∆LST was non-significant, and the coefficient between ∆LSI and ∆LST was significant at the 0.05 level. This phenomenon revealed that the number and shape of UGS may be driven more by the LST changes in the surrounding environment.

4. Discussion

4.1. Spatio−Temporal Differences between UGB Configuration and LST

With urban sprawl, the UGS configuration has been strongly disturbed by human activities, especially in regions on the periphery of the central city. The areas with high LST also gradually changed from the original spatial characteristics of high concentration in one center to multiple dispersed centers (Figure 2). This trend appears logical due to the expansion of built-up areas and the development of the five major new towns [50], where the population serves as a significant factor contributing to the elevated LST [51]. This phenomenon is not only occurring in China, but also in Europe and Africa, where large cities are now facing thermal exposure caused by heatwaves and a steadily growing urban population [52,53].

The interannual change in LPI from 2003 to 2022 was consistent with the irregular sprawl pattern of construction land, which is likely to divide the original large individual UGS patches [54], especially near large water bodies, parks, or roads. This result is consistent with the findings on Mosul city in the Iraq region by Salwan Ali Abed et al. [55]. Simultaneously, it was also related to the land use type; for example, residential areas will be generally close to being situated near an attractive environment and convenient transport [56]. This behavior will inevitably result in a decrease in UGS area and an increase in LST in regions that previously had low LST. From the central city to the suburbs, it similarly appeared that the area where the number of UGS changed was becoming larger and more dispersed. These areas have spatial overlap with units where the largest area of UGS was decreased, suggesting that UGS patches within these units are isolated and fragmented [57]. It is interesting to note that previous studies have found higher LST in areas close to barren land and lower LST in areas close to vegetation and water bodies, which is known as the ‘proximity effect’ [58]. This study further found that the probability of this phenomenon is higher in the suburbs than in the central city. Thus, UGS could be susceptible to external disturbances, such as traffic and crowding [59], which is explained by the fact that the regions with significant changes in LSI are essentially contained within the regions with the largest UGS changes in Shanghai’s suburbs.

This study indicated that in regions experiencing significant changes in the number of UGS, there are often large-scale alterations in the largest UGS. It implies a notable fragmentation of UGS configuration aligned with the process of suburbanization [60]. Additionally, causes of suburbanization may include some policies, such as the retail land use regulation in the Ile-de-France region [61], and policies for limiting land consumption in Germany [62]. According to the theory of the City Life Cycle, suburbanization was regarded as a sequential stage of metropolitan development [63]. Thus, it requires the effective use of small and medium-sized UGS to reduce LST, such as the polygonal UGS with mixed vegetation and close-to-water bodies [64], while controlling the expansion of construction land [65].

4.2. Key Driving Factor of UGS Configuration Influencing LST

By analyzing the spatial correlation between the change values of the UGS configuration indicators and LST, the result indicated that when the largest area of UGS in the region changed, it was more likely that the edges of UGS were encroached upon by construction land rather than being sliced directly into smaller patches [66]. This phenomenon is found frequently in suburban towns, such as the 18 suburban sprawl renovation projects in the United States between 2008 and 2018 [67]. This may also be related to the conversion of pre-existing intact farmland into construction land [20]. Meanwhile, the bivariate global Moran’s I implied that the research units in which the largest UGS and LST had changed were spatially separated from each other, while those in which the number of UGS and LST had changed showed aggregation. The separated distribution may be related to the expansion of villages into large areas of farmland, which is common in the sprawl of cities with flat topography [68]. Studies have proved that urban regeneration projects and continuous green spaces on the periphery of central cities prevent new sprawl [69]. The aggregation distribution can be attributed to urban regeneration in the central city, such as the construction of pocket parks and street trees, which can effectively reduce regional LST [70].

BLMI revealed that areas of an increased largest UGS area with an increase in surrounding LST were distributed between the outer ring road and the suburban ring road. This phenomenon, which could be attributed to urban sprawl and the protection policy of basic farmland in China [71], has resulted in suburbanization generating more heat than is reduced by the addition of new UGS. Areas of declining UGS numbers and regularized shape were within the outer ring road, which already has well-developed infrastructures and is less affected during urban sprawl [25,72]. In fact, with the planning of green belts, wedge-shaped UGSs, and green corridors, the focus of UGS construction in Shanghai has shifted from the central city to the suburbs [73]. It is worth noting that the L-L clustering of the number, and shape of UGS and ∆LST has spatial overlap mostly located in the central city. It is related to the regular shape of the UGS as a result of urban regeneration, and the increased number of trees within the same UGS, which could create a secondary effect of the cold island [74], or the increase in low-albedo land in the surrounding area [75].

In verifying the key factors of UGS configurations affecting the neighboring LST, previous studies have consistently confirmed a negative relationship between change in the largest UGS and ∆LST [24,76]. This is consistent with the causal association between ∆LPI xmap($\to$) ∆LST, as examined by GCCM in this study. When a region contains larger contiguous UGS patches, it is likely to decrease LST as these patches provide more shaded areas and vegetation cover, reducing the absorption of solar radiation [77,78]. For ∆NP/∆LSI $\to$ ∆LST, both $\rho$ values were smaller than ∆LPI $\to$ ∆LST. For the same number of patches, differences in vegetation types within the UGS can significantly impact the cooling effect on the surrounding environment. For instance, albedo values at the surface typically decline from grasslands with standing dead leaves, which have high reflectivity, to broadleaf forests, and finally to dark coniferous forests [79]. In other words, if a UGS patch is covered by grasslands, its causal association with ∆LST is no longer significant. Thus, reducing LST requires attention to the largest UGS in the region and the vegetation composition within it.

It is worth noting that, as explained by Suighara et al. [80], the smaller $\rho$ of ∆LST $\to$ ∆NP/∆LSI can either indicate a weaker causal influence of ∆LST on ∆NP/∆LSI, or the spatial variation in ∆LST can have a small reflection on the spatial variation of ∆NP/∆LSI (reflection, not causation). This study indicated that regions with changes in LST are often accompanied by corresponding changes in the number and shape of the UGS. The key factor of UGS configuration that has a significant effect on the neighboring LST is the size of the largest UGS. This result is consistent with a portion of existing research. For example, He, B et al. found that when the total amount of green space is fixed, the cooling effect is stronger in a few large cores than in a large number of small islands [81]. However, some studies have found the shape of the landscape/UGS was the key factor affecting the regional temperature [82,83]. These differences probably relate to the inversion methods for the thermal environment as well as to the study object, which needs to be explored in a wider range of cases.

4.3. Implications for Planning Practice and Policy

Green space is listed as one of the indicators for achieving the Sustainable Development Goals by 2030. The results of this study emphasized the UGS configuration as an effective tool to promote urban cooling. Compared to the number and shape of UGS, the area of the largest UGS patch, which exhibits strong causality and significant negative correlation with ∆LST, suggests that increasing the area of the largest green space in the region is effective in mitigating localized high temperatures. This pattern remains broadly consistent across cities with different climatic geographies [84]. While there are already many examples of using UGS to combat heat waves in the planning of developed Chinese cities, currently there is no policy to match this. For example, policies of UGS planning do not feature prominently on the ‘China Sustainable Urban Cooling Project’ (in Guangzhou). Based on this study’s findings, the challenge is not only to protect or rehabilitate the largest UGS in the region but also to avoid fragmentation of UGS. This requires policymakers to adopt different approaches to UGS planning and management according to regional characteristics. For instance, in suburbs, construction land development should be traded off against the protection of original ecological space; in central cities, reasonably distributed and sufficiently sized UGS should be preserved [85]. It also emphasizes that the construction of UGS should be treated as a long-term project, with timely adjustment of UGS restriction policies at different stages of urban development.

For planning practice, planners and designers should recognize that cities need to create wider cold islands to mitigate the UHI effect [86]. Regions located at the edge of urban sprawl show a trend of decreasing UGS number and regular shape, resulting in lower LST in the neighboring areas. Regions located in the suburbs show a trend of a decreasing largest area and increasing number, resulting in higher LST in the neighboring areas. These phenomena emphasize the importance of tailored strategies for mitigating UHI effects in different zones. In addition, a focus on policy and planning practice is a key step in preserving ecological spatial quality and connectivity in sparsely populated suburbs, compared to the central city.

4.4. Limitations and Future Research

The limitations of this study were clarified, mainly the selection of UGS configuration indicators and the limited number of independent variables. On the one hand, the construction of the index system does not contain all the landscape pattern indicators. Although the landscape pattern index, which has a significant effect on LST in general, was chosen for this study, there may be more significant factors for this study area. On the other hand, the object of causal interpretation was restricted to the UGS, and indicators describing the urban sprawl, such as land use intensity, population, and nighttime lighting, were not analyzed [87]. This study did not explicitly establish a strong spatial causal relationship between a specific type of landscape pattern or land use type with LST change.

Based on the comparative analyses between this study and the existing studies, suggestions for future research were proposed. Firstly, this study argues that more attention needs to be paid to the factors that lead to LST increase, as the results of this study are controversial with existing studies. For example, Ward et al. (2016) pointed out that urban form indicators are the most critical factors [88]. He et al. (2019) explained that the effect of urban form on UHI depends on the size of the study area [89]. Therefore, there is a need to include indicators other than UGS patterns in future studies, such as exploring the effectiveness of different UGS management strategies on LST or investigating socio-economic factors affecting UGS distribution. Secondly, with the evolution of machine learning and spatial regression models, better tools are available for analyzing the impact of UGS and socio-economic factors on LST. Instead of being limited to simple regression and linear relationships, studies can try to clarify what types and forms of UGS can effectively reduce LST from the perspective of spatial and temporal heterogeneity. It can provide more accurate recommendations for UGS planning and policy.

5. Conclusions

In this study, a novel approach was proposed to elucidate the spatial heterogeneity and its causality between LST changes and the UGS configuration. The results indicated that the high-temperature space shifted from the central city to the suburbs, expanding from 721 km² in 2003 to 3059 km² in 2022. A decreasing trend of LST was observed in the central city. There was a stronger spatial correlation between ∆LPI and ∆LST, compared to ∆LSI and ∆NP, which indicated that urban sprawl primarily impacted the UGS patches with larger areas. A significant causal association was also shown between these two variables. Conclusively, this study provided a new insight into the spatial relationship between the change in UGS configuration and LST interannually. The results of the study further confirmed that the rapidly urbanizing region, Shanghai, showed significant spatial heterogeneity in the changes in UGS configuration and LST. The LPI of UGS was more strongly causally related to the changes in LST, compared to NP and LSI. In addition, UGS configurations in adjacent research units need to be considered in UGS plans and policies for UHI mitigation. The series of results derived from this study illustrates the effectiveness of the bivariate local Moran statistics in identifying the spatial externality of UGS configuration changes to the surrounding LST. Significant spatial correlations between variables do not represent strong causality.