Multi-scale vegetation cooling assessments for urban climate adaptation: a MAUP-informed NDVI–LST analysis in Shanghai

Over the past few decades, global urbanization has accelerated markedly. The proportion of the global population living in cities has risen from approximately one-third in the mid-20th century to more than half in 2008 (Grimm et al., 2008), and is projected to reach two-thirds by 2050. This rapid urbanization, accompanied by extensive land-use changes and landscape pattern restructuring, is widely recognized as a key driver of regional and global warming (Rizwan et al., 2008; Cao et al., 2016; He and Zhu, 2018).

During urban expansion, natural and semi-natural surfaces are continuously replaced by impervious surfaces and high-density built-up areas, significantly altering surface energy exchange and thermal processes (Li et al., 2020). Consequently, urban areas typically exhibit higher land surface temperatures (LSTs) compared to surrounding rural regions (Foley et al., 2005; Stewart and Oke, 2003), a phenomenon widely referred to as the Urban Heat Island (UHI) effect (Oke, 1982; Ward et al., 2016).

Numerous studies have shown that UHI is not only a significant component of the urban climate system but also has profound impacts on public health, energy consumption and the stability of urban ecosystems. The UHI effect is closely linked to various common diseases and health risks (Patz et al., 2005) and may significantly increase mortality under extreme heat conditions (Mora et al., 2017; Di Blasi et al., 2023). In the context of climate change, the combined effect of UHI and global warming is expected to further intensify the warming trend and heat-related risks in urban areas (Hoag, 2015; Hao et al., 2025). Therefore, as the frequency of extreme heat events continues to rise, adaptation and mitigation strategies targeting the UHI effect have become a critical issue in contemporary urban climate research and policy-making (Ouria et al., 2025).

In this context, urban green spaces are widely recognized in environmental management and planning as important natural infrastructure for mitigating the UHI effect. Many studies use the Normalized Difference Vegetation Index (NDVI) as a key indicator of vegetation cover. They have confirmed that NDVI generally shows a negative correlation with LST, and this relationship is often used to evaluate the cooling potential of urban green areas (Sharmin et al., 2024). However, most existing research is conducted based on landscape pattern indices and under conditions of a single time period or season (Liu et al., 2024), with limited focus on how the choice of scale itself affects the statistical results of the NDVI–LST relationship. As urban spatial structures and land-use patterns continue to evolve, the method of aggregating analysis units and the selected scale can significantly change the observed relationship between vegetation and temperature, including its relationship and the identification of key thresholds. These differences can lead to incorrect assessments of the cooling potential of urban green spaces in environmental management practices.

Especially under rapid urbanization, where different land use types are continuously changing, the NDVI-LST relationship is often embedded in complex spatial mixing and scale dependence. If the Modifiable Areal Unit Problem (MAUP) is ignored, conclusions drawn from a single scale may be difficult to compare effectively across different spatial levels. This can reduce their functions for urban heat environment assessment and planning decisions.

Thus, this study takes Shanghai as a case city. Under the background of land use change, it constructs a multi-scale regular grid analysis framework. From the MAUP perspective, it systematically examines how different spatial aggregation scales affect the statistical form and threshold identification of the NDVI–LST relationship. It further discusses the impact of scale dependence on the assessment of urban cooling potential and environmental management decisions.

Existing studies show that vegetation plays a key role in regulating the urban thermal environment by influencing surface energy balance and hydrothermal processes (Morabito et al., 2021). Many researchers have demonstrated that increasing the area of urban green spaces can mitigate the UHI effect to varying degrees (Xiao et al., 2018; Zhou et al., 2022a, 2022b). From the perspective of mechanism, urban green spaces (UGS) primarily lower LST through evapotranspiration cooling and shading effects, creating relatively cool zones within cities, thereby alleviating the UHI effect (Wang et al., 2025; Gomez-Martinez et al., 2021; Zhou et al., 2025a, 2025b). These studies, based on various spatial scales and urban contexts, collectively verify the general regulatory capacity of vegetation coverage on the urban thermal environment.

