The performance of three machine learning (ML) methods; cubist regression (CR), random forest (RF) and generalized additive model using splines (GAM) in generating monthly air temperature grids over Costa Rica was evaluated against two heavily used geostatistical methods; ordinary kriging (OK) and kriging with external drift (KED). The skill of the interpolation methods was evaluated using a 10-fold cross-validation technique; selecting the root-mean square error (RMSE), the mean absolute error (MAE) and the Pearson correlation-coefficient (R) as agreement metrics. To this purpose, data from an irregularly-distributed observational-network comprised by 73 weather-stations were selected for the period 1950-1987. Several spatial fields derived from a high-resolution digital elevation model (DEM) were tested as covariants. Results from the 10-fold cross-validation test show that CR yielded the best individual performance followed by KED, whereas GAM performed worst. Elevation on the other hand, was the only covariant ultimately incorporated in the interpolation process, since the remaining spatial fields exhibited poor correlation with temperature or resulted in data redundancy. While the quantitative and qualitative evaluation of CR and KED can be said to be comparable, CR is considered the best approach since the method is unaffected by assumptions on data normality and homoscedasticity.