Assessing the spatial variability of soil organic carbon (SOC) is crucial for SOC monitoring and comparing management options. Topsoil (0–5 cm) SOC concentrations were surveyed in a coffee agroforestry watershed (0.9 km2) on Andisols in Costa Rica with uniform farm management. We encountered high values and large spatial variations of SOC, from 48.1 to 172 g kg- 1 in the dry combustion set (SOCref; n = 72) used for calibrating the visible-near-infrared reflectance spectroscopy (VNIRS) samples (SOCVNIRS; 350–2500 nm; n = 520). VNIRS using partial least squares regression was effective in predicting SOC (R2 = 0.85; a root mean square error (RMSE) = 12.3 g kg- 1) and proved an effective proxy measurement. We assessed several topographic, vegetation and andic soil property variables, of which only the latter (metal–humus complexes and allophanes) displayed strong correlations with SOCref concentrations. We compared Random Forest and three geostatistical approaches for the interpolation of SOC in unsampled locations. Ordinary kriging with SOCref yielded an RMSE of 28.0 g kg- 1. Random Forest was successful in incorporating many weakly and non-linearly correlated covariates with SOC (RMSE = 14.7 g kg- 1), provided Alp (the sodium pyrophosphate extractable aluminum), the best predictor of SOC (r = 0.85) but also the most costly variable to acquire. Co-kriging with Alp also showed high reduction in RMSE (16.0 g kg- 1). Co-kriging with SOCVNIRS only showed marginal reduction in RMSE to 24.2 g kg- 1 due to the presence of a high nugget effect. Local variability of SOC in this volcanic agroforestry watershed was dominated by andic properties whereas topographic or vegetation variables had very little impact. Estimation of SOC variability is recommended using inexpensive proxy measurements like VNIRS (RMSE = 12.3 g kg- 1) rather than spatial interpolation techniques.