However, while the negative correlation between vegetation and LST has been widely acknowledged, existing research has predominantly relied on average analysis at the overall or macro scale, often overlooking the spatial heterogeneity in the cooling effects of different types of green spaces. Further studies reveal significant differences in cooling range and intensity among various green space types. For instance, forest-type green spaces typically exhibit a much larger cooling influence area than scattered lawns, while the composite spatial structure of water bodies and green spaces may demonstrate stronger synergistic cooling effects (Feyisa et al., 2014; Li et al., 2013; Zhou et al., 2023). Furthermore, the cooling effect of green spaces shows distinct nonlinear characteristics across different seasons, geographical locations and climatic contexts, and may involve critical thresholds (Yu et al., 2025a, 2025b). Meanwhile, multi-city studies indicate that vegetation coverage and spatial configuration jointly shape the urban thermal environment. Their effects are modulated by climatic conditions and urban morphological backgrounds (Wang et al., 2024). These findings suggest that the thermal regulation function of urban green spaces cannot be accurately described by uniform parameters. Instead, analysis needs to integrate specific spatial scales and contextual conditions.

However, studies exploring the relationship between urban land use structure, vegetation cover and the thermal environment from different perspectives still face significant limitations in terms of temporal and spatial scale handling. On one hand, existing research predominantly focuses on single time points or short time series, which limits understanding of how vegetation-temperature relationships are conditioned by evolving urban development backgrounds (Wu et al., 2025; Liu et al., 2020). In the context of rapid urbanization, such single temporal approaches may obscure the phased differences and cumulative effects in the NDVI–LST relationship across different development stages. This could weaken the overall understanding of urban thermal environment regulation mechanisms under climate change scenarios.

On the other hand, in terms of spatial dimensions, numerous studies typically conduct analysis using a single spatial scale and have demonstrated the statistical correlation between vegetation indicators and LST across different urban contexts (Huang et al., 2024; Duan et al., 2025). However, many studies indicate that the NDVI–LST relationship exhibits significant sensitivity to spatial scale. For instance, at smaller scales, such as blocks or communities, vegetation cover characteristics often show stronger explanatory power for LST. In contrast, at regional or city-wide scales, this correlation may significantly weaken and findings across different studies are not entirely consistent (Kim et al., 2025; Chen et al., 2023).

This discrepancy is largely related to differences in scale processing methods within research designs, including the selection of spatial resolution for remote sensing imagery (Zhou et al., 2011; Li et al., 2013) and the methods for delineating and aggregating spatial analysis units (Kong et al., 2014; Jang and Jung, 2025). Relevant studies generally indicate that the statistical relationship between vegetation indicators and LST exhibits typical scale-dependent characteristics (Rahimi and Dong, 2025; Karami et al., 2025). In recent years, some studies have systematically examined thermal environment characteristics under different spatial scales and zoning methods by introducing advanced statistical approaches and machine learning models. This has led to explicit discussions on the MAUP within LST analysis (Deng et al., 2024; Lee et al., 2025). These studies provide important references for understanding the underlying drivers of thermal environment changes. In summary, systematically comparing the statistical performance of the NDVI–LST relationship across different spatial scales can improve the interpretability and comparability of findings.

While the statistical relationship between NDVI and LST has been widely explored, studies systematically examining how multi-period urban development backgrounds interact with spatial aggregation scale to affect the interpretability of the NDVI–LST relationship remain relatively scarce. By keeping the spatial resolution of remote sensing images unchanged, and investigating how changes in analysis unit scale affect the statistical representation of the NDVI–LST relationship from the perspective of the MAUP, we can gain a clearer understanding of the methodological impact of scale setting in UHI research. This study aims to reveal how scale selection influences the interpretation of the NDVI–LST relationship across different land cover types, and how this scale-dependent pattern may lead to misinterpretations in urban thermal environment management.

Based on this, this study takes Shanghai as a case study. It establishes a multi-scale regular grid analysis framework, considering the context of multi-phase land use change. A land use transition matrix is introduced to define different urban development scenarios. The structural characteristics of the NDVI–LST relationship across different spatial aggregation scales are systematically examined. The research focuses on how the spatial aggregation process, along with its interaction with land use transition scenarios, affects the statistical results.

This study aims to:

• Describe the spatial distribution characteristics of LST and NDVI under different development stages in Shanghai.

