The commonly used weather stations cannot fully capture the spatiotemporal variability of near-surface air temperature (T(air)), leading to exposure misclassification and biased health effect estimates. We aimed to improve the spatiotemporal coverage of T(air) data in Germany by using multi-stage modeling to estimate daily 1 × 1 km minimum (T(min)), mean (T(mean)), maximum (T(max)) T(air) and diurnal T(air) range during 2000-2020. We used weather station T(air) 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 T(air) from cells with both T(air) and LST available. In the second stage, we used this model to predict T(air) for cells with only LST available. In the third stage, we regressed the second stage predictions against interpolated T(air) values to obtain T(air) countrywide. All models achieved high accuracy (0.91 ≤ R(2) ≤ 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 ≤ R(2) ≤ 0.99, 0.87 °C ≤ RMSE ≤ 2.05 °C) and countrywide, for the T(mean) model (0.71 ≤ R(2) ≤ 0.99, 0.79 °C ≤ RMSE ≤ 1.19 °C). Annual T(mean) 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 T(air) variability representation for exposure assessment of participants in health cohorts. Our results indicate the proposed models as suitable for estimating nationwide T(air) at high resolution. Our product is critical for temperature-based epidemiological studies and is also available for other research purposes.