Spatiotemporal Hotspot Analysis of Forest Thermal Variability and Vegetation Stability Using Time Series LST and NDVI -Focusing on the Namyangju, Republic of Korea

1. Introduction

Forests identified in land-cover maps are increasingly recognized as a Nature-based Solution (NBS) that moderates climate through tree-dominated processes [1,2,3]. In particular, forest canopies buffer energy exchange between the atmosphere and the land surface and regulate surface temperature via radiative shielding and evapotranspiration, thereby playing a pivotal role in maintaining local climate stability [4,5,6]. Forests also store atmospheric carbon dioxide (CO₂) in biomass, reducing greenhouse gases (GHGs) and mitigating climate change. These functions position forests, within the NBS framework, as critical natural capital that delivers both mitigation and adaptation [7,8,9,10]. Sustaining these NBS benefits requires continuous, wall-to-wall monitoring and adaptive forest management. [11,12,13]. However, approximately 63% of South Korea is forested, making comprehensive field surveys across the entire territory impractical [14,15]. To overcome this constraint, the forest sector increasingly relies on remote sensing for synoptic, repeatable, and spatially explicit monitoring beyond sample-based assessments.

To operationalize monitoring of these NBS functions, we employ remote-sensing indicators that jointly capture canopy condition and surface thermal states—most notably the normalized difference vegetation index (NDVI) and land surface temperature (LST). The NDVI is a canonical observation-derived index that quantifies vegetation activity and cover [16,17]. Generally, higher NDVI indicates denser and healthier vegetation, supporting observation-based monitoring [18,19,20]. Land Surface Temperature (LST) is derived from thermal-infrared emissive properties to estimate surface temperature [21,22]. LST enables monitoring of thermal-related phenomena such as urban heat islands, drought, and crop or vegetation stress [23,24,25]. Both NDVI and LST can be obtained from satellite observations at commensurate acquisition times, for example from the United States Geological Survey (USGS) Landsat series [26,27]. NDVI and LST are tightly coupled: in summer, areas with high NDVI typically exhibit lower LST because canopy shading and evapotranspiration depress surface temperatures, yielding a negative association; in winter, the relationship can invert and become positive [28,29,30,31,32]. These seasonal reversals, driven by acquisition timing and phenology, are a critical consideration for data use [18,33,34,35], especially in strongly seasonal climates such as South Korea, where phenological dynamics render the NDVI–LST association non-stationary throughout the year.

Accordingly, treating the NDVI–LST relationship as a single global correlation is inadequate; approaches that reflect spatiotemporal structure and variability are required. In forested landscapes, both variables exhibit pronounced spatial autocorrelation, with similar values clustering among neighboring locations [36,37]. Conventional regression models that ignore this spatial structure can yield biased or misleading inferences, and mean-based statistics alone are ill-suited to capture such patterns. To characterize the spatial distribution of NDVI and LST more precisely, spatial statistics, particularly hotspot analysis capable of identifying local clusters, are warranted [38,39]. While hotspot analysis is informative, it is typically univariate, and thus additional procedures are needed to incorporate temporal variation.

This study addresses those limitations by integrating hotspot analysis with time-series statistics to elucidate the spatiotemporal relationship between LST and NDVI and, in turn, to explore forests’ climate-buffering function. We compute the annual coefficient of variation (CV) for time-series NDVI and LST and conduct hotspot analysis on these CV fields using Getis–Ord Gi* hotspot analysis [40,41]. The CV, a dimensionless ratio of the annual standard deviation to the mean, expresses fluctuation amplitude relative to the mean and thereby enables direct comparison of the relative variability of two variables with different units. We then assess their association by regressing the Getis–Ord Gi* Z-score (GiZScore) fields. In addition, we designate the upper 5% of the GiZScore distribution as “management-priority sites” and derive policy implications for spatially explicit forest management. Beyond purely numeric summary statistics, this framework provides a basis for understanding, in space and time, how forests respond to climate variability and for designing targeted responses.

This study offers decision-relevant evidence by analyzing NDVI–LST interactions across space and time and by objectively identifying forest areas where climate-mitigation and buffering functions are most pronounced. By overcoming limitations of conventional spatial analyses and leveraging simple transformations of readily available inputs, the approach proposes a new mode of using satellite-derived data with strong temporal dependence in seasonal regions such as South Korea.

2. Materials and Methods

2.1. Study Area

