Urban Heat Island (UHI) effects are among the most visible climatic consequences of rapid urbanization, driven by the expansion of impervious surfaces, increasing building density, and changes in vegetation cover. Satellite remote sensing provides consistent spatial coverage and repeated observations over long time periods, making it a valuable tool for analyzing these processes. In particular, land surface temperature (LST) derived from multispectral satellite missions enables systematic monitoring of thermal conditions in urban environments. Satellite observations have therefore become an important source of information for examining how urban thermal patterns evolve over time.
However, interpreting multi-year satellite-derived land surface temperature time series is not always straightforward. LST observations often show strong interannual variability influenced by meteorological conditions, seasonal differences in data availability, and changes in land cover characteristics. Distinguishing a stable long-term urban warming signal from short-term fluctuations therefore remains a methodological challenge in many urban heat island studies. This difficulty becomes particularly relevant when relatively short satellite records are analyzed, as individual anomalous years can influence the interpretation of long-term trends.
This study investigates how stable UHI signals can be detected in multi-year satellite-derived land surface temperature time series using a reproducible open-source geospatial workflow. The analysis is demonstrated for the city of Zagreb, Croatia, using a ten-year satellite record (2015–2024) derived from Landsat observations (Landsat 8/9). Instead of analyzing the entire urban area as a single unit, the study focuses on several representative urban neighborhoods used as analytical units for evaluating temporal patterns in the temperature time series. Using several analytical units also allows temporal patterns to be compared across different parts of the city rather than relying on a single aggregated urban value.
Zagreb represents a typical Central European urban environment characterized by a mixture of dense built-up areas, residential neighborhoods, urban green spaces, and peri-urban zones. To capture this diversity, several representative neighborhoods were selected as analytical units. The selected areas reflect different urban development patterns within the city, including contrasts in vegetation cover and building density. This allows the analysis to compare thermal behavior across distinct types of urban environments within the same metropolitan area. Such variation in land cover and urban structure provides a useful setting for examining how temperature dynamics differ between neighborhoods with different surface characteristics.
Annual warm-season composites of land surface temperature are generated to reduce noise associated with individual satellite scenes and to provide a consistent representation of thermal conditions for each observation year. To place temperature dynamics in the context of surface characteristics, two commonly used spectral indices are analyzed alongside LST: the Normalized Difference Vegetation Index (NDVI) and the Normalized Difference Built-up Index (NDBI). Analyzing these indices together with temperature values helps interpret whether observed thermal changes are associated with variations in vegetation or built-up surfaces.
Satellite data processing follows a reproducible workflow that integrates image preprocessing, time series aggregation, and statistical analysis. Individual Landsat scenes are first filtered according to cloud cover and seasonal criteria. Annual warm-season composites are then generated to reduce scene-level variability and provide consistent yearly observations. LST, NDVI, and NDBI layers are subsequently extracted for each selected neighborhood, forming the basis for the temporal analysis. Using a consistent processing workflow also helps ensure that results from different years remain comparable.
Particular attention is given to the temporal behavior of the data through the analysis of interannual differences between consecutive observations. This allows year-to-year variations in land surface temperature to be examined together with corresponding changes in vegetation and built-up surfaces. Interannual differences help reveal short-term dynamics and provide additional insight into whether temperature changes occur in parallel with changes in surface characteristics. These comparisons provide an additional perspective on the variability present in the satellite record.
To assess the robustness of the detected trends, several analytical approaches are used. Linear trends and non-parametric slope estimates are calculated for the LST time series, and the results are compared across different parts of the observation period. This makes it possible to see whether similar trends appear regardless of the method used or the selected time interval. In this way, it can be evaluated whether the observed temperature signal reflects a persistent long-term pattern or mainly short-term variability in the satellite record. Examining the consistency of trend estimates therefore provides an additional check on the stability of the detected patterns.
The results help interpret satellite-derived temperature time series in urban heat island research more carefully. Examining temporal changes, interannual differences, and the consistency of trend estimates helps identify whether similar warming patterns appear across different parts of the time series. This comparison makes it possible to separate persistent temperature trends from short-term variability in the satellite observations.
Such an analysis also shows that time plays a key role when interpreting urban heat island dynamics. The use of reproducible open-source geospatial tools allows the same workflow to be applied to satellite temperature time series in other urban areas. This makes the approach potentially applicable to other cities where comparable satellite datasets are available.