Visible-near infrared reflectance spectroscopy (VNIR) is a promising method for rapid spatial assessment of a range of soil properties. Raw hyperspectral VNIR data may be combined with first and second-derivatives to estimate information about soil constituents. In this study, we evaluated whether VNIR or geostatistical methods, or their combined use, could provide more efficient approaches to spatial soil assessment. We tested the prediction accuracy of VNIR reflectance spectroscopy and kriging methods using soil samples from a grid (800x400m) sampling design of 32 ha, established in Northern Turkey. Soil variables considered were CaCO3, organic matter, clay, sand, silt (%), pH, electrical conductivity (EC, dsm-1) and cation exchange capacity (CEC), Ca, Mg, Na and K (cmolkg-1) content. In order to compare the two methods, systematically selected 13, 25, 50 and 100 % of the grid data (n=512) were used for calibration and the remaining 87, 75 and 50 %, respectively, were used for validation. Partial Least Square regression was used for calibrating soil properties to both raw and first derivative reflectance spectra, while ordinary and co-kriging methods were used for spatial interpolation. In order to evaluate accuracy of the predicted results, R2 values were first obtained by plotting predicted versus observed values, and then R2 and RMSE of predictions were compared for each method. Results show that the raw spectra and first derivatives provided similar results. The VNIR method provided better prediction results than ordinary kriging for soil organic matter and slightly better for clay and sand from the use of either raw spectra or first derivative data, producing higher R2 (range of 0.55-0.73, 0.71-0.81, 0.57-0.78, respectively) and lower RMSE (range of 0.27-0.47, 3.74-5.30, 4.69-6.92, respectively). EC, pH, Na, K and silt content were poorly predicted using both methods, either due to the low range or the random nature of spatial variability. Overall, the VNIR prediction accuracy was not much impacted by sample size as much as the one for ordinary kriging, whose results improved with increased sample size of the calibration model. We also performed co-kriging-based prediction after principal component (PC) analyses of the VNIR data using 25 % of the data to predict the remaining samples for soil variables of clay, sand, CEC, soil organic matter, and Ca, which have significant correlations with PC2, PC3 and PC4 (r > 0.60). According to cokriging analysis cross validation results, clay was best predicted, ( r2 = 0.74 RMSE = 4.92) and others moderately; sand (R2= 0.62, RMSE=6,82), CEC (R2 = 0.65, RMSE= 2.12), Ca ( R2=0.66, RMSE= 1.73), SOM (R2= 0.43, RMSE=0.54) and less accurately Mg (R2=0.30, RMSE= 0.55) and silt (R2= 0.16, RMSE= 5.40). The combined use of VNIR spectroscopy and geostatistics provides considerable opportunities for spatial assessment of soil properties at reduced cost compared to methods that are solely based on conventional sample analysis. In this approach, a small fraction of soil samples is analyzed in the lab, and the combination of VNIR and geostatistics is applied to estimate soil properties at other sample locations.