• Compare the statistical differences in the NDVI–LST relationship across various spatial aggregation scales, based on different transition scenarios defined by the land use transition matrix.

• Examine how scale selection (MAUP) and land-use context shape the interpretability of NDVI-LST evidence for urban thermal management under land-use change.

Shanghai (121°28′E, 31°14′N) is located on the eastern coast of China. It is a typical large-scale, high-density city in the country, covering a total area of approximately 6340.5 km², with a permanent population of about 24.85 million (Figure 1). The region has a subtropical monsoon climate, with an annual average temperature of around 16 °C–20 °C and an annual precipitation of about 1200 mm. During the study period, Shanghai experienced significant urban expansion and land use restructuring. Built-up areas continuously expanded, and the spatial pattern of vegetation was continuously reorganized. This ongoing urban development process provided a dynamic background for analyzing the statistical performance of the relationship between vegetation cover and LST across different spatial scales (Huang et al., 2017). Therefore, this study selected Shanghai as the case area. Against the backdrop of multi-period land use changes, it explored the influence of spatial scale and aggregation methods on the interpretation of the NDVI–LST relationship from the perspective of the MAUP.

This study used multi-temporal Landsat remote sensing images, land use/land cover (LULC) data and auxiliary topographic data. The remote sensing images were sourced from the USGS Earth Explorer platform, using Landsat 7 and Landsat 8 data with a spatial resolution of 30 m. To minimize the impacts of cloud cover and seasonal variation, images from four summer periods (July–August) in the years 2000, 2010, 2017 and 2024 were selected, all with cloud cover controlled below 10%, to reflect thermal environment characteristics during high-temperature urban periods.

LULC maps were created using Landsat imagery from the same periods as the LST data. A combination of automated classification and manual correction was applied to ensure the spatial and temporal consistency of the multi-temporal data. Based on the research needs, land cover was grouped into four categories: vegetation land, built-up land, water body (refers to inland waters within the Shanghai administrative area, sea is excluded from analysis), and other land uses(including bare land, marshes and tidal flats). Land use transition matrices for multiple periods were subsequently constructed to describe the structural changes in land use during the study period.

Through the above processing, a multi-temporal remote sensing data set with spatial and temporal alignment was established. This data set provides the data foundation for subsequent analysis of the statistical characteristics of the NDVI–LST relationship under different land use transitions.

LST was used to characterize the urban thermal environment at different development stages of the study area. LST was retrieved from the thermal infrared bands of multi-temporal Landsat images using a single-channel algorithm. This method has been widely applied in long-term LST studies based on Landsat data (Artis and Carnahan, 1982; Qin et al., 2001; Sobrino et al., 2004). It has also been proven to provide stable and comparable results across multi-source remote sensing data sets.

The NDVI was used to characterize vegetation cover conditions and spatial distribution patterns in the study area. It served as an important indicator reflecting urban green space conditions. NDVI was calculated from remote sensing images corresponding to the same periods as the LST data (Rouse et al., 1974). This ensured temporal consistency between the vegetation information and the thermal environmental conditions. Through a unified data pre-processing workflow, the NDVI data maintained consistent spatial resolution and projection systems across different years. This consistency supported cross-temporal and cross-scale statistical comparisons. In this study, NDVI was used as a representative variable of vegetation status to analyze its statistical relationship with LST.

Land use/land cover (LULC) data were used to identify different land use types and their changes during the study period. This provided the basis for contextual classification when analyzing the NDVI–LST relationship under the background of urban development. A land use transition matrix was constructed from the multi-temporal LULC data. This matrix identified major transition pathways between land use types, thereby defining different land use transition scenarios. To ensure consistency in scale analysis, the LULC data were uniformly projected and spatially aligned with the NDVI and LST data. In subsequent multi-scale analyses, they were used only for spatial grouping and result comparison.

Statistical and spatial analysis biases induced by the MAUP typically manifest in two aspects: the scale effect and the zoning effect (Openshaw, 1984; Wong, 2004). The scale effect refers to the statistical variability and heterogeneity of geographic phenomena observed at different spatial scales. It is closely related to the size of the spatial units used during aggregation processes (De Andrade et al., 2021; Wong, 2004). In urban thermal environment studies, changes in the scale of analysis units can significantly affect LST and its statistical relationships with environmental variables.