The study area is Namyangju, Gyeonggi-do, Republic of Korea (Figure 1). Situated on the eastern periphery of the Seoul Capital Region, Namyangju is an administrative district in which approximately 70% of the total area is forested, substantially higher than the provincial average of 52% [42]. This high forest fraction provides favorable conditions for evaluating the ecological and climatic functions that forests perform in urban-adjacent settings. Owing to its mosaic of urban fabric interspersed with forests, Namyangju is well suited to observe how multiple forest functions, including climate regulation, mitigation of urban thermal environments, and the provision of ecosystem services, are connected to everyday living spaces.

Additionally, Namyangju features steep slopes and mountainous terrain spanning a wide elevational range, making it well suited to analyze topographically driven spatial heterogeneity in vegetation and the attendant variation in land surface temperature. The area exhibits a mid-latitude continental climate with pronounced seasonal temperature contrasts, allowing observation of both summer heat and winter cold, which is advantageous for elucidating the seasonal relationship between LST and NDVI [43,44,45]. Among the relevant administrative units from which cloud-free, multi-temporal Landsat-8 and Landsat-9 imagery can be acquired to construct the time-series LST and NDVI inputs, including the neighboring, heavily forested jurisdictions of Pocheon and Yangpyeong, only Namyangju met the data-availability criteria for this study.

2.2. Datasets for Analysis

Both LST and NDVI are canonical observation-derived indices computed from satellite imagery. The two metrics respectively reflect the thermal state of the land surface and the vigor of vegetation activity. LST is retrieved from thermal-infrared radiance, typically by inverting a radiative transfer/Planck function to brightness temperature and then converting to land surface temperature with emissivity correction [46]. NDVI, defined from the contrast between red and near-infrared reflectance, captures vegetation density and condition [47]. Meaningful interpretation is enabled by comparing these datasets at coincident acquisition times and harmonized spatial resolution. Accordingly, quantifying the LST–NDVI relationship requires image data acquired by the same (or cross-calibrated) sensor platform. In this context, we used imagery from the Landsat-8 and Landsat-9 missions provided by the USGS.

The two satellites carry nearly identical and cross-calibrated payloads, Operational Land Imager (OLI/OLI-2) and Thermal Infrared Sensor (TIRS/TIRS-2), yielding high data continuity and interoperability [48]. OLI supplies the visible and near-infrared bands used to compute NDVI, while TIRS provides the thermal bands needed to estimate LST. Spatial sampling is 30 m for OLI, and 100 m (delivered resampled to 30 m) for TIRS, with each satellite on a 16-day repeat cycle; phased orbits provide ~8-day effective revisit when combined. For this study, we assembled eight cloud-free, multi-temporal Landsat-8/9 scenes over Namyangju during 2023; the acquired data are provided as Figure 2.

2.3. Data Processing

LST was retrieved from the thermal-infrared Band 10 of Landsat-8 and Landsat-9 Level-1 scenes, and NDVI was computed from the red (Band 4) and near-infrared (Band 5) bands. Although the USGS Level-2 product provides an LST layer, we found that the accompanying cloud mask was not adequate for our purposes and led to loss of valid observations across the time series. Accordingly, we derived both LST and NDVI directly from the Level-1 inputs.

To quantify the temporal stability and variability of LST and NDVI within forests, we spatially masked satellite scenes for Namyangju using the Ministry of Environment’s 2023 micro class land-cover map and retained only areas classified as “forest.” [49]. This restriction was intended to isolate and analyze climate–vegetation interactions specifically within forested landscapes.

Over the extracted forest mask, we generated a regular lattice of grid points at 30-m spacing, co-registered to the imagery and consistent with the native spatial resolution, yielding a total of 307,886 points. For each point, we extracted the complete time-series of LST and NDVI values for the study period. The resulting per-point series were then used to compute the mean, standard deviation (SD), and CV for each variable.

The CV, defined as SD divided by the mean, is a dimensionless measure that enables direct comparison of relative variability across variables with different units and ranges. Reliance on SD alone can over- or under-state variability as a function of absolute magnitude; introducing the CV mitigates this bias and provides a more interpretable basis for comparing the relative stability or sensitivity of LST and NDVI within forest areas.

We next performed hotspot analysis on the CV fields of LST and NDVI separately, using the Getis–Ord Gi* statistic and its standardized Z-score (GiZScore) to identify clusters of significantly high (hot) and low (cold) variability[40,41,50]. To assess the spatial coupling of variability between the two indices, we regressed the NDVI GiZScores on the LST GiZScores. The data flow of the study is summarized in Figure 3. Based on the hypothesis that areas with high thermal variability also experience elevated vegetation variability, we expected a positive regression coefficient. This analysis evaluates the strength and significance of that hypothesis at spatial scales relevant to management, providing evidence on how the forest climate-buffering function varies across space. In doing so, it empirically links thermal sensitivity and vegetation stability and offers spatially explicit information to support climate-adaptation and forest-management decision-making.

