Accurate spatial interpolation (SI) of climate data is vital for the management and supervision of natural resources and agriculture. Owing to the lack of an adequate number of meteorological stations, meteorological-station-data-based SI methods may not always reflect the real climatic conditions of an interpolated point. Land surface temperature (LST) data obtained from satellite sensors enable the characterization of meteorological conditions of areas without meteorological stations. The aim of this article is to present a new modified inverse distance weighting (M-IDW) SI method for air temperature (T-a), total precipitation (P-t), and relative humidity (RH) by integrating Landsat LST data with meteorological station data for the interpolation process. The M-IDW approach is based on the correlation relationship between the climate data and LST at each meteorological station, which is incorporated into the traditional IDW to improve the estimation of the climate data at an interpolation location of interest. The proposed method, M-IDW, is applied for the interpolation of long years' (i.e. long term) monthly average (LYMA) T-a, P-t, and RH climate data from meteorological stations in the Eastern Thrace region, which is 23,764 km(2), located in southeast Europe. The LYMA of the T-a, P-t, and RH has been constructed using data obtained from 27 meteorological stations that had functioned at least 10 years between 2000 and 2012 and from the corresponding satellite data. The outputs of the interpolation are in the form of LYMA, so are the analysed climate data. The spatial resolution of the predicted surface was taken as 30 m, similar to the original data presented by United States Geological Survey. The results were compared with those of the standard IDW, ordinary kriging (OK), and ordinary cokriging (OCK) methods to analyse the performance and accuracy of the proposed method. The results show that the proposed M-IDW method has the potential for SI of climate data, if enough number of images and cloudless pixels are incorporated in the LYMA LST computation. The proposed method, in general, yields better results than standard IDW and OK methods, especially during spring, summer, and partially in autumn for the interpolation of T-a (with 0.72 degrees C, 0.53 degrees C, and 0.66 degrees C root mean square error (RMSE) values, respectively) and Pt (with 11.07 mm, 7.64 mm, and 4.85 mm RMSE values, respectively). OCK and M-IDW results were comparable in spring, summer, and autumn where M-IDW was slightly better for T-a in autumn and spring and was slightly better for P-t in summer. For the RH interpolation, although M-IDW results were found to be close to the results of IDW, OK, and OCK in spring, summer, and autumn, for the overall seasonal interpretation, the RMSE values of M-IDW were worse than the others. In general, M-IDW yields worse results for the winter months, which in turn is related to cloudiness and availability of satellite images.