High soil variability necessitates a large number of samples, which poses disadvantages in terms of labor, time, and economic and environmental impacts. Utilizing spectroradiometers and geostatistical methods can lead to significant savings in chemical inputs and time. In this study, sixty surface soil samples from basaltic parent material areas were analyzed in the laboratory for their physical (clay, silt, sand), chemical (pH, EC, exchangeable cations: Ca, Mg, Na, K, CEC, percent CaCO3) and biological (soil organic matter; OM100mμ, OM2mm) properties. Spectral and geostatistical methods were employed to estimate and map these properties. Spectral reflectance were obtained within the 350 to 2500 nm wavelength range. Modeling the relationships between laboratory measurements and spectral readings were performed using Partial Least Squares Regression (PLSR). Additionally, geostatistical techniques such as Inverse Distance Weighting (IDW), Ordinary Kriging (OK), and Cokriging (COK) were utilized to generate maps illustrating the spatial distribution of soil parameters. The accuracy of the predictions were evaluated using RMSE (Root Mean Square of Estimation) parameter. The predictive success of prediction techniques varied depending on the specific soil property under investigation. The VNIRS-PLSR method achieved the highest accuracy and the lowest RMSE values for parameters such as organic matter, sand, clay contents, cation exchange capacity (CEC), and electrical conductivity (EC). Conversely, geostatistical methods yielded the lowest RMSE results for parameters such as lime (CaCO3), pH, silt, exchangeable Ca, exchangeable K, exchangeable Na, and exchangeable Mg. The application of the COK technique using a secondary variable resulted in a 1 % to 19 % increase in prediction success compared to OK and IDW techniques. Overall, each estimation technique has its own advantages and disadvantages, which should be taken into consideration in the selection of the technique for prediction of soil variables.