Mapping Urban Heat Risk and Climate-Adaptation Priorities in Saudi Cities: An Explainable Geospatial Machine-Learning Framework Using Open Satellite, Climate, and Population Data

The management of urban heat in Saudi Arabia needs decision-support methods to separate persistent thermal patterns from city-specific climate and morphology effects. The study presents a reproducible open-data framework for Riyadh, Jeddah, and Dammam based on observations during the summer months (June–September) from 2014 to 2025. We derived land surface temperature (LST), NDVI and NDBI from Landsat 8/9 Collection 2 Level 2 imagery, background meteorology from ERA5-Land, built-up surface and time-varying population from GHSL, and a 2020 age-structure vulnerability layer from WorldPop. Each city was studied within a standard 30 km radial study window with permanent water masked. The final database included 545,027 spatial records and a Landsat manifest of 720 scenes, with no missing values for the main analytical variables. Annual median summer LST P90 varied significantly at the city level (Kruskal–Wallis H = 20.20, p < 0.001). The highest 12-year median was observed in Dammam (56.80 °C) followed by Jeddah (54.90 °C) and Riyadh (53.66 °C). No statistically significant monotonic trend was detected in the annual spatial median of summer LST P90 in any of the three cities over 2014–2025. NDVI was negatively correlated with LST P90, while NDBI positively correlated with LST P90 in all three cities. Under 0.05° spatial-block cross-validation, histogram gradient boosting achieved R2 = 0.679 ± 0.044 and RMSE = 1.54 ± 0.06 °C. Because ERA5-Land predictors are spatially coarser than the validation blocks, these performance estimates may be optimistic for independently sampled locations. Leave-one-city-out tests showed restricted transferability under a two-city training design, potentially reflecting environmental distribution shifts and limited model extrapolation. Permutation importance and SHAP identified background air temperature, relative humidity, elevation, NDBI and population intensity as important predictors. A 2025 adaptation-priority index combining thermal hazard, population exposure, and age-related sensitivity showed similar rankings under geometric and arithmetic aggregation (Spearman ρ = 0.90–0.95), although this agreement does not constitute independent validation, and the resulting classes are strongly differentiated by population exposure. The results support city-specific and spatially targeted heat-resilience planning, rather than a single national thermal model.