An accurate estimation of soil organic matter (SOM) content for spatial non-point prediction is an important driving force for the agricultural carbon cycle and sustainable productivity. This study proposed a hybrid geostatistical method of extreme learning machine-ordinary kriging (ELMOK), to predict the spatial variability of the SOM content. To assess the feasibility of ELMOK, a case study was conducted in a regional scale study area in Shaanxi Province, China. A total of 472 topsoil (0-20 cm) samples were collected. A total of 14 auxiliary variables (predictors) were obtained from remote sensing data and environmental factors. The proposed method was compared with the ability of traditional geostatistical methods such as simple kriging (SK) and ordinary kriging (OK), in addition to hybrid geostatistical methods such as regression-ordinary kriging (ROK) and artificial neural network-ordinary kriging (ANNOK). The results showed that the extreme learning machines (ELM) model used principal components (PCs) as input variables, and performed better than both multiple linear regression (MLR) and artificial neural network (ANN) models. Compared with geostatistical and hybrid geostatistical prediction methods of SOM spatial distribution, the ELMOK model had the highest coefficient of determination (R 2 = 0.671) and ratio of performance to deviation (RPD = 2.05), as well as the lowest root mean square error (RMSE = 1.402 g kg −1 ). In conclusion, the application of remote sensing imagery and environmental factors has a deeper driven significance of a non-linear and multi-dimensional hierarchy relationship for explaining the spatial variability of SOM, tracing local carbon sink and high quality SOM maps. More importantly, it is possibly concluded that the sustainable monitoring of SOM is a significant process through the pixel-based revisit sampling, an analysis of the mapping results of SOM, and methodological integration, which is the primary step in spatial variations and time series. The proposed ELMOK methodology is a promising and effective approach which can play a vital role in predicting the spatial variability of SOM at a regional scale.