High-resolution air temperature data is indispensable for analysing heatwave-related non-accidental mortality. However, the limited number of weather stations in urban areas makes obtaining such data challenging. Multi-source data fusion has been proposed as a countermeasure to tackle such challenges. Satellite products often offered high spatial resolution but suffered from being temporally discontinuous due to weather conditions. The characteristics of the data from reanalysis models were the opposite. However, few studies have explored the fusion of these datasets. This study is the first attempt to integrate satellite and reanalysis datasets by developing a two-step downscaling model to generate hourly air temperature data during heatwaves in London at 1 km resolution. Specifically, MODIS land surface temperature (LST) and other satellite-based local variables, including normalised difference vegetation index (NDVI), normalized difference water index (NDWI), modified normalised difference water index (MNDWI), elevation, surface emissivity, and ERA5-Land hourly air temperature were used. The model employed genetic programming (GP) algorithm to fuse multi-source data and generate statistical models and evaluated using ground measurements from six weather stations. The results showed that our model achieved promising performance with the RMSE of 0.335 °C, R-squared of 0.949, MAE of 1.115 °C, and NSE of 0.924. Elevation was indicated to be the most effective explanatory variable. The developed model provided continuous, hourly 1 km estimations and accurately described the temporal and spatial patterns of air temperature in London. Furthermore, it effectively captured the temporal variation of air temperature in urban areas during heatwaves, providing valuable insights for assessing the impact on human health.