3. Results

3.1. Result of Getis–Ord Gi* Hotspot Analysis for CVs of LST and NDVI

The spatial distributions of hotspots in annual CV of LST and NDVI, derived via Getis–Ord Gi* hotspot analysis, are shown in Figure 4a, b. In both variables, the Getis–Ord Gi* maps of annual CV reveal a coherent spatial organization. Despite differences in the areal extent of clusters across confidence thresholds, LST and NDVI exhibit broadly congruent patterns: hotspots form continuous bands along high-elevation ridges and dissected mountain slopes, whereas coldspots dominate valley floors and low-lying terrain. The close correspondence with the digital elevation model (DEM) indicates that terrain imposes a first-order control on intra-annual variability (Figure 1b).

We delineated management-priority areas by selecting grid cells that simultaneously fell within the top 5% of Getis–Ord Gi* Z-scores for both LST CV and NDVI CV—i.e., locations where intra-annual thermal and greenness variability are each significantly clustered at or above the 95th percentile. Operationally, we intersected the two hotspot layers to isolate zones of co-clustered variability. Comparison with the DEM (Figure 1b) shows that these co-hotspot belts are preferentially situated in higher-elevation terrain, with elevation distributions shifted upward relative to the forest-wide baseline and spatial concentrations along upland ridges and dissected, steep slopes. This topographic alignment indicates that orographic controls—enhanced thermal contrasts, aspect-driven insolation, wind exposure, and thinner soils—coincide with stronger phenological amplitude and canopy structural heterogeneity, jointly amplifying both LST and NDVI variability. From a management standpoint, these high-elevation co-hotspots constitute priority surveillance and intervention zones for climate-adaptation actions.

In mountainous zones, stronger seasonal thermal contrasts, aspect-driven insolation gradients, wind exposure, thinner soils, and intermittent water limitation jointly amplify LST variability; concomitantly, phenological amplitude and canopy structural heterogeneity enhance NDVI variability. Conversely, valley bottoms—benefiting from hydrologic buffering, gentler slopes, and more continuous canopy cover—show relatively muted variability in both fields.

A few local mismatches between the two maps suggest partial decoupling of thermal and greenness dynamics. For example, LST hotspots extend more broadly than NDVI hotspots in some ridge–shoulder complexes, consistent with thermal variability not always translating into proportional changes in vegetation greenness (e.g., where evergreen cover, shading, or soil moisture mitigate phenological swings). Taken together, the patterns support our hypothesis that spatial clustering of thermal variability tends to co-occur with clustering of vegetation variability, especially in complex terrain, and they identify upland belts as priority zones for monitoring climate–vegetation sensitivity.

3.2. Quantification of LST–NDVI Linkage Using Ordinary Least Squares

We estimated two ordinary least squares (OLS) specifications to quantify the LST–NDVI linkage across 307,886 forest grid points: (i) a variability model relating the coefficients of variation (CV–CV) and (ii) a hotspot model relating the Getis–Ord Gi* Z-scores (GiZ–GiZ), in each case regressing NDVI on LST. The CV specification yielded NDVI_CV = 0.982 × LST_CV − 0.066 with R² = 0.464 (Pearson r = 0.681), indicating a moderate–strong positive association in relative variability. The hotspot specification produced NDVI_GiZ = 0.615 × LST_GiZ + 0.025 with R² = 0.623 (Pearson r = 0.789), demonstrating a substantially tighter fit and stronger spatial concordance. Rank-based (Spearman) correlations showed the same ordering, confirming that the advantage of the hotspot formulation is not driven by a few extreme observations (Figure 5a, b). Interpreting the slopes, a one-unit increase in the LST GiZScore (i.e., one standard-unit increase in local clustering intensity of thermal variability) is associated with a ~0.62-unit increase in the NDVI GiZScore, consistent with co-located clusters of thermal and greenness variability (Figure 4a, b).

To facilitate direct visual comparison of the two relationships, we standardized each variable to z-scores (mean 0, standard deviation 1) and overlaid the CV–CV and GiZ–GiZ pairs in a single panel with identical axes. Specifically, both axes in Figure 5c are constrained to [−2, 5] and use a 1:1 scale; CV–CV is plotted in orange and GiZ–GiZ in dark green, with OLS lines drawn over point clouds. This normalization removes unit and scale effects, allowing the slopes to be interpreted purely as standardized effect sizes. Consistent with the regression results, the GiZ–GiZ line exhibits a steeper standardized slope than the CV–CV line, visually reinforcing that local clustering of thermal variability co-occurs more strongly with local clustering of vegetation variability.