Specifically, the study area was divided into three regular grid systems: 30 × 30, 60 × 60 and 90 × 90 grids (Figure 2). This created a multi-scale analysis framework from coarse to fine. Among them, the 30 × 30 grid corresponded to the largest individual analysis unit, with an area of approximately 7.05 km², reflecting the spatial aggregation characteristics at the sub-district or block-group scale. The 60 × 60 grid had a unit area of about 1.76 km², corresponding to a community-block mixed scale. The 90 × 90 grid formed the finest-scale analysis unit, with an area of approximately 0.78 km², approaching the neighborhood scale. These grid systems were generated as independent regular fishnet units, because the study focused on MAUP-related sensitivity in areal aggregation units. NDVI, LST and LULC were then overlaid with each grid system, and zonal statistics were calculated consistently for each grid cell at each scale. Regular fishnet grids were used to isolate MAUP-related scale effects under a controlled aggregation framework, whereas irregular zoning effects were beyond the scope of the present study.

To analyze the performance of MAUP-related scale effects under different land use contexts and how these effects impact NDVI-based cooling indicators, the research was conducted in two steps. First, we examined the variation in the response curves between NDVI and LST across different spatial scales for various land use types (LULC), using land use as a contextual control variable. Second, for each grid scale (30 × 30, 60 × 60 and 90 × 90), we trained a Random Forest regressor with NDVI as the predictor and LST as the response. To ensure reproducibility, samples were split into training and testing sets (70 / 30) using a fixed random seed (random_state = 42). The model was implemented in scikit-learn with n_estimators = 200 and max_depth = 6 (other parameters kept as default). Model performance was evaluated on the held-out test set using R², RMSE and MAE. Robustness was further assessed using five-fold cross-validation (KFold with shuffle = True, random_state = 42), and the mean ± SD of cross-validated R² across folds is reported (Table 1).

The RF + SHAP framework was used to capture potentially nonlinear NDVI-LST responses across scales and to identify interpretable turning points within a consistent analytical framework. Compared with using separate simple nonlinear fits, this framework allowed the response pattern and its contribution structure to be examined in a unified way across different aggregation scales.

Figure 3 shows the spatial distribution of LST in Shanghai across four representative years. During the study period, LST exhibited an overall intensification, with a clear spatial expansion of high-temperature areas. In 2000, elevated LST values were mainly concentrated within the urban core area, while surrounding areas generally maintained lower temperature levels. Starting from 2010, high-temperature areas gradually expanded toward suburban areas. By 2017 and 2024, the increase in LST was no longer confined to the city center but became increasingly evident in newly developed peri-urban regions. In contrast, water bodies and coastal areas consistently exhibited lower LST values throughout all periods, highlighting their role as stable cooling backgrounds. Furthermore, the spatial heterogeneity of LST increased over time, with later periods characterized by more fragmented and heterogeneous thermal patterns. This indicates the growing complexity of the urban thermal environment alongside urban development.

The multi-temporal LULC results show that Shanghai’s land use pattern experienced significant changes from 2000 to 2024 (Figure 4). Overall, urban development was characterized by the continuous expansion of built-up areas, accompanied by structural adjustments and spatial reorganization of vegetation land.

In 2000, built-up areas were primarily concentrated in the central urban core and its immediate surroundings, exhibiting a relatively compact spatial distribution. By 2010, built-up areas had expanded noticeably outward, forming continuous urban expansion belts and gradually extending into peripheral regions. After 2017, the expansion pace of built-up areas slowed down, but features of infill development and functional intensification became more pronounced. By 2024, the spatial pattern of built-up areas had stabilized, showing an evolution characteristic of slowed outward expansion and intensified internal densification.

The total area of vegetation land remained relatively stable during the study period, but its spatial distribution changed significantly. On one hand, some original vegetation land was converted to built-up areas, mainly located along the urban expansion frontier. On the other hand, certain degrees of green space restoration and reallocation occurred within the city and along waterfront areas, reflecting dynamic adjustments in green space structure during urban development. Water bodies exhibited minimal overall change, mainly distributed along coastal and river channel areas, demonstrating strong spatial stability.

