Spatial climate datasets currently available for Bhutan are limited by weather station data availability, spatial resolution or interpolation methodology. This article presents new datasets for monthly maximum temperature, minimum temperature, precipitation and vapour pressure climate normals interpolated for the 1986–2015 reference period using trivariate smoothing splines. The inclusion of standardized day time Moderate Resolution Imaging Spectroradiometer (MODIS) land surface temperature (LST) values as partial spline dependencies reduced cross validated root mean square error (RMSE) for maximum temperature by up to 16.0% and was most effective between March and September. Using both a topographic index of relative elevation and standardized night time MODIS LST values as partial spline dependencies reduced monthly mean minimum temperature RMSE by up to 23.4%. Neither variable was effective for minimum temperature interpolation between June and September. High humidity, extensive cloud cover and heavy precipitation occur during these months, which are likely to suppress the formation of temperature inversions that typically form under clear, calm conditions. These new temperature and precipitation surfaces show distinct differences from the WorldClim and CRU CL 2.0 datasets, which do not use weather stations within Bhutan for calibration. New precipitation surfaces better describe the heavy rainfall experienced in the southern foothills while retaining the effect of orography throughout the central valleys and ranges. The development of vapour pressure surfaces also allow for the calculation of ecologically important variables such as vapour pressure deficit, and may also be useful for solar radiation modelling in the region. The different datasets presented in this article will facilitate ecological and agricultural research in Bhutan and provide high quality surfaces needed for future climate change scenarios.