<p>Understanding the spatial variability of precipitation is essential for water resource management and climate adaptation, especially in arid and semi-arid regions with strong spatiotemporal heterogeneity. Traditional geostatistical methods, such as ordinary kriging, often struggle to capture nonlinear relationships between rainfall and spatial coordinates. This study focuses on comparing ML–RK methods for spatial interpolation using only latitude and longitude as predictors, rather than developing a full rainfall prediction model. As, machine learning techniques integrated with regression kriging (RK) have wide applications for capturing complex spatial patterns. Therefore, this study evaluates RK combined with six regression models including Random Forest (RF), Support Vector Machine (SVM), K-Nearest Neighbors (KNN), Neural Network (NN), Elastic Net (EN), and Polynomial Regression (PR). In this research, we used monthly and decadal averages of precipitation from 42 meteorological stations (2001–2021) of Pakistan. For assessing optimal spatial structure, four theoretical variogram models including exponential, circular, spherical, and linear–were tested using Leave-One-Out Cross-Validation. Here, the performance of the variogram was assessed using RMSE and MAE. Outcomes associated with this research show that RF-RK consistently outperformed other combinations of ML-RK. Consequently, the combination of ensemble learning and geostatistical interpolation effectively captured both nonlinear relationships and spatial dependencies. The resulting high-resolution rainfall maps can support climate adaptation planning, irrigation scheduling, and sustainable management of water resources in data-scarce regions such as Pakistan.</p>