The land use transition matrix further quantitatively reveals the transformation pathways and their scale characteristics between different LULC types across various periods [Tables 2 (a)–(c)].

During the 2000–2010 stage, the conversion of vegetation land to built-up areas was the most significant. The converted area reached 1389.98 km², making it the primary source of built-up area expansion during this phase. Meanwhile, the original built-up areas maintained high stability. Approximately 2390.89 km² of built-up area remained unchanged in this stage. The conversion scale of water bodies and other land uses was relatively small, with overall limited changes.

During the 2010–2017 stage, the expansion of built-up areas continued, but the transformation pathways became more diverse. The area of vegetation converted to built-up land decreased to 416.72 km², indicating a slowdown in large-scale outward expansion. At the same time, the stable area within existing built-up zones increased significantly (3418.22 km²). This suggests a gradual shift in urban spatial development from outward expansion to internal integration and consolidation.

During the 2017–2024 stage, the overall scale of LULC conversion further decreased, and the land use pattern stabilized. The scale of bidirectional conversion between vegetation and built-up areas remained relatively small. The stability of built-up areas strengthened further(3347.63 km²). The transformation characteristics of this stage indicate that urban development entered a phase primarily focused on stock adjustment and structural optimization.

Overall, the three-stage transition matrices clearly reflect Shanghai’s evolution from rapid expansion to relatively stable development. The vegetation to-built-up land conversion remained the dominant pathway throughout. However, its intensity changed significantly over time. This provided important contextual support for the subsequent analysis of NDVI-LST response characteristics under different LULC transition scenarios.

To reveal the regulatory mechanism of vegetation cover on LST under different land use backgrounds and its scale sensitivity, NDVI-LST response curves were constructed at three grid scales: 30 × 30, 60 × 60 and 90 × 90. These curves were then compared and analyzed across four land use types: Built-up, Vegetation, Water and Other.

In the built-up area context [Figure 5(a)], a consistent nonlinear structure between NDVI and LST was observed across all scales. As NDVI gradually increased from negative values to approximately the 0.10–0.20 range, LST rose significantly, reaching a local peak near this interval. When NDVI further increased, LST declined slowly afterward. The overall pattern showed an inverted U-shaped response of “rising first, then falling.”

Regarding scale effects, the curve at the 30 × 30 coarse scale was the smoothest, with a broader variation range near the peak. This indicates that under larger spatial aggregation, local thermal environment differences were somewhat weakened. As the scale progressively refined to 60 × 60 and 90 × 90, the slope of the curve in the medium-to-low NDVI range increased noticeably, and the peak position became more concentrated. This suggests that at finer scales, NDVI responds more sensitively to local thermal environment variations. Furthermore, within the high NDVI range, the declining trend of LST was clearer at finer scales. This indicated that the cooling effect of vegetation within built-up areas is more easily identified in finer-scale spatial structures.

Within vegetation land [Figure 5(b)], the NDVI–LST relationship exhibited a monotonic increase that gradually leveled off across all three scales. When NDVI rose from low to moderate levels, LST correspondingly increased. However, once NDVI exceeded approximately 0.20–0.30, the rate of LST increase weakened significantly, entering an approximately saturated state.

The shape of the curves was generally consistent across different scales, with only slight variations in the low NDVI range. At the 30 × 30 coarse scale, LST values in the low NDVI region were relatively higher. In contrast, LST values were slightly lower at the finer 60 × 60 and 90 × 90 scales. This suggests that scale refinement helped reduce the high-temperature bias introduced by spatial mixing at coarser scales. Overall, the NDVI-LST response in vegetation land demonstrated strong stability under scale changes.

The NDVI-LST response curve for water bodies [Figure 5(c)] exhibited noticeable fluctuation characteristics. In the low NDVI range, LST increased with rising NDVI, forming a local peak near zero. Subsequently, LST changes were minimal in the high NDVI range, with the overall trend stabilizing.

In scale comparison, the curve showed the largest amplitude of fluctuation at the 30 × 30 coarse scale. At the finer 60 × 60 and 90 × 90 scales, the curves became generally smoother, with extreme fluctuations significantly reduced. These results indicated differences in the thermal response structure related to water bodies across scales.

