The commonly used weather stations cannot fully capture the spatiotemporal variability of near-surface air temperature (Tair), leading to exposure misclassification and biased health effect estimates. We aimed to improve the spatiotemporal coverage of Tair data in Germany by using multi-stage modeling to estimate daily 1 × 1 km minimum (Tmin), mean (Tmean), maximum (Tmax) Tair and diurnal Tair range during 2000–2020. We used weather station Tair observations, satellite-based land surface temperature (LST), elevation, vegetation and various land use predictors. In the first stage, we built a linear mixed model with daily random intercepts and slopes for LST adjusted for several spatial predictors to estimate Tair from cells with both Tair and LST available. In the second stage, we used this model to predict Tair for cells with only LST available. In the third stage, we regressed the second stage predictions against interpolated Tair values to obtain Tair countrywide. All models achieved high accuracy (0.91 ≤ R2 ≤ 0.98) and low errors (1.03 °C ≤ Root Mean Square Error (RMSE) ≤ 2.02 °C). Validation with external data confirmed the good performance, locally, i.e., in Augsburg for all models (0.74 ≤ R2 ≤ 0.99, 0.87 °C ≤ RMSE ≤ 2.05 °C) and countrywide, for the Tmean model (0.71 ≤ R2 ≤ 0.99, 0.79 °C ≤ RMSE ≤ 1.19 °C). Annual Tmean averages ranged from 8.56 °C to 10.42 °C with the years beyond 2016 being constantly hotter than the 21-year average. The spatial variability within Germany exceeded 15 °C annually on average following patterns including mountains, rivers and urbanization. Using a case study, we showed that modeling leads to broader Tair variability representation for exposure assessment of participants in health cohorts. Our results indicate the proposed models as suitable for estimating nationwide Tair at high resolution. Our product is critical for temperature-based epidemiological studies and is also available for other research purposes.