Predicting and Mapping of Soil Organic Carbon Using Machine Learning Algorithms in Northern Iran

In plain words

The paper asks how much carbon is stored in the topsoil of a region in northern Iran, and it tries to draw a map of that amount from places nobody has sampled. Soil carbon matters for soil health and for how much carbon the land holds. Measuring it means digging and lab work, so the authors combine about 1900 soil samples with many cheap map layers: satellite images, terrain shape, rainfall and temperature, and soil and land-use maps.

The map layers are too many to use all at once, so a search procedure, inspired by natural selection, first picks the most useful ones. Then six standard prediction methods are compared. A deep neural network does best, and it comes with a sense of how uncertain each predicted value is. The result is a baseline map of the region's soil carbon.

TL;DR

• Predicts topsoil organic carbon (SOC, %) in Mazandaran province, Iran, from 1879 samples and 105 covariates §2.2, §2.3.
• A genetic algorithm (GA) selects 35 covariates §2.4, §3.2. Six regressors are compared under 10-fold cross-validation §2.6, tab. 4.
• The deep neural network (DNN) has the lowest error and the highest R² and CCC tab. 4. The paper recommends it for mapping with uncertainty §3.3.
• Author-reported results on one study area. No official code is linked.

Problem

Soil organic carbon varies strongly and is far from normally distributed (skewness 2.33, kurtosis 8.2) tab. 2, §3.1. Mapping it needs many correlated auxiliary variables, and the authors argue that classical methods struggle with that complexity §1. Feature selection is needed to remove redundant covariates §2.4. The authors state that no best-performing algorithm is known worldwide and that results depend on the local area §1.

Method

Data §2.2. 1879 composite topsoil samples, 0–20 cm deep. Each composite comes from at least 10 cores within a 20 m radius. After air-drying and sieving to 2 mm, SOC was measured by wet oxidation (Walkley–Black, cited in the paper as [65]). About 1055 samples come from five M.Sc. projects, mostly uncultivated areas. The rest come from a 2000 m grid survey of cultivated areas §2.2, fig. 1.

Covariates §2.3, tab. A1. 105 candidates in four groups:

• 60 remote-sensing variables: Landsat 8 median of 8 scenes from 2012–2016, plus MODIS red and near-infrared at 250 m and day and night land surface temperature at 1 km §2.3;
• 30 terrain attributes from a 30 m SRTM DEM, computed in SAGA GIS §2.3;
• 10 WorldClim v2 climate variables at 1 km §2.3;
• 5 categorical maps (soil and land use, among others) §2.3.

All layers were resampled to 30 m by nearest-neighbour or bilinear resampling. The paper does not say which method was used for which layer §2.3. The full list of 105 with their formulas is in tab. A1.

Feature selection with a GA §2.4, §3.2. Each individual is a binary mask over the 105 covariates. Fitness is the RMSE of a random forest, scored with 10-fold cross-validation. Each generation keeps the two fittest individuals as parents, then applies crossover and mutation. Population size is 50, crossover rate 0.6 and mutation rate 0.001 §2.4. The generation with the best external RMSE gave 35 predictors: 13 terrain, 18 remote-sensing, 2 climate and 2 categorical §3.2, tab. 3.

# Genetic-algorithm covariate selection  [§2.4]
# Input:  covariates (1879, 105); soc (1879,)
# Output: binary mask over 105 covariates

def select_covariates(covariates, soc):
    # Random starting population of binary masks  [§2.4]
    population = random_binary_masks(count=50, n_features=105)
    for generation in range(n_generations):  # count not stated, see §2.4 vs §3.2
        # Fitness: RMSE of a random forest, 10-fold CV  [§2.4]
        scores = [cv_rmse_random_forest(covariates, soc, mask)
                  for mask in population]
        # Two fittest (lowest RMSE) are the parents  [§2.4]
        parents = two_lowest_rmse(population, scores)
        # Crossover rate 0.6, mutation rate 0.001  [§2.4]
        population = crossover_and_mutate(parents, p_cross=0.6, p_mut=0.001)
    # Keep the mask from the generation with best external RMSE  [§3.2]
    return best_external_mask(population)

Regressors §2.5, tab. 1. Six models are fitted on the 35 selected covariates §1:

• SVM with an RBF kernel §2.5.1, tab. 1;
• ANN, one hidden layer, sigmoid activation, in the nnet package §2.5.2;
• Cubist regression tree, with committees and neighbours §2.5.3;
• random forest, tuned for Mtry and Ntree §2.5.4;
• XGBoost with gbtree boosters §2.5.5;
• DNN with the rectifier activation, trained in H2O §2.5.6.

DNN §2.5.6, eq. 1–4. Hidden layers use the form z = σ(W·z_prev + b). The output layer is linear, Y = W^(L+1) z^L + b^L eq. 2. The loss is the mean squared error over the training set eq. 4. The paper states that a GA is used to minimise the loss §2.5.6. The final configuration in the abstract has 7 hidden layers and 50 neurons per layer. Table 1 gives the search ranges instead.

# Feed-forward DNN regressor  [§2.5.6] [eq. 1–3]
# Input:  x (B, 35) selected covariates
# Output: predicted soc (B,)

def dnn_forward(x, weights, biases):
    # Hidden layers, rectifier activation  [§2.5.6] [eq. 1]
    h = x
    for layer in range(n_hidden):  # 7 in the final model  [abstract]
        h = relu(h @ weights[layer] + biases[layer])  # (B, 50)
    # Linear output layer  [eq. 2]
    return h @ weights[n_hidden] + biases[n_hidden]  # (B,)