In the Other land category [Figure 5(d)], only the 30 × 30 coarse scale contained sufficient samples to support analysis. The NDVI-LST curve in this category showed that LST increased with rising NDVI in the low NDVI range, reached a local peak near medium NDVI levels and then gradually declined. Due to the limited sample size at the 60 × 60 and 90 × 90 scales, cross-scale comparative analysis was not conducted for this land type.

In summary, significant differences existed in the morphological characteristics and scale sensitivity of NDVI-LST response curves across different land use types. Built-up areas were most sensitive to scale refinement, with local nonlinear features becoming more prominent at finer scales. Vegetation land exhibited strong stability across scales. Water bodies were significantly affected by spatial mixing effects at coarse scales, while their response relationship became clearer at finer scales. These findings indicated that the role of scale effects in the NDVI–LST relationship has a clear dependency on land use type.

To further quantify how NDVI influences LST across different spatial scales, we used the RF + SHAP framework to examine scale-specific contribution patterns and turning points (Figure 6).

At the coarse scale (30 × 30), the distribution range of NDVI SHAP values was relatively wide, showing a clear bidirectional contribution pattern. Low NDVI samples (blue points) were mainly concentrated on the negative side of the SHAP axis, whereas higher NDVI samples gradually extended toward the positive side. This indicates that NDVI-related contributions to model output varied substantially across samples at this scale. However, SHAP values at this scale showed strong dispersion, with positive and negative contributions highly mixed. This suggests that, under coarse-scale aggregation, heterogeneous land surface types and micro-environmental conditions were mixed within the same grid cells, weakening the separability and interpretability of NDVI-related effects on LST. Consequently, the contribution pattern of NDVI was susceptible to spatial averaging effects.

At the medium scale (60 × 60), the distribution structure of NDVI SHAP values became clearer. As NDVI increased, SHAP values showed a smoother transition across the contribution gradient. The concentration of high NDVI samples in the positive SHAP region was notably enhanced. This suggests that NDVI-related contributions became more structured and stable at this scale. Compared to the coarse scale, the dispersion of SHAP values decreased. The boundary between positive and negative contributions became more distinct.

At the fine scale (90 × 90), the distribution pattern of NDVI SHAP values became the most concentrated and consistent. High NDVI almost entirely corresponded to positive SHAP values, while low NDVI was predominantly concentrated in the negative SHAP area. This indicated that under fine-scale conditions, NDVI-related contribution patterns to LST exhibited high stability and consistency.

Overall, as the spatial scale transitioned from coarse to fine, the regulating effect of NDVI on LST shifted from being discrete and unstable to becoming concentrated and consistent. This indicated that scale refinement helped reduce spatial mixing effects and enhanced the interpretability of the NDVI-LST relationship. However, the SHAP summary plot mainly reflected the overall distribution of contribution direction and intensity. It was insufficient for interpreting the continuous response structure across different NDVI intervals. Therefore, the next section further analyzed the SHAP-NDVI relationship to identify the critical transition ranges of the cooling effect.

Based on the overall contribution distribution revealed by the SHAP summary plot, the SHAP-NDVI relationship was further analyzed to identify critical transition points in the response structure (Figure 7). The results showed that the response structure of NDVI’s cooling contribution exhibited systematic changes across different spatial scales, reflecting a clear scale effect.

At the coarse grid scale (30 × 30), a relatively stable structural turning point was identified near NDVI≈−0.12. The magnitude of change in marginal contribution weakened notably near NDVI≈−0.19, and the curve leveled off. This indicated that under strong spatial aggregation, the NDVI–LST relationship exhibited a more stable and identifiable segmented structure. The larger grid units reduced local fluctuations through statistical averaging effects.

In comparison, at the medium (60 × 60) and fine (90 × 90) scales, the initial structural turning points did not exhibit stable and consistent numerical characteristics. However, near NDVI≈-0.32 and −0.31, segments where the marginal contribution changes slowed down could still be identified. As the scale refined, the local fluctuation amplitude of the response curve increased significantly. This indicated that reduced spatial aggregation strengthened the expression of local variability, causing the contribution structure of NDVI-LST to display more pronounced nonlinear characteristics.