We regard the GiZ-based formulation as the more appropriate representation of the LST–NDVI relationship. GiZScores are standardized local cluster statistics that explicitly incorporate neighborhood structure through a spatial weights matrix W. This construction reduces sensitivity to pixel-level noise (e.g., residual cloud shadow, geolocation jitter) and mitigates variance heterogeneity induced by spatial non-stationarity (e.g., terrain-driven shifts in local means and variances). As a result, the GiZ–GiZ regression yields a clearer and more stable signal of spatial co-clustering—the phenomenon of primary interest in this study—than a regression on raw CVs. Put differently, whereas CV–CV models can be confounded by large-scale gradients (elevation, aspect, soil moisture) that inflate or depress variability in different subregions, the GiZ transformation benchmarks each location against its local context, thereby sharpening the inferential contrast between hot and cold belts of joint variability.

For these reasons, we adopt the GiZ–GiZ OLS as our primary indicator and interpret it as a statistical characterization of the spatial association between thermal and greenness dynamics. We did not fit additional spatial-lag or spatial-error regressions because the GiZ transformation already embeds key spatial features, local averaging, standardization, and neighbor weighting, into the dependent and independent variables. This parsimonious approach avoids over-parameterization while directly targeting the mechanism of interest: the co-location of hotspots and coldspots in LST and NDVI.

4. Discussion

This study demonstrates that spatiotemporal clustering of thermal variability (LST CV) is closely mirrored by clustering of vegetation variability (NDVI CV) across a mountainous, forest-dominated urban fringe. The Getis–Ord Gi* results show that hotspots of intra-annual variability align with high-elevation ridges and dissected slopes, whereas coldspots concentrate in valley floors. This concordance with topography indicates that terrain imposes a first-order control on both thermal dynamics and phenology: stronger seasonal thermal contrasts, aspect-driven insolation gradients, wind exposure, and thinner soils in uplands amplify LST fluctuations; simultaneously, phenological amplitude and canopy structural heterogeneity increase NDVI variability. Conversely, hydrologically buffered lowlands with gentler slopes exhibit muted variability in both fields. Together, these patterns support the interpretation of forests as climate buffers whose effectiveness varies systematically with terrain and landscape position [51,52,53].

A key contribution of this study is the comparison between two ways of formalizing the LST–NDVI linkage using observation based remote sensing data: a direct regression of variability magnitudes (CV–CV) and a regression of local clustering intensities (GiZ–GiZ) [29,54,55]. Although both specifications yielded positive and statistically meaningful associations, the hotspot formulation provided a substantially tighter fit (R² = 0.623 vs. 0.464) and stronger rank concordance. We attribute this improvement to the properties of GiZScores: as standardized local statistics that explicitly embed neighborhood structure via a spatial weights matrix, GiZScores benchmark each pixel against its neighbors, attenuating pixel-scale noise and reducing variance heterogeneity caused by spatial non-stationarity[56,57]. In practical terms, the GiZ–GiZ model targets the phenomenon of interest—the co-location of thermal and greenness hotspots/coldspots—more directly than CV–CV and therefore serves as a more appropriate indicator of spatial coupling in complex terrain [58,59].

These findings carry several implications for forest monitoring and management in urban-adjacent mountain regions. First, the co-clustering of high thermal and high greenness variability in uplands suggests that upland belts function as sensitivity zones where climate anomalies are more likely to propagate into pronounced phenological responses [60,61,62]. Second, because our GiZ approach identifies statistically concentrated belts rather than isolated pixels, it can inform operational prioritization—for example, targeting the top 5% of GiZScores as “management-priority sites” for intensified field checks, fuel-load management, microclimate restoration (e.g., windbreaks, understory moisture conservation), or the placement of in situ temperature/soil-moisture loggers. Third, the persistence of coldspots in valley floors highlights the role of hydrological buffering; safeguarding riparian corridors and cold-air drainage pathways may help maintain local thermal stability and phenological regularity within the broader urban matrix.

