The study compares two approaches to geostatistical analysis (kriging, cokriging) and ANN for evaluating and forecasting the spatial distribution of chemical elements in the topsoil by the example of a small area in natural geological landscape. The size of the selected test plot (1.0 x 1.0 m) allows assuming its homogeneity, thus providing suitable conditions for the application of geostatistical analysis. However, the authors believed that even under such conditions the trained ANN will give prediction models comparable in accuracy with the geostatistical methods. The test plot was split into 100 cells. Soil cores were sampled to a depth 0.05 m by stainless steel cylindrical sampling device with a diameter of 0.05 m. For the analysis of soil specimens, the X-ray fluorescence spectrometry was chosen. First, the network model for constructing the distribution of each element contained in the specimen was selected. A multilayer perceptron with Levenberg-Marquardt training algorithm was chosen for this study. Next, concentration distributions for each element using ordinary kriging and ordinary cokriging based on data of X-ray fluorescence analysis by Esri ArcGIS geostatistical software were built. Second, in order to compare the predictability of two methods, the sample was randomly divided using the “create subset” function of Geostatistical Analyst in Esri ArcGIS in two sub-sets, i.e., the test and the training in the ratio of 30 : 70. The training sub-set (70 specimens) was used as the learning set. Then the concentration value of each element was returned by ANN, kriging, and cokriging by the test sub-set (30 specimens). Three validation indices were used to evaluate the performance of different interpolation methods: the mean absolute error (MAE), root mean square error (RMSE), and relative root mean square error (RRMSE).