In summary, the structural turning points identified by SHAP differed across grid scales. This showed that the nonlinear response structure of NDVI-LST exhibited a clear scale dependency. These results suggested that the segmented response characteristics identified within a nonlinear modeling framework were also influenced by the method of spatial aggregation. Therefore, scale effects were not only reflected in parameter estimation differences within traditional linear models, but also in the stability of the nonlinear response structure and its threshold identification. In this study, the turning points are better understood as context-specific structural features under changing urban development backgrounds, rather than as temporally invariant values confirmed for each period.

This study finds that the scale effect in the NDVI–LST relationship does not exist in isolation. Instead, it is embedded within the broader context of urban land use evolution. This aligns with existing findings that rapid urbanization enhances and restructures the cooling capacity of vegetation (Wang et al., 2023). As urbanization progresses, the expansion of built-up areas and the restructuring of land cover significantly alter the spatial configuration of vegetation. Consequently, NDVI carries different physical and functional meanings across various land use types (Weng et al., 2004; Zhou et al., 2011; Li et al., 2013). In built-up areas, NDVI often corresponds to scattered or fragmented green patches with weak spatial continuity. In vegetation areas, NDVI is more likely to reflect relatively continuous vegetation cover (Zhang et al., 2026; Khan et al., 2025). Therefore, when analyzing the regulating effect of NDVI on LST, the scale effect is not merely a statistical phenomenon. It is also closely linked to the reorganization of land use structure during urban development.

This context explains why the shape and scale sensitivity of the NDVI-LST response curve varied significantly across different land use types. Land use change, by reshaping the spatial organization of vegetation, provided a physical basis for the emergence of scale effects. As a result, the same NDVI value could correspond to different thermal environment responses under varying spatial and developmental contexts.

This study reveals that the NDVI-LST relationship is highly sensitive to the scale of the analysis unit. As the grid scale changes from coarse to fine, the NDVI-LST relationship does not simply show a change in correlation strength. Instead, a systematic restructuring occurs in its response characteristics and stability. At the coarse scale, the spatial aggregation process smooths local heterogeneity. This results in a relatively continuous and stable nonlinear structure in the relationship between NDVI and LST. At the medium and fine scales, the influence of local land cover composition, boundary effects and mixed pixels is significantly amplified. This causes the response relationship to become more dispersed and variable. These findings are consistent with previous studies emphasizing the close relationship between LST and spatial scale effects (Fu and Weng, 2016; Wang et al., 2016). When the form of spatial units changes, attribute data may be re-aggregated or dis-aggregated. This leads to changes in the attribute data themselves and their correlations (Josselin and Louvet, 2019; Luan et al., 2020).

This result suggests that scale selection does not merely affect the numerical magnitude of statistical outcomes. More fundamentally, it shapes how researchers understand the form of the NDVI–LST relationship. If the scale effects introduced by MAUP are ignored, scale-specific response structures may be misinterpreted as general laws (Openshaw, 1984; Wong, 2004). This undermines the interpretability of research findings in cross-scale comparisons.

This study further demonstrates that the key segmented characteristics within the nonlinear NDVI-LST response exhibit significant scale dependency. This finding aligns with previous research on the nonlinear nature of the NDVI–LST relationship and its sensitivity to scale (Zhou et al., 2022a, 2022b; Zhou et al., 2025a, 2025b; Yu et al., 2020). The cooling onset turning point only manifests as a relatively stable and identifiable structural feature at the coarse scale. At the medium and fine scales, it does not form a consistent numerical interval. This indicates that the identification of the turning point itself is influenced by spatial aggregation conditions.

At the coarse scale, spatial heterogeneity is statistically averaged out. This allows the NDVI-LST response to exhibit a relatively continuous segmented structure. At the fine scale, local variations amplify the volatility of the response curve. This makes it difficult to stably identify a single numerical interval. In contrast, the later segment where the response levels off can still be identified across different scales. However, its position shifts systematically with scale. This further reflects the sensitivity of the nonlinear response structure to scale conditions (Zong et al., 2025; Yu et al., 2025a, 2025b; Luo et al., 2026). This finding helps explain the variability in NDVI cooling thresholds reported in existing studies and suggests that such thresholds are better understood as scale-dependent and context-specific rather than universally fixed. More importantly, the contribution of this study lies in showing that threshold-like NDVI-LST transitions are shaped by spatial aggregation rather than year-by-year threshold stability. This interpretation is consistent with recent evidence showing that vegetation cooling efficacy varies with climatic background, urban morphology and vegetation traits (Rahman et al., 2020; Li et al., 2024). A brief year-specific supplementary comparison also showed that the overall response pattern remained broadly comparable across the four years, while exact turning-point locations varied by year and scale (see supplementary material, Figure S1).