Methodologically, the study contributes a parsimonious yet spatially informed pipeline for coupling time-series variability metrics with local cluster statistics. Computing annual CVs provides a unitless, comparable measure of intra-annual fluctuation; applying Gi* to the CV fields transforms those fluctuations into statistically vetted clusters; and a simple OLS on GiZScores yields an interpretable slope (co-clustering strength) without resorting to heavily parameterized spatial-lag or spatial-error models [63,64]. We deliberately did not fit SAR/SEM/SDM models because the GiZ transformation already embeds key spatial features—local averaging, standardization, and neighbor weighting—into both the predictor and response, thereby mitigating spatial dependence at the modeling stage while keeping the analysis transparent and reproducible. In settings where explicit process-based spatial spillovers are hypothesized (e.g., fire spread), spatial regression would remain valuable; however, for co-clustering of variability, GiZ–GiZ OLS strikes a favorable balance between robustness and simplicity.

This study has four main limitations. (1) Temporal scope—analyses are limited to a single year (2023), so interannual stability and trends remain to be tested. Additionally, satellite imagery from the summer season (June to August) was unavailable. This is because, due to the extreme cloud cover characteristic of Korea's climate, images with no cloud cover were excluded during the screening process. Therefore, strictly speaking, this does not cover all phenolgy for the entire year. (2) LST retrieval uncertainty—Level-1–based thermal estimates may retain emissivity and atmospheric correction errors and warrant in situ validation. (3) Index choice—NDVI saturation can limit sensitivity in dense canopies; alternative or complementary indices (e.g., Enhanced vegetation index (EVI), red-edge, Solar-induced chlorophyll fluorescence (SIF)) should be evaluated. (4) Neighborhood design—GiZ outcomes depend on the spatial weights matrix and neighborhood size. Future work can address these by extending to multi-year datasets, refining LST retrievals with ground truth, testing additional vegetation indices, and conducting sensitivity and multiple-testing assessments.

From a policy perspective, integrating GiZ-based hotspot belts into municipal forest plans could help align routine surveillance with areas of greatest climate–vegetation sensitivity. Because the method is sensor-agnostic and computationally light, it is readily transferable to other counties in Gyeonggi-do and to other seasonal, mountainous regions where rapid, wall-to-wall screening is needed. More broadly, mapping the spatial concordance of thermal and phenological variability provides an observation-based line of evidence for climate-buffering by forests, supporting nature-based solutions in local climate-adaptation and carbon-neutrality strategies [65,66].

In sum, by combining time-series variability metrics with standardized local clustering, this study provides a statistically grounded yet operationally usable view of how forests modulate microclimate across complex terrain. The stronger performance of the GiZ formulation underscores the value of explicitly leveraging neighborhood structure when translating pixel-level satellite observations into management-relevant insights about spatial coupling between thermal regimes and vegetation dynamics.

5. Conclusions

This study shows that intra-annual thermal variability (LST CV) and vegetation variability (NDVI CV) co-cluster systematically in mountainous, forested urban fringes, with hotspots aligned to upland ridges and coldspots concentrated in valleys. By pairing time-series variability metrics with Getis–Ord Gi* hotspot statistics, we quantified not only the magnitude of variability but also its spatial organization. Regressions confirmed a positive LST–NDVI linkage in both formulations, with the GiZ–GiZ specification yielding a substantially tighter fit than CV–CV, indicating that local clustering of thermal dynamics more strongly co-occurs with clustering of greenness dynamics than raw variability alone would suggest.

Methodologically, the workflow—annual CV computation → Gi* hotspot mapping → OLS on GiZScores—offers a transparent, lightweight, and spatially informed alternative to heavily parameterized spatial models. Because GiZScores embed neighborhood structure via the spatial weights matrix, they mitigate pixel-scale noise and variance heterogeneity, producing an interpretable slope that directly indexes co-clustering strength. The standardized overlay (Figure 3) further demonstrates, on a common scale, the stronger standardized slope for GiZ–GiZ, reinforcing the spatial coupling between thermal and phenological variability.

Practically, these results support using GiZ-based hotspot belts to prioritize monitoring and intervention in upland “sensitivity zones,” while recognizing valley coldspots as hydrologically buffered areas that stabilize local microclimate. The approach is sensor-agnostic, computationally efficient, and readily transferable to other seasonal, mountainous regions, making it suitable for routine, wall-to-wall screening in municipal forest plans and for integrating observation-based evidence into nature-based climate-adaptation and carbon-neutrality strategies.

Future work should generalize beyond a single year, refine thermal retrievals with ground validation, and test complementary vegetation indices (e.g., EVI, red-edge, SIF) to probe canopy dynamics more sensitively. Even so, the present findings provide decision-relevant evidence that explicitly leveraging neighborhood structure clarifies where and how forests buffer climate across complex terrain—turning pixel-level satellite signals into spatial intelligence for management.