2 Library of salinization models in soilassessment package
2.1 Harmonization models
The soilassessment package contains three groups of harmonization models: 1) parametric models from the literature; 2) expressions for developing pedotransfer functions; and 3) expressions for developing own harmonization models (Omuto et al., 2023). Equation \(\ref{eq:1}\) is the general format for the harmonization models.
\[ y = f ( x, \theta ) \tag {1} \label{eq:1} \]
where y is the harmonized soil indicator, x is the indicator value to harmonize, and \(\theta\) is a set of coefficients (regression parameters) of the f harmonization model. The soilassessment contains more than 30 harmonization models (Omuto et al., 2023). Information about these models can be queried using the function SASmodels. The query needs two arguments: the type of soil indicator the models target (e.g., “ec” for electrical conductivity) and soil-water ratio for the extract solution in which the indicator was measured (e.g., “1:1”). The following example shows how the query is implemented using the function SASmodels:
library(soilassessment)
SASmodels("ec","1:1")
#> extract model texture
#> 1 1:1 FAO yes
#> 2 1:1 sonmez yes
#> 3 1:1 hogg yes
#> 4 1:1 landon no
#> 5 1:1 kargas no
#> 6 1:1 ozkan no
#> 7 1:1 USDA no
#> 8 1:1 hogg no
#> 9 1:1 zhang noThe query returns information about the extract solutions, names of the harmonization models, and whether the models require soil texture to harmonize the target soil indicator. Soil texture information is included in the query output because some models use it to model the harmonization.
2.1.1 Harmonization using parametric models from literature
Most salinization models in the soilassessment package are collections of parametric equations from the literature. They are grouped as either purely fixed-effect or mixed-effects models (Omuto et al., 2023). Fixed-effects models are average regression models of the form given in Equation \(\ref{eq:2}\).
\[ y_{i} = f (x_{i},\theta) + \epsilon_{i}\tag{2} \quad\text{for}\ {1 \leq i \leq n} \label{eq:2} \]
where yi is a vector of target variable, xi is a vector of soil indicator to harmonize, \(\epsilon_i\) is the random residual, n is the number of samples, and \(\theta\) is a set of coefficients (also known as fixed effects).
Fixed-effects models are regression models depicting the average relationship between x and y. They were developed with the assumption that a regression model is sufficient to represent the harmonization relationship between x and y without regards to any potential differences between natural soil groupings. Natural soil groupings such as texture or soil types affect soil properties and should not be ignored in soil property modelling. Mixed-effects models differ from the fixed-effects models in the way they treat grouping factors like texture, soil depth, climatic zones etc., which are known to influence salt characteristics in the soil (Seo et al., 2022; Shahid et al., 2018). Mixed-effects models not only recognize the grouping factors but also incorporate them in the modelling algorithms (Pinheiro and Bates, 2000). The format for mixed-effects models is given in Equation \(\ref{eq:3}\) and contains \(\phi\) model-fitting parameters comprising \(\theta\) fixed-effects and bj random-effects for each of the m grouping variables, and \(\epsilon_i\) random residuals.
\[ y_{ij} = f (x_{ij}, \phi_{j}) + g (\epsilon_{ij}) , \qquad\phi_{j} = \theta + b_{j} \qquad\text{ for}\ {1 \leq i \leq n} \qquad\text{and} \qquad\text{for}\ {1 \leq j \leq m} \qquad \tag{3} \label{eq:3} \]
where g is a function for the residual model.
Mixed-effects or fixed-effects harmonization models in the soilassessment package are implemented by specifying the function for the model, type of soil indicator to harmonize, model name, and the extract solution involved. Soil textural classes are optional for the models that require their use. There are five groups of functions in the soilassessment package for implementing harmonization models from the literature:
- The function ECconversion 1: for fixed-effects models in different soil textural classes
- The function ECconversion 2: for fixed-effects models which do not use soil texture
- The function ECconversion 4: for harmonizing EC and TSS or TDS
- The function ME_ECharm : for mixed-effects models to harmonize electrical conductivity (EC)
- The function ME_PHharm : for mixed-effects models to harmonize pH
The first group of functions are for models in the library which use fixed-effects equations. The equations were developed separately for each soil textural classes. The second group of functions are for fixed-effects models in the library where soil textural classes are not required. The fourth and the fifth groups use mixed-effects approach. The third group is for library models that harmonize TSS or TDS with EC of saturated soil paste extract (ECse). The following example shows how these harmonization models can be implemented in the soilassessment package. The example uses EC2 data (EC in 1:2 extract solution) which need to be harmonized into ECse.
# Get the data
EC2=c(1,1.4,0.8,0.3,0.3,1.3,2.8,0.7,0.6,3.5,1.5,2.5,0.9,0.5,0.8)
ECse=c(1.5,2.9,1.7,0.7,0.6,3,4.1,1.1,2.1,6,1.7,4.7,3,1.1,2)
texture=c("Cl","Cl","SiLo","Si","SiLo","Si","Lo","Lo","SaLo", "SaLo","LoSa",
"LoSa","LoSa","SiClLo","SiClLo")
soild = data.frame(EC2,ECse,texture)
# Harmonize using a fixed-effects model without use of soil texture
soild$ECse1=ECconversion2(soild$EC2,"USDA","1:2")
predAccuracy(soild$ECse,soild$ECse1)$Rsquared
#> [1] 0.8611294
# Harmonize using a mixed-effects model
soild$ECse2=ME_ECharm(soild$EC2,soild$texture,"linear","1:2")
predAccuracy(soild$ECse,soild$ECse2)$Rsquared
#> [1] 0.84345372.1.2 Harmonization using pedotransfer function
Pedotransfer function (PTF) is a mathematical model between soil properties in a soil database. PTFs are often developed between readily available and easily measurable soil properties and soil properties which are cumbersome and expensive to measure. They have been popularly used to estimate missing soil properties in a database (Van Looy et al., 2017). Readily available soil properties which are often used in PTFs to estimate target indicators of salts in the soil are textural proportions (percent sand, silt, and clay contents), organic carbon content, pH, cation exchange capacity (CEC), and exchangeable ions. Equation \(\ref{eq:4}\) is a general PTF format where the function f represents the PTF.
\[ y_{i} = f (x_{i},\theta) + \epsilon_{i}\tag{4} \quad\text{for}\ {1 \leq i \leq n} \label{eq:4} \] where yi is the vector of target variables, xi is a set of the predictor soil properties, \(\epsilon_i\) is the random residual, and n is the number of samples. The term \(\theta\) is optional for parametric models and represents a set of model-fitting parameters.
The soilassessment package allows for the development of PTFs using algorithms for linear regression (here denoted as linear), support vector machine (svm), random forest (randomforest), neural network (neuralnetwork), etc. These algorithms are also popularly used in soil science (Wadoux et al., 2020). PTFs are implemented in the soilassessment package using the function pedoTransfer. Arguments for this function include the target algorithm, data frame containing the target soil indicator to be harmonized, and a set of independent soil properties (predictors) to predict the target soil indicator. The PTFs need appropriate evaluation before using them. Table 2.1 is an example data for demonstrating the potential of PTF in gap-filling and harmonization. It contains gaps (NAs) in EC5 data (EC in 1:5 extract solution). The data is an extract from the global soil database for an area in southern Nigeria (Batjes et al., 2020).
The code below shows how PTFs are built in the soilassessment package using pedoTransfer function. Model building strategy where the data is split into two is recommended to provide data for model building and testing (Joseph, 2022). One part of the split data is used in developing the PTF and the other part is held-out for evaluating the PTF. Different PTF algorithms in the soilassessment package may be tested to arrive at the most suitable one for the dataset.
| EC5 | Sand | Silt | Clay | pH | OC | Ca | Mg | K | Na | CEC | BS |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.081 | 68.0 | 14.00 | 18.00 | 5.800 | 1.800 | 2.800 | 0.700 | 0.070 | 0.040 | 4.600 | 78.48 |
| 0.090 | 59.0 | 28.00 | 13.00 | 6.300 | 4.600 | 2.800 | 0.800 | 0.110 | 0.210 | 5.500 | 71.27 |
| 0.093 | 34.0 | 28.00 | 38.00 | 6.200 | 3.670 | 4.600 | 1.900 | 0.180 | 0.110 | 2.500 | 70.00 |
| 0.100 | 50.0 | 24.00 | 26.00 | 7.400 | 7.620 | 5.000 | 2.000 | 0.340 | 0.230 | 10.200 | 74.20 |
| NA | 72.0 | 8.00 | 20.00 | 6.100 | 2.400 | 1.200 | 0.300 | 0.420 | 0.280 | 2.000 | 50.00 |
| 0.102 | 54.0 | 12.00 | 34.00 | 5.600 | 0.600 | 5.600 | 2.300 | 0.100 | 0.070 | 8.600 | 70.60 |
| 0.130 | 59.0 | 32.00 | 9.00 | 6.800 | 6.100 | 3.200 | 1.300 | 0.180 | 0.170 | 6.400 | 75.78 |
| 0.150 | 52.0 | 24.00 | 24.00 | 7.300 | 3.990 | 4.600 | 1.700 | 0.300 | 0.100 | 7.400 | 90.50 |
| 0.160 | 58.0 | 36.00 | 6.00 | 7.000 | 4.100 | 2.500 | 0.900 | 0.360 | 0.350 | 4.700 | 87.06 |
| NA | 82.0 | 10.00 | 80.00 | 7.500 | 5.750 | 1.400 | 0.300 | 0.150 | 0.070 | 2.600 | 73.80 |
| 0.180 | 44.0 | 30.00 | 26.00 | 7.300 | 6.220 | 9.600 | 3.400 | 1.130 | 0.100 | 15.100 | 94.20 |
| 0.200 | 43.0 | 44.00 | 13.00 | 6.100 | 15.600 | 3.200 | 1.600 | 0.320 | 0.180 | 6.800 | 77.94 |
| NA | 48.7 | 21.20 | 30.10 | 5.775 | 0.605 | 7.195 | 3.095 | 1.235 | 0.890 | 12.810 | 33.17 |
| NA | 85.2 | 11.20 | 3.60 | 5.090 | 0.400 | 1.075 | 0.325 | 0.160 | 0.190 | 5.840 | 29.96 |
| NA | 84.6 | 8.10 | 7.30 | 5.355 | 0.565 | 1.750 | 1.465 | 0.585 | 0.360 | 7.270 | 69.55 |
| NA | 61.9 | 16.85 | 21.25 | 6.460 | 0.950 | 4.915 | 2.290 | 1.090 | 0.285 | 9.310 | 92.16 |
| NA | 58.1 | 19.20 | 22.70 | 6.425 | 0.945 | 3.860 | 1.675 | 1.675 | 0.455 | 12.275 | 62.44 |
| NA | 85.0 | 9.40 | 5.60 | 5.490 | 0.420 | 1.200 | 0.575 | 0.205 | 0.235 | 3.080 | 71.92 |
# select part of the data without NAs
psoil1=subset(psoil,!is.na(psoil$EC5))
# Randomly split data into two: 75% for training and 25% held-out for testing
set.seed(40)
bound = floor((nrow(psoil1)/3)*2)
randomsample = sample(nrow(psoil1))
Train = psoil1[randomsample, ][1:bound, ]
Test = psoil1[randomsample, ][(bound+1):nrow(psoil1[randomsample, ]), ]
ptf=pedoTransfer("randomforest",Train,EC5,Sand,Silt,Clay,pH,OC,Mg,K,Na,CEC,BS)
Test$EC5p=predict(ptf,newdata = Test)
predAccuracy(Test$EC5,Test$EC5p)$Rsquared
#> [1] 0.1311012PTF harmonization is used in the soilassessment package in two ways: 1) to gap-fill missing entries of a target variable, and 2) to develop a harmonization model between a standard salt indicator and other soil properties. It may be necessary to evaluate the output of the PTF harmonization using a graphical summary such as the histogram to check whether the harmonization model preserved the empirical distribution of the data (Figure 2.1).
Figure 2.1: Test Histogram comparison in PTF harmonization.
2.1.3 Developing own harmonization models using functions in soilassessment
More harmonization models are needed to complement the models in the soilassessment package. These additional models may be necessary to handle extract solutions other than those considered in the soilassessment package or to harmonize salt indicators which have not been adequately covered by the current models in the soilassessment package. New models may also be needed in areas where the soilassessment package models have inadequacies such as in poorly sampled parts of the world. All these cases may prompt the need to develop own harmonization models. The soilassessment package contains the function ECconversion 3 which can facilitate development of own harmonization model. ECconversion 3 allows model building using machine learning algorithms such as random forest, cart, svm, cubist, etc. The following code example illustrates how to develop own harmonization model using ECconversion 3(). In the example, data on apparent electrical conductivity (ECa) from an electromagnetic induction (EM38) and the corresponding EC of saturated soil paste extract (ECse) are used. The data was obtained from a study of soil salinity in an irrigation project in Taita Taveta County in Kenya in 1992 (Dooremolen, 1992). It contained 13 observations on soil surface (0-30 cm) in a 50m by 50m farm. The harmonization model is shown in Figure 2.2.
# Get the data and convert into a dataframe
x=c(0.78,1.071,0.048,0.967,0.91,0.429,0.282,0.743,0.499,0.197,
0.922,0.209,0.364) # ECa data
y=c(10.411,12.887,1.769,9.962,11.18,5.452,3.574,6.443,7.952,4.258,
9.774,7.311,9.984) # ECse data
xy=data.frame(cbind(x,y))
names(xy)=c("ECa", "ECse")
# Fit an exponential model
strt=c(A=0.1, B=0.8)
EC3.ml=nls(ECse~ECconversion3(ECa,A,B, "exponential"),
start=strt, data=xy)
predAccuracy(fitted(EC3.ml),xy$ECse)$Rsquare
#> [1] 0.6683013
plot(ECse~ECa,xy, xlab="Apparent EC (ECa, dS/m)",
ylab = "Measured EC (ECse, dS/m)")
lines(xy$ECa[order(xy$ECa)], predict(EC3.ml,
newdata=xy[order(xy$ECa),]),col="red")
Figure 2.2: Example of own harmonization model
2.2 Models for using remote sensing images and other spatial data
The soilassessment package also contains functions to accommodate use of factors influencing the distribution of salt in the soil. These factors are topography, climate, geology, hydrogeology, land use, proximity to salty water bodies such as seawater, groundwater, etc. The literature recommends considering these factors when assessing salt accumulation in the soil (Omuto et al., 2020). Spatial data representing these factors are remote sensing images and maps of soil types, geology, hydrogeology, elevation, climate, and land cover. They can be integrated into salt assessment through modelling. Besides modelling, there are also applications which have used reflectance from remote sensing images to directly assess the distribution of salt levels in the soil (Asfaw et al., 2018). Reflectance of remote sensing images contain integral information of land surface characteristics which can represent some of the factors influencing salt distribution in the soil. The reflectance can be directly correlated with soil surface characteristics or indirectly through image indices to assess occurrence and distribution of salts in the soil (Nguyen et al., 2020; Omuto et al., 2020). Many remote sensing image indices have been developed in the literature to support assessment of salts in the soil (Table 2.2). The soilassessment package contains the function imageIndices for deriving these image indices.
| Image Bands | Image Index | Estimation Model |
|---|---|---|
| swir1(b6),nir(b5),swir2(b7) | NSI | \((b6-b7)/(b6-b5)\) |
| green(b3),red(b4) | SI1 | \(\sqrt{(b4*b3)}\) |
| blue(2),red(4) | SI2 | \(\sqrt{(b4*b2)}\) |
| green(b3),red(b4) | SI3 | \(\sqrt{(b4^2*b3^2)}\) |
| swir1(b6),nir(b5) | SI4 | \((b5*b6-b6^2)/b5\) |
| blue(2),red(b4) | SI5 | \(b2/b4\) |
| red(b4),nir(b5),green(b3) | SI6 | \(b4*b5/b3\) |
| nir(b5),red(b4) | SAVI | \((b5-b4)/((b5+b4+0.5)*1.5\)) |
| green(b3),red(4),nir(b5) | VSSI | \(2*b3-5*(b4+b5)\) |
| red(b4),nir(b5) | NDSI | \((b4-b5)/(b4+b5)\) |
| red(b4),nir(b5) | NDVI | \((b5-b4)/(b5+b4)\) |
| blue(2),green(b3),red(4) | SR | \((b3-b4)/(b2+b4)\) |
| blue(2),green(b3),red(4),nir(5) | CRSI | \(\sqrt{(b5*b4-b3*b2)/(b5*b3+b3*b2)}\) |
| green(b3),red(b4),nir(b5) | BI | \(\sqrt{b3^2+b4^2+b5^2}\) |
The following code shows how to use the function imageIndices to derive remote sensing image indices for assessing the presence of salts in the soil. Inputs for the function are image bands of reflectance (blue, green, red, near infrared (nir), and shortwave infrared - (swir1 and swir2)) and the target image index for soil salt (Table 2.2). The input image bands can be supplied individually or as a stack of images. It is important to ensure that the input image bands are compatible in terms of temporal and spatial resolutions as well as the geographic extent and projection. Examples of outputs from the imageIndices function are shown in Figure 2.3.
library(raster)
#> Loading required package: sp
library(rasterVis)
#> Loading required package: lattice
green=raster("green.tif")
blue=raster("blue.tif")
red=raster("red.tif")
nir=raster("nir.tif")
swir1=raster("swir1.tif")
swir2=raster("swir2.tif")
SI1=imageIndices(blue,green,red,nir,swir1,swir2,"SI1")
NDSI=imageIndices(blue,green,red,nir,swir1,swir2,"NDSI")
SR=imageIndices(blue,green,red,nir,swir1,swir2,"SR")
names(SI1)=c("SI1")
names(SR)=c("SR")
names(NDSI)=c("NDSI")
levelplot(stack(SI1, SR, NDSI), scales =list(draw=TRUE),
col.regions = terrain.colors(255),
colorkey = list(space = "right"))
Figure 2.3: Examples of image indices for salt distribution in soil
#> Warning in sqrt(getValues(x)): NaNs producedRemote sensing indices can also be integrated with other factors to model the spatial distribution of salt in the soil. Several modelling algorithms are available in the literature for this purpose. The soilassessment package has a function for searching for suitable models and guide selection of the appropriate algorithm. It contains the function regmodelSuit which can search and report modelling statistics for top performing algorithms for the dataset. It tests potential algorithms on a dataset and reports their performance statistics such as mean error (ME), root mean square error (RMSE), correlation (r2) and Nash-Sutcliffe coefficient of efficiency (NSE). These statistics provide the guideline for selecting the most suitable algorithm. Models with the lowest ME and RMSE and the highest r2 and NSE should be preferred.
It is important to note the following points when preparing data for use with the function regmodelSuit (Omuto et al., 2020):
- The dataset containing the modelling variables need to be an R data frame
- It may be necessary to transform the target variable to enforce positivity where necessary
- In spatial assessment, predictors may be pixel values extracted from a raster grid of factors influencing salt distribution in the soil. These factors are popularly known as environmental covariates.
The following code example shows how the data is first prepared and organized. Environmental covariates used are remote sensing indices of salt in soil, climate variables (such as maximum and minimum temperature and rainfall), land cover, and topographic terrain parameters (such as slope, length of slope, plan curvature, channel network to base level - cnbl, etc.) (Omuto et al., 2020).
#Prepare and organize environmental covariates as predictors
dem=raster("dem.asc")
ls=raster("ls.asc")
slope=raster("slope.asc")
valley=raster("valley.asc")
cnbl=raster("cnbl.asc")
plan=raster("plan.asc")
maxtemp=raster("maxtemp.tif")
mintemp=raster("mintemp.tif")
lcover=raster("lcover.tif")
rain=raster("rain.tif")
dem=raster("dem.asc")
#stack the predictors
predictors=stack(BI,CRSI,NDSI,NDVI,NSI,SAVI,SI1,SI2,SI3,SI4,SI5,SI6,SR,
VSSI,dem,ls,slope,cnbl,plan,valley,maxtemp,mintemp,
rain,lcover)
names(predictors)=c("BI","CRSI","NDSI","NDVI","NSI","SAVI","SI1","SI2",
"SI3","SI4","SI5","SI6","SR","VSSI","dem","ls",
"slope","cnbl","plan","valley","maxtemp","mintemp",
"rain","lcover")Some times it is necessary to reduce the data dimensionality where there are more predictors than the number of cases of the target soil variable to be modeled. Principal Component Analysis (PCA) is one of the methods for reducing predictors to a few principal components (PCs) (Greenacre et al., 2022). The following code example shows how PCA can be used to reduce 14 image indices to a few PCs. A scree plot of the PCs (Figure 2.4) can aid selection of a suitable number of PCs which account for more than 95% of the information content in the original set of image indices.
library(factoextra)
#> Loading required package: ggplot2
#> Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
predictors2=(predictors[[1:14]])
pred.val = getValues(predictors2)
pr.id = which(!is.na(pred.val))
pca = princomp(pred.val, cor=T)
fviz_eig(pca) # scree plot
Figure 2.4: Scree plot for choosing a suitable number of PCs
ncomp = 5 # Selected number of PCs
pr.pca <- predictors2[[1:ncomp]]
for(i in 1:ncomp) { pr.pca[[i]][pr.id] <- pca$scores[,i] }
names(pr.pca)=c("PC1","PC2","PC3","PC4","PC5")
predicters=stack(pr.pca$PC1,pr.pca$PC2,pr.pca$PC3,pr.pca$PC4,pr.pca$PC5,
dem,ls,slope,cnbl,plan,valley,maxtemp,mintemp,rain,lcover)In this documentation, a sample soil data was used with the environmental covariates to show how to build the spatial relationship for mapping indicators of salts in the soil. The sample soil data contains measured soil properties from 78 sampling locations in southern Kenya. Each sampling location had more than four samples at varied depth intervals in a soil profile. A spatial distribution of the sample soil data is shown in Figure 2.5.
Figure 2.5: Example of locations of soil pH data from southern Kenya
The pixel values of the covariates are extracted and appended onto the soil data as shown in the following example.
#Import the point shapefile and extract predictor pixel values
soil1=extract(predictors, soil)
soil2 =cbind(soil,soil1)
#organize the target and predictor into dataset for modelling
soilp=soil2[,c("EC5","TEXCLASS","ECse","pH","ESP","BI","CRSI","NDSI",
"NDVI","NSI","SAVI","SI1","SI2","SI3","SI4","SI5","SI6",
"SR","VSSI","dem","ls","slope","cnbl","plan","valley",
"maxtemp","mintemp","rain","lcover")]
soilp=data.frame(soilp)When using regmodelSuit function, the first argument is the dataset containing modelling variables, the second argument is the target variable to model, and the subsequent variables are the predictors. The names of the variables in the dataset must correspond to the names of the arguments supplied to the regmodelSuit function. The function may take time to complete the process especially where a large dataset is involved. The output is a table showing evaluation statistics for top suitable models for the dataset (Table 2.3).
soilec =soilp[, colnames(soilp)[c(3, 6:29)]]
regmodelSuit(soilec,ECse,BI,CRSI,NDSI,NDVI,NSI,SAVI,SI1,SI2,SI3,SI4,SI5,SI6,
SR,VSSI,dem,ls,slope,cnbl,plan,valley, maxtemp,mintemp,rain,
lcover)# NB: This part can take time| Algorithm | ME | RMSE | R2 | NSE |
|---|---|---|---|---|
| Linear | 0.4980970 | 0.9221965 | 0.4131760 | 0.0545453 |
| RandomForest | 0.3200139 | 0.8874233 | 0.6415919 | 0.4811995 |
| SVM | 0.4051999 | 1.0922583 | 0.3448001 | -5.9226998 |
| BayesianGLM | 0.4860159 | 0.9464917 | 0.4202483 | -0.3776671 |
| BaggedCART | 0.4032268 | 0.8978902 | 0.4444890 | 0.0577982 |
| Cubist | 0.3010037 | 0.8355491 | 0.5804149 | 0.5111329 |
| CART | 0.4544668 | 1.0954523 | 0.2593884 | -0.0213372 |
| Ranger | 0.3285088 | 0.9232545 | 0.5128640 | 0.4577277 |
| QuantRandForest | 0.3236330 | 1.0494950 | 0.4848026 | -0.6961260 |
| QuantNeuralNT | 0.3827372 | 1.1198034 | 0.4613748 | -4.2935941 |
Analysis of the output of the function regmodelSuit in Table 2.3 shows that the randomForest, cubist, and ranger were likely the best algorithms for spatial modelling of salt levels in the sample data. They had low ME and RMSE and high r2 and NSE. Any of these three algorithms can be selected for modelling the spatial distribution of the target variable as shown in the following example code. Example output prediction is shown in Figure 2.6.
library(caret)
#summary(soilp$EC5) #Check for missing values
#soilp =subset(soilp, !is.na(soilp$EC5)) #Remove missing data
# Mixed-effects polynomial model
soilp$ECse_pl=ME_ECharm(soilp$EC5,soilp$TEXCLASS, "polynomial","1:5")
# Mixed-effects linear model
soilp$ECse_lin=ME_ECharm(soilp$EC5,soilp$TEXCLASS, "linear","1:5")
trCl=trainControl( method = "cv",number=5,returnResamp = "all",
savePredictions = TRUE, search = "random",
verboseIter = FALSE)
ec.m=train(sqrt(ECse_pl)~(BI+CRSI+NDSI+NSI+SAVI+SI1+SI2+SI3+SI4+SI5+SI6+SR+
VSSI+ dem+ls+slope+cnbl+plan+valley+maxtemp+
mintemp+rain+lcover),soilp, method = "rf",
trControl=trCl) # Repeat for linear harmonization model
ECse=(predict(predictors, ec.m))^2
Figure 2.6: Spatial distribution of EC, pH, and ESP
2.3 Models for classification of levels of salt content in soil
Classification of soil attempts to put different soils into groups of similar characteristics. It provides the general guidelines for estimating characteristics of any soil in the group since all soils in each class are assumed to have the same behavior or characteristics. Most salt classification schemes in the literature use EC, pH, and ESP or SAR to group soils into different classes. Popularly used classes are saline, saline-sodic, and sodic, because they represent majority of problematic and widespread salt-affected soils worldwide (Abrol et al., 1988; Szabolcs, 1989). Classification schemes identify these classes using certain thresholds on pH, EC, and ESP or SAR. The most popular schemes are those proposed by FAO and USDA (Omuto et al., 2020). Figure 2.7 is an example of a classification scheme by FAO (FAO, 2006).
Figure 2.7: FAO salt classification scheme
The soilassessment package contains the function saltSeverity for implementing popularly used salt classification schemes. Its implementation requires a list of the soil indicators for the classification and the preferred classification scheme. The following example illustrates how to implement saltSeverity function for estimating salt levels in the soil using FAO classification scheme in Figure @ref(fig:SAStree0.
#> terra 1.8.5
Saltlevel=saltSeverity(values(ec1$ECse),values(ph1$PH),values(esp1$ESP),
"FAO")
esp1$salt=Saltlevel
esp1$salta=classCode(values(esp1$salt),"saltseverity")
#
par(mfrow=c(1,2))
pal=c("magenta","brown","yellow","white","gray","cyan","green"
,"peachpuff1", "orange","red")
plot(esp1$salta2,cex.main=1, cex=2,col=pal, main="(a) Linear model")
plot(esp1$salta, cex.main=1,cex=2, col=pal2,
main="(b) Linear mixed-effects model",legend=FALSE)
Figure 2.8: Example of distribution of salt levels
The distribution of salt classes in Figure 2.8 demonstrate the impact of not selecting appropriate harmonization models. It shows results of a linear harmonization model that did not consider soil texture and that of a mixed-effects harmonization model which considered soil texture. Comparison of the performance of the harmonization models revealed that the linear model did not only poorly harmonize high EC but also had low overall performance compared to the linear mixed-effects model. Consequently, the linear model produced more areas with high salinity than the case with the linear mixed-effects (Figure 2.8).
