Intensifying Cooling Service of Vegetation Under Urban Warming: A Multi-Temporal Google Earth Engine Analysis of the NDVI–Land Surface Temperature Relationship in Rajshahi City, Bangladesh (2001–2025)

Rapid urbanisation in mid-sized South Asian cities has intensified the surface urban heat island (SUHI), yet long-term, satellite-based assessments of the vegetation–temperature relationship for Rajshahi — one of the hottest cities in Bangladesh — remain scarce. This study quantifies the spatial and temporal relationship between vegetation health and land surface temperature (LST) across a 546 km2 study area centred on Rajshahi City for four pre-monsoon periods (2001, 2009, 2017 and 2025) using cloud-based Google Earth Engine (GEE) processing of Landsat 5, 8 and 9 Collection-2 surface-reflectance imagery. The Normalised Difference Vegetation Index (NDVI) and LST were derived per pixel at 30 m resolution; permanent water was masked using the Normalised Difference Water Index (NDWI), and a physically-based quality-control (QC) filter (LST ≥ 20 °C) was applied uniformly to remove residual cloud and cloud-shadow contamination. Pearson and Spearman correlation, ordinary-least-squares (OLS) regression, non-parametric bootstrap confidence intervals and Fisher r-to-z tests were applied to stratified random pixel samples (n ≈ 9,500–9,700 per year). A consistent, highly significant negative NDVI–LST correlation was found in 2001 (r = −0.583), 2009 (r = −0.774) and 2025 (r = −0.609), all at p < 0.001, confirming the cooling effect of vegetation. The strongest gradient occurred in the drought year 2009. Critically, the apparent positive correlation reported for the cloud-affected 2017 pre-monsoon (raw r = +0.359) was shown to be an artefact of ~250 contaminated low-temperature pixels: after QC, the relationship collapsed to statistical insignificance (r = −0.008, p = 0.42), demonstrating a climate-driven breakdown rather than a true inversion. The OLS cooling slope steepened from −10.82 °C per NDVI unit in 2001 to −12.77 °C per NDVI unit in 2025 — an 18.1 % intensification that a Fisher r-to-z test confirms is statistically significant (z = 2.80, p = 0.005). Concurrently, the share of land experiencing extreme surface heat (≥ 40 °C) rose from 1.1 % in 2001 to 8.0 % in 2025. These results provide quantitative evidence that the marginal climate-regulating value of vegetation is increasing as Rajshahi warms, and underscore the urgency of protecting the city’s remaining mango-orchard belt and peri-urban cropland as a measurable, low-cost climate-adaptation strategy.