From the perspective of urban thermal environment management, our findings indicate that the cooling effect of vegetation does not occur in isolation. Instead, it is embedded within a thermal regulation system shaped by the combination of urban physical structure and land cover. Its performance is highly sensitive to spatial scale. Specifically, the NDVI–LST relationship is more stable at coarser scales and more context-dependent at finer scales, indicating that thermal management strategies should adopt differentiated configuration and design approaches across scales.

Based on the nonlinear NDVI-LST turning points identified in this study, these turning points can be used as scale-specific screening references to identify areas that remain below minimum effective greening levels under a given land use context. Accordingly, a hierarchical management strategy aligned with the MAUP multi-scale grid framework can be formulated:

• At the macro scale, the focus should be on the continuous layout of structural green spaces, including green belts and waterfront corridors, as well as ventilation corridors. Priority should be given to built-up areas with high heat exposure that remain below the screening thresholds, so that they can be targeted for green network reinforcement and cold-source connectivity. This approach enhances the stability of overall thermal regulation.

• At the meso scale, the thresholds can be used to screen communities where green coverage remains below the minimum effective level and where additional continuous vegetation cover and better connectivity are needed. In communities near or above this level, efforts should focus on optimizing green-space morphology and configuration, including connectivity, aggregation and boundary complexity, to improve cooling efficiency per unit of green coverage and avoid excessive fragmentation.

• At the micro scale, in areas characterized by stronger local heterogeneity, refined measures should be adopted to improve local cooling efficiency. These measures include street trees and shading, pocket parks, as well as vertical greening and green roofs. Such interventions should be coordinated with heat-adaptive design strategies, including high-albedo materials and permeable pavements.

This hierarchical strategy helps avoid policy inefficiencies and goal conflicts caused by scale mismatches. It also provides a reusable framework for cross-city threshold validation and scale transfer.

This study has several limitations. First, NDVI, as a two-dimensional indicator, cannot fully capture the thermal effects of three-dimensional vegetation structure, urban morphology or key meteorological variables such as wind speed. Second, because the observations are spatially structured, the ordinary random train/test split and conventional K-fold cross-validation used in this study may not fully remove spatial dependence. Third, although four periods from 2000 to 2024 were used to characterize changing urban backgrounds, and a brief year-specific supplementary comparison was conducted, systematic year-by-year quantification of threshold locations and response curves was not included in the main analytical framework. The findings should therefore be interpreted primarily as methodological and structural evidence of scale dependence rather than as proof of temporally invariant threshold stability. In addition, cooling capacity can vary across climate contexts, so the threshold-like responses identified in Shanghai should be interpreted as context-dependent rather than directly transferable (Ren et al., 2024). Future research could further test these responses across different years and with spatially explicit validation strategies.

This study takes Shanghai as a case study. Based on multi-temporal remote sensing data, it analyzes the relationship between NDVI and LST across different spatial scales and land use contexts from the perspective of the MAUP. The results show that the NDVI–LST relationship exhibits significant nonlinear characteristics and differences in stability across scales. Its response structure and key threshold features show clear scale dependency. Furthermore, the regulatory effect of vegetation varies across different land use types. This indicates that the thermal environment implications of NDVI are strongly dependent on spatial context.

From a management perspective, this study emphasizes that scale selection and land use context critically influence urban thermal environment assessment. The NDVI–LST relationship manifests differently across spatial scales. Therefore, thermal environment evaluation and cooling potential assessments need to align with the specific scale context. Overall, this study reveals how spatial scale and land use context jointly shape the nonlinear response structure of the vegetation and thermal environment relationship. It provides a theoretical basis for a scale-dependent understanding of urban thermal regulation mechanisms.

The supplementary material for this article can be found online.