Exploring the spatiotemporal patterns of urban thermal environments is crucial for mitigating the detrimental effects of urban heat islands (UHI). However, the long-term and fine-grained monitoring of UHI is limited by the temporal and spatial resolutions of various sensors. To address this limitation, this study employed the Google Earth Engine (GEE) platform and a multi-source remote sensing data fusion approach to generate a densely time-resolved Landsat-like Land Surface Temperature (LST) dataset for daytime observations spanning from 2001 to 2020 in Shanghai. A comprehensive analysis of the spatiotemporal patterns of UHI was conducted. The results indicate that over the past 20 years, the highest increase in average LST was observed during spring with a growth coefficient of 0.23, while the lowest increase occurred during autumn (growth coefficient of 0.12). The summer season exhibited the most pronounced UHI effect in the region (average proportion of Strong UHI and General UHI was 28.73%), while the winter season showed the weakest UHI effect (proportion of 22.77%). The Strong UHI areas gradually expanded outward over time, with a noticeable intensification of heat island intensity in the northwest and coastal regions, while other areas did not exhibit significant changes. Impervious surfaces contributed the most to LST, with a contribution of 0.96 °C, while water had the lowest contribution (−0.42 °C). The average correlation coefficients between LST and NDVI, NDWI, and NDBI over 20 years were −0.4236, −0.5128, and 0.5631, respectively.

A good understanding of the processes of land surface temperature (LST) change is important for assessing regional climate change. In the present study, we obtained the MODIS MOD11A2 LST products over the Yangtze River Delta (YRD) from 2001 to 2020. In order to comprehensively assess the spatial and temporal variability of LST in the YRD region over the past two decades, the Theil–Sen Median trend analysis and Mann–Kendall test, BFAST01 trend decomposition, and landscape pattern analysis were used in this study. We show that the rate of linear change in LST in the YRD ranges from −0.019 °C/month to 0.046 °C/month. The BFAST01 trend decomposition identifies more details of LST change and monotonic increases, reversal increase, and interruption increase are the main warming trends. The distribution of the different trend types shows strong aggregation with high spatial heterogeneity. The LST breakpoints are mainly located in the northern and southern YRD, which frequently occurred during 2010–2013. Of the various land types, breakpoints occur most frequently in cropland and high NDVI (0.5–0.7) areas, and the intensity of most of them is within 2 °C. In addition, much stronger warming occurs in urban areas than in other land types. Our study provides a better understanding of the dynamics of LST in the YRD region over the past 20 years and highlights that breakpoints cannot be circumvented in regional temperature assessment.