Validation §2.6, eq. 5–8. The 1879 samples are split into 10 sets. Each set serves once as test data. Each model is refitted for every fold, so each point is predicted at least once. Four metrics are computed per fold and averaged: R², Lin's CCC, MAE and RMSE eq. 5–8. The best model has the lowest MAE and RMSE and the highest R² and CCC §2.6.

Uncertainty §2.7. For each pixel, the mean and standard deviation of the predictions across the 10 realisations give a normal-approximation 90% interval (±1.64 SD). The paper also says the 5th and 95th percentiles were mapped §2.7. Three maps are produced for the best model: mean, lower and upper limit §2.7, fig. 7.

Data transforms §3.1, tab. 2. SOC is anchored at 1.00 and then log-transformed. Skewness drops from 2.33 to 0.59 and kurtosis from 8.2 to 1.50. Predictions are back-transformed to the original scale §3.1.

Hyperparameters

Results

Author-reported, 10-fold cross-validation over 1879 samples tab. 4. SOC in %. Means ± SD across folds.

Headline numbers, from the abstract §abstract:

• The DNN has MAE 0.59%, RMSE 0.75%, R² 0.65 and CCC 0.83.
• Precipitation explains 14.9% of SOC variability, followed by NDVI (12.5%), MODIS day temperature (10.6%), MrVBF (8.7%) and land use (8.2%) fig. 5.

Prediction intervals, 90% tab. 5: the fraction of the 1879 points inside the interval is SVM 79.30%, ANN 81.11%, Cubist 84.09%, RF 82.97%, XGBoost 84.46%, DNN 87.81%. The DNN is the best calibrated, but its interval coverage is still below 90%.

Summary statistics tab. 2: min 0.02%, max 11.48%, mean 2.19%, SD 1.27%, CV 58.23%.

Group means §3.5–3.6, tab. 6, fig. 8–9: the highest SOC is under dense forest (3.77%) and the lowest is on seashore land (0.75%). Ultisols (4.04%) and Mollisols (3.20%) have the highest mean SOC by soil order. These figures describe the study area only.

Paper vs code

No official code. The paper gives no repository link and no project page § title, abstract, footnotes. Searches: GitHub repository search for "soil organic carbon machine learning Emadi" and "soc mazandaran" returned zero results. A GitHub user search for "taghizadeh-mehrjardi" returned zero results.

Software named in the paper: R with the caret package §2.4, ref. 72, Cubist §2.5.3, ref. 84, nnet §2.5.2 and H2O §2.5.6. Versions and scripts are not stated.

Pitfalls

• Units in the arXiv abstract. The arXiv abstract says "59 percent" MAE and "75 percent" RMSE. The PDF abstract and tab. 4 give 0.59% and 0.75% of SOC. Use the PDF values, which are in SOC percentage units.
• GA iteration count conflicts. §2.4 says 100 iterations. §3.2 says the GA ran 1000 times and shows 1000 generations in fig. 4. Termination rule and count are unclear.
• Final DNN configuration is only in the abstract. The abstract gives 7 hidden layers and 50 neurons. tab. 1 gives searched ranges (2–10 layers, 15–200 neurons). fig. 3 labels 15–100 neurons. The search procedure and selection criterion are not stated.
• Training details are missing. Optimiser, learning-rate schedule, epochs, batch size, dropout value, weight initialisation and random seed are not stated for the final DNN. The paper says the GA minimises the loss §2.5.6, but H2O trains with its own optimiser. Whether a GA was used for the DNN weights is therefore unclear.
• Evaluation may be optimistic. The GA picks the generation with the best external RMSE §3.2 and reports 10-fold CV on the same samples. The paper describes no separate held-out test set.
• Fold construction is not stated. It is not said whether folds are random or spatially blocked. Spatial autocorrelation could make random-fold CV optimistic.
• Log-transform detail. "Anchored at 1.00" is not defined. The back-transform bias correction, if any, is not stated §3.1.
• Uncertainty method is internally inconsistent. §2.7 gives both a ±1.64 SD normal interval and 5th/95th percentiles. It is not stated whether the 10 realisations are the same models as the CV folds.
• Equation typos. Eq. 4 is typeset as |·|² rather than (·)², and eq. 5 uses p in the numerator where P_i is meant. Both read as ordinary MSE and R² once corrected.
• Citation errors. H2O is cited as [92], which is a saline-soil paper. The correct reference is [79] (Candel et al., Deep Learning with H2O). The nnet package is also cited as [79].
• Table 1 label. The ANN row labelled "decay" is defined as a learning rate. The ANN "size" range is 1–10 in tab. 1 but 2–10 in the text.
• Sample composition. About 56% of samples come from M.Sc. projects, mostly in uncultivated areas, and the rest from a grid survey of cultivated areas §2.2. Sampling bias is not discussed.
• Predictor resampling. Coarse layers (250 m MODIS, 1 km WorldClim) were resampled to 30 m. Which method applies to which layer is not stated §2.3.
• Variable importance method. fig. 5 uses a sensitivity analysis from ref. [45], per §3.2. The exact procedure is not given in the paper.
• Single study area. The paper claims the DNN is a promising algorithm for province-scale mapping [abstract], but it tests only Mazandaran province. The authors themselves recommend testing it elsewhere in Iran §4.

Lineage

• Builds on: the DNN formulation of Wang et al. (ref. 94), which the paper cites for eq. 1–3 §2.5.6. The GA-for-covariate-selection approach follows Welikala et al. (ref. 52) §2.4.
• Compared against: the five other regressors in tab. 1 and tab. 4, all fitted in this study. No external baseline is reported.
• Superseded by: not known from content in this wiki.