Spectral Method for Liming Recommendation in Oxisol Based on the Prediction of Chemical Characteristics Using Interval Partial Least Squares Regression

: Thousands of chemical analyses are carried out annually with the aim of recommending soil correction; however, these analyses are expensive, destructive, time-consuming, and can be harmful to the environment. As an alternative to conventional analysis methods, diffuse reﬂectance spectroscopy has been proposed as an option for evaluating the chemical characteristics of soil. The selection of variables has also emerged as an alternative to improve the performance of PLSR (partial least squares regression), as it decreases the root mean square error (RMSE) and increases the accuracy of the models. However, few studies have used a previous selection of variables for the construction of PLSR models to estimate the chemical characteristics of soil. In this context, the hypothesis in this study was that it is possible to calculate the liming recommendation in Oxisol based on the chemical characteristics estimated by PLSR, with a previous selection of variables using iPLS (Interval PLS). The objective was to calculate the need for liming based on chemical characteristics estimated via iPLS selection and PLSR modeling of speciﬁc wavelengths of soil reﬂectance. The experimental area was treated with different application rates of limestone, with and without incorporation, and phosphogypsum was applied in additional treatments. Soil assessments were carried out 5, 12, 24, and 36 months after the application of the treatments, using six layers: 0.00–0.05, 0.05–0.10, 0.10–0.20, 0.20–0.30, 0.30–0.40 and 0.40–0.60 m. Samples were subjected to conventional laboratory analyses, and spectral readings (400–2500 nm) were obtained with a spectroradiometer. The spectral curves were subjected to the iPLS variable selection method to generate PLSR models of the chemical characteristics used to calculate the liming recommendation. The chemical characteristics of the soil, such as Ca 2+ , sum of bases (SB), effective cation exchange capacity (CTCe), cation exchange capacity (CTC), and base saturation (BS), could be estimated, with values of R 2 ranging from 0.83 to 0.92 in the calibration and validation steps, and from 0.84 to 0.90 for the prediction step (in the fourth assessment). The liming recommendation calculated based on the chemical characteristics predicted from the PLSR models showed a strong correlation (r > 0.86) with the liming recommendation calculated by conventional laboratory techniques. The fourth soil assessment yielded the best correlation coefﬁcient (r = 0.95).


Introduction
A large proportion of the world's agriculture in tropical and subtropical environments is carried out on oxisols. Most Brazilian oxisols impose limitations on the establishment and development of the main crops due to the acidity of these soils, the presence of toxic elements (such as Al 3+ and Mn 2+ ), and low levels of basic cations (such as Ca 2+ and Mg 2+ ), which restrict the development of the plant root system and reduce their ability to absorb water and nutrients, especially in regions or periods of low water availability [1][2][3][4][5]. The adverse effects of soil acidity on plants occur both in the surface layer (0.0-0.2 m) and in the subsurface (depth > 0.2 m) of acidic soils. Thus, to promote greater efficiency of water and nutrient absorption by plants, and hence to enlarge yields, the surface acidity must be reduced and the chemical environment in the subsurface improved.
The application of limestone (liming) aims to correct acidity (increasing the pH), which reduces the toxic effects of high concentrations of Al 3+ and Mn 2+ and increases the levels of Ca 2+ and Mg 2+ in the soil. Liming also increases the effective CTC and decreases the leaching of bases, and is, therefore, considered one of the practices that most contributes to the increase of fertilizer efficiency, and as one of the pillars for obtaining higher yield and increased agricultural profitability and sustainability [2,[6][7][8]. Since liming mainly improves the characteristics of the soil surface, especially in a no-till system, the application of phosphogypsum is an alternative to surface liming, as this improves the root environment in the subsurface of the soil and, under certain conditions, makes limestone incorporation unnecessary [9,10].
In 2020, Brazilian agribusiness used about 45 million tons of limestone, representing an increase in consumption of almost 60% over the previous 20 years and 45% over the previous 10 years [11]. Despite this increase in limestone consumption in Brazil, it is estimated that the ideal consumption would be around 70 million tons per year [11]. This increase in limestone consumption is directly related to soil analysis. The characteristics used to estimate liming recommendations are obtained through conventional laboratory analyses; although these are reliable techniques, they are expensive, destructive, time-consuming, and require reagents which can cause harm to the environment. Furthermore, the sample density is limited by the costs of laboratory tests and operational logistics [12][13][14][15].
As an alternative to conventional soil analysis methods, diffuse reflectance spectroscopy has been used to quickly, non-destructively, and cost-effectively assess soil characteristics [16,17]. The use of this technique to determine the chemical, physical and mineralogical properties of soil has shown favorable results in the regions of the visible (Vis) (350-750 nm), near-infrared (NIR) (750-1100 nm), and shortwave infrared (SWIR) (1100-2500 nm) ranges. The availability of such a large set of spectral information allows researchers to estimate multiple soil characteristics from a single reflectance spectrum [16,[18][19][20][21].
Several studies have demonstrated the potential of diffuse reflectance spectroscopy to estimate exchangeable soil ions, such as Ca 2+ , Mg 2+ , K + , P, and SO 4 2− , and pH [22][23][24][25][26][27][28]. This is made possible by the interaction between electromagnetic radiation and soil particles, which causes chemical bonds to vibrate at specific frequencies. This specificity allows for the evaluation of different soil characteristics once the reflectance is known [29,30].
Multivariate statistical techniques can be used to build prediction models for specific soil characteristics based on soil reflectance. However, these techniques need to be robust and accurate, given the large volumes of data. Among the different statistical methods, partial least squares regression (PLSR) is one of the most used [27,[31][32][33][34].
Initially, it was supposed that the inclusion of more spectral regions would improve the accuracy of PLSR models, but several studies have found that the use of the full spectrum can cause overfitting. The selection of variables technique has emerged as an alternative to minimize this problem, as it improves the performance of regression, decreases the root mean square error (RMSE), and increases the accuracy of the models [35][36][37][38]. In this context, an interactive algorithm called Interval PLS (iPLS) was developed to generate local PLS models in subintervals of the spectral curve of interest. iPLS selects variables with relevant Remote Sens. 2022, 14, 1972 3 of 22 information in different spectral subdivisions, thus focusing on critical spectral regions and removing interference from other regions [39].
The iPLS algorithm selects a relevant subset of intervals (a subset of wavelengths) for each of the soil characteristics individually, which results in superior models compared to the use of the entire Vis-NIR spectrum, and hence decreases the complexity of the models by removing redundant information from the spectra [40]. When iPLS has split the spectra into different ranges, a PLSR validation model is applied to each of them. The interval that gives the lowest prediction error is modeled again by adding the remaining intervals, one by one, until the error no longer decreases [41].
Although many studies have reported the use of the Vis-NIR-SWIR range to predict soil characteristics, studies that have used a previous selection of variables for constructing PLSR models to estimate the necessary chemical characteristics are still rare in terms of calculating liming recommendations. These could make it possible in the future to develop cheaper sensors (or even cameras using specific wavelengths) which would not need to operate in the full Vis-NIR-SWIR range.
Based on current progress in this area, the hypothesis of this study is that it is possible to calculate liming recommendations for oxisols based on the chemical characteristics estimated by PLSR, with previous application of variable selection by iPLS. Our objective was therefore to calculate the need for liming based on the chemical characteristics estimated by iPLS selection and PLSR modeling of specific wavelengths of soil reflectance.

Experimental Area
The experiment (Figure 1)  with relevant information in different spectral subdivisions, thus focusing on critical spectral regions and removing interference from other regions [39]. The iPLS algorithm selects a relevant subset of intervals (a subset of wavelengths) for each of the soil characteristics individually, which results in superior models compared to the use of the entire Vis-NIR spectrum, and hence decreases the complexity of the models by removing redundant information from the spectra [40]. When iPLS has split the spectra into different ranges, a PLSR validation model is applied to each of them. The interval that gives the lowest prediction error is modeled again by adding the remaining intervals, one by one, until the error no longer decreases [41].
Although many studies have reported the use of the Vis-NIR-SWIR range to predict soil characteristics, studies that have used a previous selection of variables for constructing PLSR models to estimate the necessary chemical characteristics are still rare in terms of calculating liming recommendations. These could make it possible in the future to develop cheaper sensors (or even cameras using specific wavelengths) which would not need to operate in the full Vis-NIR-SWIR range.
Based on current progress in this area, the hypothesis of this study is that it is possible to calculate liming recommendations for oxisols based on the chemical characteristics estimated by PLSR, with previous application of variable selection by iPLS. Our objective was therefore to calculate the need for liming based on the chemical characteristics estimated by iPLS selection and PLSR modeling of specific wavelengths of soil reflectance.

Experimental Area
The experiment (Figure 1) was carried out at the experimental farm of the Cooperativa Agroindustrial Mourãoense (COAMO), located in Campo Mourão, State of Paraná, Brazil (24°05′28″ S, 52°21′31″ W; altitude 605 m) and was evaluated between 2016 and 2019. According to the Köppen classification, the climate in this region is classified as (Cfa), a humid subtropical climate with hot summers. The annual average temperature is 20-21 °C (in the coldest quarter, the average temperature is 16-17 °C, and in the warmest, 27-28 According to the Köppen classification, the climate in this region is classified as (Cfa), a humid subtropical climate with hot summers. The annual average temperature is 20-21 • C (in the coldest quarter, the average temperature is 16-17 • C, and in the warmest, 27-28 • C), and the annual precipitation ranges from 1600 to 1800 mm [42]. The weather data (minimum and maximum temperatures and rainfall) for the experimental area between June 2016 and Remote Sens. 2022, 14, 1972 4 of 22 May 2019 are shown in Figure 2. During this period, the accumulated precipitation was approximately 5484 mm. °C), and the annual precipitation ranges from 1600 to 1800 mm [42]. The weather data (minimum and maximum temperatures and rainfall) for the experimental area between June 2016 and May 2019 are shown in Figure 2. During this period, the accumulated precipitation was approximately 5484 mm.
Before the experiment, the area had been managed with a no-till approach since 2008, with the following crops: oats (2008)  The soil within the experimental area was classified as a Latossolo Vermelho distroférrico típico [43], corresponding to a ferralsol [44] and an oxisol [45] of very clayey texture (740 g kg −1 clay). This soil analysis was carried out on May 18, 2016 (Table 1).

Implementation of the Experiment
The experimental design used was that of complete blocks with treatments distributed at random in a 2 × 4 crossed factorial scheme with three additional treatments and four blocks (giving a total of 11 treatments). These treatments included two types of limestone scale application (superficial and incorporated) and four base saturation levels (V%), corresponding to 44% (without correction), 60%, 70%, and 90%. Three additional treatments were applied that involved the use of phosphogypsum: V% 60 + 3.71 Mg ha −1 of phosphogypsum (60G50), V% 70 + 3.71 Mg ha −1 of phosphogypsum (70G50), and V% 70 + 7.42 Mg ha −1 of phosphogypsum (70G100). Each experimental plot was 12 m in length and 7 m in width, giving a total of 44 experimental plots of 84 m 2 .
The limestone application rates were calculated based on the elevation of base saturation recommended for the state of Paraná [5]. The application rates of phosphogypsum were calculated using the formula proposed by Sousa et al. [46], where the clay content was multiplied by 50 to obtain the recommended application rate and by 100 for twice the The soil within the experimental area was classified as a Latossolo Vermelho distroférrico típico [43], corresponding to a ferralsol [44] and an oxisol [45] of very clayey texture (740 g kg −1 clay). This soil analysis was carried out on 18 May 2016 (Table 1).

Implementation of the Experiment
The experimental design used was that of complete blocks with treatments distributed at random in a 2 × 4 crossed factorial scheme with three additional treatments and four blocks (giving a total of 11 treatments). These treatments included two types of limestone scale application (superficial and incorporated) and four base saturation levels (V%), corresponding to 44% (without correction), 60%, 70%, and 90%. Three additional treatments were applied that involved the use of phosphogypsum: V% 60 + 3.71 Mg ha −1 of phosphogypsum (60G50), V% 70 + 3.71 Mg ha −1 of phosphogypsum (70G50), and V% 70 + 7.42 Mg ha −1 of phosphogypsum (70G100). Each experimental plot was 12 m in length and 7 m in width, giving a total of 44 experimental plots of 84 m 2 .
The limestone application rates were calculated based on the elevation of base saturation recommended for the state of Paraná [5]. The application rates of phosphogypsum were calculated using the formula proposed by Sousa et al. [46], where the clay content was multiplied by 50 to obtain the recommended application rate and by 100 for twice the recommended application rate. The applied dolomitic limestone had a relative total neutralizing power (RTNP) of 80%. The treatments were applied on 10 June 2016, and the limestone application rates used are given in Table 1. Limestone was incorporated using a moldboard plow with an effective depth of incorporation of 0.20 m, followed by harrowing with a plow harrow containing 20 28" discs (Ø 0.71 m) and then light harrowing with a harrow containing 42 20" discs (Ø 0.51 m) to level the ground.

Soil Assessment and Chemical Analysis
Four soil assessments were performed at 5, 12, 24, and 36 months after the treatments were applied on 10 June 2016. The sampled profile was stratified into six layers: 0.00-0.05, 0.05-0.10, and 0.10-0.20 m (with a cutting shovel) and 0.20-0.30, 0.30-0.40, and 0.40-0.60 (with the aid of a Dutch auger). Two subsamples were collected within each experimental plot to obtain a composite sample. On each soil assessment date (5, 12, 24, and 36 months after applying the treatments), 264 samples were collected, giving a total of 1056 soil samples.
The soil samples were oven-dried at 40 • C, and passed through a 2 mm sieve [5]. Exchangeable Ca 2+ , Mg 2+ , and Al 3+ were extracted using a KCl 1 mol L −1 (1:10 v/v soil/solution). Exchangeable Ca 2+ and Mg 2+ contents were determined by atomic absorption spectrophotometry with an air-acetylene flame. Exchangeable Al 3+ was determined by titration (NaOH 0.0125 mol L −1 ). Exchangeable K + was extracted using Mehlich-1 solution (H 2 SO 4 0.0125 mol L −1 and HCl 0.05 mol L −1 at a ratio of 1:10 (v/v) soil/solution) and determined by flame emission spectroscopy. Phosphorus (P) was extracted by Mehlich-1 (H 2 SO 4 0.0125 mol L −1 and HCl 0.05 mol L −1 at a ratio of 1:10 (v/v) soil/solution) and was determined by UV-Vis spectrophotometry. Extraction of S-SO 4 2− was performed using calcium phosphate (500 mg L −1 P) in acetic acid 2.0 mol L −1 (1:2.5 v/v), and the concentration was determined using a UV/Vis spectrophotometer. Soil pH was determined in a CaCl 2 0.01 mol L −1 suspension (at a soil/solution ratio of 1:2.5 v/v). Following the procedures described by Pavan et al. [47], the potential acidity was determined using SMP buffer solution. Except for the potential acidity, the other methodologies were those of Silva et al. [48].
The sum of bases (SB) was calculated by adding the Ca 2+ , Mg 2+ , and K + contents. In naturally acidic and intensely leached soils, the values for Na + and NH 4 + are very low (practically negligible) and these were therefore excluded from the calculation of BS. The effective cation-exchange capacity (CECe = SB + Al 3+ ), cation exchange capacity at pH 7 (CEC = SB + H+Al), base saturation (BS = (SB/CEC) × 100) and aluminum saturation (m% = (Al 3+ /CECe) × 100) were also calculated. The chemical characteristics of the soil were subjected to a descriptive statistical analysis.

Acquisition and Processing of Spectral Data
Soil samples were placed in Petri dishes 9 cm in diameter, and spectral readings were acquired in the spectral range 350-2500 nm (Vis-NIR) with a Fieldspec 3 Jr spectroradiometer (ASD Inc, Boulder, CO, USA), with a spectral resolution of 3 nm between 350 and 700 nm and 10 nm between 700 nm and 2500 nm.
The geometry used to collect spectral data was adopted from Demattê et al. [49] and Nanni et al. [50]: a white Spectralon plate was used as a reflectance standard with 100% reflectance, as per Labsphere Reflectance Calibration Laboratory (LABSPHERE Inc., North Sutton, NH, USA); and a 650 W halogen lamp was used as a light source, positioned 0.35 m from the sample and with a zenith angle of 30 • . The optical sensor was placed in an upright position 0.08 m away from the sample. The soil reflectance was measured at the center of the sample, as shown in Figure 3.
trophotometry with an air-acetylene flame. Exchangeable Al 3+ was determined by titration (NaOH 0.0125 mol L −1 ). Exchangeable K + was extracted using Mehlich-1 solution (H2SO4 0.0125 mol L −1 and HCl 0.05 mol L −1 at a ratio of 1:10 (v/v) soil/solution) and determined by flame emission spectroscopy. Phosphorus (P) was extracted by Mehlich-1 (H2SO4 0.0125 mol L −1 and HCl 0.05 mol L −1 at a ratio of 1:10 (v/v) soil/solution) and was determined by UV-Vis spectrophotometry. Extraction of S-SO4 2− was performed using calcium phosphate (500 mg L −1 P) in acetic acid 2.0 mol L −1 (1:2.5 v/v), and the concentration was determined using a UV/Vis spectrophotometer. Soil pH was determined in a CaCl2 0.01 mol L −1 suspension (at a soil/solution ratio of 1:2.5 v/v). Following the procedures described by Pavan et al. [47], the potential acidity was determined using SMP buffer solution. Except for the potential acidity, the other methodologies were those of Silva et al. [48].
The sum of bases (SB) was calculated by adding the Ca 2+ , Mg 2+ , and K + contents. In naturally acidic and intensely leached soils, the values for Na + and NH4 + are very low (practically negligible) and these were therefore excluded from the calculation of BS. The effective cation-exchange capacity (CECe = SB + Al 3+ ), cation exchange capacity at pH 7 (CEC = SB + H+Al), base saturation (BS = (SB/CEC) × 100) and aluminum saturation (m% = (Al 3+ /CECe) × 100) were also calculated. The chemical characteristics of the soil were subjected to a descriptive statistical analysis.

Acquisition and Processing of Spectral Data
Soil samples were placed in Petri dishes 9 cm in diameter, and spectral readings were acquired in the spectral range 350-2500 nm (Vis-NIR) with a Fieldspec 3 Jr spectroradiometer (ASD Inc, Boulder, CO, USA), with a spectral resolution of 3 nm between 350 and 700 nm and 10 nm between 700 and 2500 nm.
The geometry used to collect spectral data was adopted from Demattê et al. [49] and Nanni et al. [50]: a white Spectralon plate was used as a reflectance standard with 100% reflectance, as per Labsphere Reflectance Calibration Laboratory (LABSPHERE Inc., North Sutton, NH, USA); and a 650 W halogen lamp was used as a light source, positioned 0.35 m from the sample and with a zenith angle of 30°. The optical sensor was placed in an upright position 0.08 m away from the sample. The soil reflectance was measured at the center of the sample, as shown in Figure 3. For each soil sample, three readings were performed by rotating the Petri dish through 120 • [50,51], to avoid any bias due to the arrangement of particles. The reflectance factor was obtained, which expresses the ratio between the radiant flux of the soil sample and the radiant flux reflected by the reference standard under the same geometry and lighting conditions [52].
The collected spectra were corrected using the "splice correction" routine of the ViewSpecPro 6.0 software package (ASD Inc, Boulder, CO, USA). This routine aims to eliminate steps (or offsets) in the original spectra caused by phase differences between the measurements of the three detectors used with the spectroradiometer: VNIR (350-1000 nm), SWIR 1 (1000-1800 nm), and SWIR 2 (1800-2500 nm). Parts of the spectral curve generated for each Remote Sens. 2022, 14,1972 7 of 22 sample were removed due to noise, mainly at the beginning and end of the spectral curve. The final spectrum band used for analysis was 400-2400 nm, with a total of 2001 bands.
Several techniques were tested for the preprocessing of the curves, to eliminate or smooth the curves and reduce interference related to the sensor, in order to give better calibration and prediction results. This preprocessing is usually performed in chemometrics with spectral data to reduce optical interference in spectra unrelated to the chemical composition of the sample [22,28,50]. Unscrambler X 10.4.1 software (Camo Software, Oslo, Norway) was used for this step.
Of the preprocessing methods tested at this stage, the Savitzky-Golay approach [53] resulted in the best performance of the model, with a first-order polynomial and a five-point segment about the central point. This technique is a digital polynomial filter in which each value in the signal series is replaced with a new value obtained from a polynomial fit to successive subsets of adjacent data points. The fit is performed using linear least squares (with 2n + 1 neighboring points), where n can be equal to or greater than the order of the polynomial. The more neighbors used in the averaging process, the smoother the signal becomes. Least squares smoothing suppresses noise while maintaining signal information [54].
SAS software (SAS Institute Inc., Cary, NC, USA) was used for this discriminant analysis, using the STEPDISC and DISCRIM procedures. The STEPDISC procedure was executed by the FORWARD method at 5% probability, and began with no variables in the model. PROC STEPDISC inserted the variable that contributed the most to the model's discriminatory power at each step, as measured by Wilks' lambda. The process ended when none of the unselected variables met the entry criteria [55].
The variables that best explained the differences among the treatments and the sampled layers were selected to compose the discriminant model. This procedure generated discriminant models through linear combinations of predictor variables (wavelengths) [55].
Using SAS software, a simulation was performed to test the quality of the models obtained. In this simulation, 70% of the spectral curves were randomly selected and used to construct a new discriminant model, which was later tested using the remaining 30% of the spectral curves. This simulation process was repeated 50 consecutive times [50,56,57].

Partial Least Squares Regression
After the discriminant analysis, variables were selected based on the iPLS method, using PLS Toolbox 8.9.2 (Eigenvector Research Inc., Wenatchee, WA, USA), a plugin for MatLab software (MathWorks Inc., Natick, MA, USA). PLS Toolbox was operated in "forward" mode, where intervals are successively included in the analysis. The program was set to automatic mode, in which the iPLS algorithm continued to add (or remove) gaps until the RMSE no longer improved. The size of the intervals was defined as "one" (1 = one variable per interval), meaning that the number of selected intervals was equal to the number of selected variables [58].
PLSR was performed after selecting the wavelengths, using Unscrambler X 10.4.1 software (Camo Software, Oslo, Norway) with the NIPALS algorithm [59]. PLSR was conducted in three steps, consisting of calibration, validation, and prediction. The model quality parameters obtained in the validation steps indicated whether the models were suitable for use in predicting the chemical characteristics of the soil.
The database was randomized using Unscrambler X 10.4.1 software (Camo Software, Oslo, Norway), and was divided into two groups: the calibration and validation set (containing 75% of the samples) and the prediction set (containing the remaining 25% of the samples). PLSR models were created to quantify the following chemical characteristics of the soil: Ca 2+ , Mg 2+ , K + , P, SO 4 2− , pH, Al 3+ , H+Al, SB, CECe, CEC, m%, and BS. In this study, we used the following statistical parameters as indicators of the goodness of fit of the regression models: the coefficient of determination (R 2 ), root mean square error (RMSE), systematic error (BIAS), relative percentage deviation (RPD), and the ratio of performance to interquartile range (RPIQ).
R 2 values were classified according to Wenjun et al. [60], Williams [61], Viscarra Rossel et al. [34], and Saeys et al. [62], where: R 2 < 0.50 represented models with poor predictions, which were only able to distinguish between high and low values; R 2 between 0.50 and 0.65 indicated models with moderate predictions that could be used for evaluation and correlation; R 2 between 0.65 and 0.80 represented good prediction models that allowed for quantitative predictions; R 2 between 0.80 and 0.90 indicated very good quantitative prediction models; and R 2 ≥ 0.90 represented excellent quantitative prediction models.
The RPD values were interpreted following Chang et al. [63] and the results were divided into three categories: highly reliable models with RPD > 2 (which can be used for quantitative analysis and making predictions); reliable models with 1.4 < RPD < 2; and unreliable models for predictions with RPD < 1.4.
Four categories of models were defined based on the RPIQ values: an RPIQ score of between 2.02 and 2.70 indicates a poor model, which can only distinguish between high and low; an RPIQ of between 2.70 and 3.37 represents a regular model, in which quantitative predictions are possible; an RPIQ of between 3.37 and 4.05 represents a good quantitative model; and RPIQ > 4.05 implies an excellent model [64].

Spectral-Based Method for Liming Recommendation
The chemical characteristics predicted by the PLSR calibration models were used for liming recommendation (LR) using the base saturation method, as follows: where BS 2 is the desired base saturation (70%), BS 1 is the current base saturation; RPTN is the relative power of total neutralization (100%); and CEC is the cation exchange capacity at pH 7 (cmol c dm −3 ), according to SBCS/NEPAR [5].
A correlation analysis was then performed between the estimated liming recommendation using the characteristics obtained by the PLSR models (based only on the prediction dataset), and the calculated liming recommendation, using the characteristics obtained by conventional laboratory techniques.  (Table 2). This was as expected, due to the large variability of the chemical characteristics in the soil profile, particularly when compared to the surface layers (which had higher SB, CECe, and BS content) and subsurface layers (with lower SB, CECe, and BS content). The smallest coefficients of variance were obtained for H+Al and CEC, as this oxisol showed only slight mineralogical and texture variations along the profile.

Soil Chemical Characteristics
The highest values of the coefficient of variance (CV) were observed for P (118.68% in the second soil assessment), S-SO 4 2− (160.89% in the first soil assessment), and Al 3+ (93.89% in the first soil assessment). The lowest CVs were observed for pH (7.44% in the fourth soil assessment), H+Al (16.83% in the fourth soil assessment), and CEC (18.32% in the first soil assessment). * The data follow the normal Kolmogorov-Smirnov distribution (p-value ≤ 0.01). Ca 2+ , Mg 2+ and Al 3+ contents were extracted using KCl 1.0 mol L −1 ; K + and P contents were extracted using Mehlich-1; pH CaCl 2 (0.01 mol L −1 ) with a soil:solution ratio of 1:2.5; H+Al was extracted using the SMP method; CECe is the effective cation exchange capacity; CEC is the cation exchange capacity at pH 7; m is the aluminum saturation; BS is the base saturation; CV is the coefficient of variance; std. dev. is the standard deviation.
The minimum and maximum levels of Ca 2+ increased after the application of the treatments (5 to 24 months after application). In the fourth soil assessment (36 months after application), the Ca 2+ content was found to have decreased. The Mg 2+ content did not vary significantly among the different soil assessments. The SB increased after the treatments were applied, but the values in the fourth soil assessment were similar to those in the first soil assessment. This can be explained by the extraction of bases by the cultures and their leaching. Similar behavior was observed for CECe, CEC, and BS (Table 2).

Spectral Modeling: Discriminant Analysis
The selection of variables for the discriminant analysis performed among treatments, using PROC STEPDISC at 5% probability, did not select any variable. Thus, the discriminant analysis between treatments was interrupted. The discriminant analysis performed to discriminate the different soil layers sampled was completed. The accuracy of the discriminant analysis of the four soil assessments is illustrated in Figure 4.
It was observed that the 0.00-0.05 m layer had the highest accuracy in terms of the classification of samples among soil assessments, ranging from 96.27-100% in the calibration step and 92.46-100% in the validation phase. The discriminant model of the 0.00-0.05 m layer in the fourth soil assessment reached 100% accuracy in both the calibration and validation steps. This separation of the 0.00-0.05 m layer was probably due to the higher content of organic matter and SB, which are usually higher in the first centimeters of the soil profile, thus distinguishing it from the other layers.
The 0.40-0.60 m layer also obtained a high accuracy percentage, ranging from 84.58-90.63% in the calibration step and 74.06-86.49% in the validation step for the different soil assessments. This difference was probably due to the lower content of organic matter and SB in the 0.40-0.60 m layers.  The results from the discriminant analysis indicate that the differences among th sampled soil layers had a larger influence on the PLSR models than the differences be tween the applied treatments.

Spectral Modeling: iPLS
An example of the selection of iPLS variables from the reflectance curves (400-240 nm) for CEC prediction in the second soil assessment can be seen in Figure 5. A validatio PLSR model was applied to each of the intervals of the reflectance curves, and after var able selection, a decrease was observed in the RMSE from 0.608 (red line, using all varia bles) to 0.34 cmolc dm −3 (blue line, using only the selected variables). The intermediate soil layers between 0.05 and 0.40 m had an accuracy ranging from 65.59-95.22% in the calibration step and 51.74-85.14% in the validation step. The lower accuracy of these discriminant models is due to the larger variation observed among the spectral curves of these layers.
It was observed that the models had a high percentage accuracy, indicating that this statistical technique is valid and sensitive in terms of the differentiation of the layers. In the fourth soil assessment, even in the intermediate layers, the percentage of success was higher than 89.1% in the calibration step and 68.4% in the validation step ( Figure 4).
The results from the discriminant analysis indicate that the differences among the sampled soil layers had a larger influence on the PLSR models than the differences between the applied treatments.

Spectral Modeling: iPLS
An example of the selection of iPLS variables from the reflectance curves (400-2400 nm) for CEC prediction in the second soil assessment can be seen in Figure 5. A validation PLSR model was applied to each of the intervals of the reflectance curves, and after variable selection, a decrease was observed in the RMSE from 0.608 (red line, using all variables) to 0.34 cmol c dm −3 (blue line, using only the selected variables).
Remote Sens. 2022, 14, x FOR PEER REVIEW 12 Figure 5. Example of selection of iPLS variables for CEC in the second soil assessment. The Y unit is in reflectance factor for the mean sample, and for RMSECV is in cmolc dm −3 .
The wavelengths selected by the iPLS algorithm for each soil characteristic are sho in Figure 6. Note that for each chemical characteristic, different wavelengths were sele in the four assessments due to the variability in these characteristics (Table 2) over assessment period. For characteristics with lower variability, such as CEC, the sele wavelengths were more similar among the soil assessments. Furthermore, for the Ca 2+ Mg 2+ content, the selected wavelengths of close to 1.690 and 2.300 nm are related to applied dolomitic limestone (CaMg (CO3)2) [65]. The wavelengths selected by the iPLS algorithm for each soil characteristic are shown in Figure 6. Note that for each chemical characteristic, different wavelengths were selected in the four assessments due to the variability in these characteristics (Table 2) over the assessment period. For characteristics with lower variability, such as CEC, the selected wavelengths were more similar among the soil assessments. Furthermore, for the Ca 2+ and Mg 2+ content, the selected wavelengths of close to 1.690 and 2.300 nm are related to the applied dolomitic limestone (CaMg (CO 3 ) 2 ) [65].

Spectral Modeling: PLSR
The statistical parameters of the PLSR models developed in the first soil assessment are shown in Table 3. The characteristics of Ca 2+ and K + content, P, SB, CECe, CEC, and BS had values for R 2 ranging from 0.70 to 0.86 in the calibration and validation steps, and were therefore classified as good or very good prediction models. In the prediction step, the values of R 2 ranged from 0.63 to 0.82, giving models with moderate to very good quantitative prediction capabilities. Our values were similar to those found by Rodrigues et al. [28] for Ca 2+ , K + , and P (with R 2 above 0.76 in the calibration and validation steps), and were higher than those found by Veum et al. [66] for Ca 2+ , K + , and P.
Of the PLSR models developed in the first soil assessment (Table 3), those for K + , P, and CECe stood out, with values of R 2 of above 0.81. For P, the prediction R 2 was 0.82, while for K + , the value was 0.81, both of which were higher than those reported by Zornoza et al. [22]. With R 2 values of above 0.80, these models provide a very good estimate of the characteristics. For Mg 2+ , it was not possible to perform the prediction step, possibly due to the low value of R 2 in both the calibration and the validation steps. For S-SO 4 2− , Al 3+ , and m%, it was not possible to develop good models, as they could not explain the variability of the data. Remote Sens. 2022, 14, x FOR PEER REVIEW 13 of 25

Spectral Modeling: PLSR
The statistical parameters of the PLSR models developed in the first soil assessment are shown in Table 3. The characteristics of Ca 2+ and K + content, P, SB, CECe, CEC, and BS had values for R 2 ranging from 0.70 to 0.86 in the calibration and validation steps, and were therefore classified as good or very good prediction models. In the prediction step, the values of R 2 ranged from 0.63 to 0.82, giving models with moderate to very good quantitative prediction capabilities. Our values were similar to those found by Rodrigues et al. [28] for Ca 2+ , K + , and P (with R 2 above 0.76 in the calibration and validation steps), and were higher than those found by Veum et al. [66] for Ca 2+ , K + , and P.
Of the PLSR models developed in the first soil assessment (Table 3), those for K + , P, and CECe stood out, with values of R 2 of above 0.81. For P, the prediction R 2 was 0.82, while for K + , the value was 0.81, both of which were higher than those reported by Zornoza et al. [22]. With R 2 values of above 0.80, these models provide a very good estimate of the characteristics. For Mg 2+ , it was not possible to perform the prediction step, possibly due to the low value of R 2 in both the calibration and the validation steps. For S-SO4 2− , Al 3+ , For the RMSE value in the validation step (Table 3), there was a slight variation compared to the calibration step, indicating models with good predictive capacity. These low RMSE values agree with the coefficients of determination, which were classified as good and very good. The same results were found in the other evaluations (Tables 4-6).
Based on the RPD results (Table 3), the PLSR models in the calibration phase were classified as highly reliable for Ca 2+ , K + , SB, CECe and CEC, reliable for Mg 2+ , P, pH, and BS, and unreliable for S-SO 4 2− , Al 3+ , H+Al, and m%. Our RPD values were higher than those reported by Veum et al. [66] for Ca 2+ , Mg 2+ , K + , P, and pH.
Based on the RPIQ values from the first evaluation (Table 3), the PLSR models were classified as poor for Ca 2+ , pH, and BS (in which only high and low values are distinguishable), regular for P, SB, and CECe (in which quantitative predictions are possible), and good for CEC (representing a good quantitative model). These results show that the RPIQ parameter proved to be more rigorous than the other parameters for the classification of the developed models; even so, it was still possible to obtain regular and good models. In particular, the PLSR model for CEC was classified as good based on the RPIQ. Both the R 2 and the RPD demonstrated that this was a good and highly reliable model. Ca 2+ , Mg 2+ , and Al 3+ contents were extracted using KCl 1.0 mol L −1 ; K + and P contents were extracted using Mehlich-1; pH CaCl 2 (0.01 mol L −1 ) at a soil:solution ratio of 1:2.5; H+Al was extracted using the SMP method; CECe is the effective cation exchange capacity; CEC is the cation exchange capacity at pH 7; m indicates the aluminum saturation; BS is the base saturation; R 2 is the coefficient of determination; RMSE is the root mean square error; BIAS is the systematic error; RPD represents the relative percentage deviation; RPIQ is the ratio of performance to interquartile range. Ca 2+ , Mg 2+ , and Al 3+ contents were extracted using KCl 1.0 mol L −1 ; K + and P contents were extracted using Mehlich-1; pH CaCl 2 (0.01 mol L −1 ) at a soil:solution ratio of 1:2.5; H+Al was extracted using the SMP method; CECe is the effective cation exchange capacity; CEC is the cation exchange capacity at pH 7; m represents the aluminum saturation; BS is the base saturation; R 2 is the coefficient of determination; RMSE is the root mean square error; BIAS is the systematic error; RPD is the relative percentage deviation; and RPIQ is the ratio of performance to interquartile range.
In general, an improvement was seen in most PLSR models with iPLS selection (Table 3) compared to PLSR models without iPLS selection (Table A1), although no improvement was observed for S-SO4 2− , pH, Al +3 , H+Al, and m%. The R 2 values in the calibration step were higher by about 16% for Mg 2+ and P, 9% for CECe and CEC, and 4.5% for SB and BS (Table A1). For Ca 2+ and K + , no increase in R 2 was observed. There was an average reduction of 16% in the calibration RMSE with iPLS for Ca 2+ , Mg 2+ , K + , P, SB, CECe, CEC, and BS compared to the PLSR models without iPLS selection; in particular, the RMSE values for K+, P and CEC showed reductions of 20%, 33% and 24%, respectively. In the validation phase, there was an average increase of 10% in the R 2 values for Ca 2+ , Mg 2+ , K + , P, SB, CECe, CEC, and BS, where the values for Mg 2+ , P, CECe, and CEC showed increases of 16%, 16%, 11% and 11%, respectively. In the same way as for the calibration step, there was an average reduction of 16% in the RMSE for the validation step with iPLS, and in particular, the RMSE values for P, SB, CECe, and CEC showed reductions of 31%, 22%, 18%, and 27%, respectively. K+ and P stood out in the prediction phase, as they showed increases in R 2 of 15% and 27% and reductions in the RMSE of 17% and 53%, respectively. The prediction RMSE for CEC was also noteworthy, with a reduction of 33%. Ca 2+ , Mg 2+ , and Al 3+ contents were extracted using KCl 1.0 mol L −1 ; K + and P contents were extracted using Mehlich-1; pH CaCl 2 (0.01 mol L −1 ) at a soil:solution ratio of 1:2.5; H+Al was extracted using the SMP method; CECe is the effective cation exchange capacity; CEC is the cation exchange capacity at pH 7; m represents the aluminum saturation; BS is the base saturation; R 2 is the coefficient of determination; RMSE is the root mean square error; BIAS is the systematic error; RPD is the relative percentage deviation; and RPIQ is the ratio of performance to interquartile range. Ca 2+ , Mg 2+ , and Al 3+ contents were extracted using KCl 1.0 mol L −1 ; K + and P contents were extracted using Mehlich-1; pH CaCl 2 (0.01 mol L −1 ) at a soil:solution ratio of 1:2.5; H+Al was extracted using the SMP method; CECe is the effective cation exchange capacity; CEC is the cation exchange capacity at pH 7; m represents the aluminum saturation; BS is the base saturation; R 2 is the coefficient of determination; RMSE is the root mean square error; BIAS is the systematic error; RPD is the relative percentage deviation; and RPIQ is the ratio of performance to interquartile range.
In the second soil assessment (Table 4), the PLSR models of Ca 2+ , K + , P, SB, CECe, CEC, and BS achieved the highest values of R 2 , ranging from 0.73 to 0.90 in the calibration and validation steps; in the prediction step, the values of R 2 ranged from 0.67 to 0.91. This range of R 2 values indicates that the models have good to excellent quantitative prediction capabilities. These values were higher than those reported by Pinheiro et al. [26] and Veum et al. [66]. The results from the second soil assessment were also superior to those of the first assessment (Table 3).
In addition, the PLSR models developed for the prediction of SB, CECe, and CEC achieved R 2 values of above 0.85; the value for CEC reached 0.91 (Table 4), and the model had lower numbers of variables (93) and factors (three) compared to those of SB and CECe. For S-SO 4 2− , Al 3+ , H+Al, and m%, it was not possible to develop good models in the second soil assessment (Table 4), as they could not explain the data variability.
Based on the RPD values (Table 4), the PLSR models of the calibration phase were classified as highly reliable for Ca 2+ , P, SB, CECe, CEC and BS, reliable for Mg 2+ , K + and pH, and unreliable for S-SO 4 2− , Al 3+ , H+Al, and m%. Note that there was an improvement in the RPD for BS compared to the first evaluation, which was accompanied by an increase in R 2 values in all phases.
Based on the RPIQ values for the second evaluation (Table 4), the PLSR models were classified as poor for pH (meaning that only high and low values were distinguishable), regular for Ca 2+ , P, SB, and BS (meaning that quantitative predictions were possible), good for CECe, and excellent for CEC. These RPIQ values show an improvement in the predictive capacity of the models compared to the first evaluation, mainly for Ca 2+ , CECe, CEC, and BS. The PLSR model for CEC, in particular, was classified as excellent based on the RPIQ, and was shown to be a very good and highly reliable model based on both the R 2 and the RPD.
The better accuracy in the second evaluation can be explained by the longer reaction time for the limestone and phosphogypsum, which allowed for higher values of SB, CECe, CEC, and BS ( Table 2). As a result of this greater amplitude, more robust models were obtained with a better fit.
In general, there was an improvement in most PLSR models with iPLS selection (Table 4) compared to models without iPLS selection (Table A1). However, again, no improvements were observed for S-SO 4 2− , pH, Al +3 , or H+Al, and these results were repeated in the other evaluations. After iPLS selection, the R 2 values in the calibration step were increased by 22% for Mg 2+ , 18% for K + , and 13% for SB and BS. Regarding the calibration RMSE, there was an average reduction of 28% for Ca 2+ , Mg 2+ , K + , P, SB, CECe, CEC, and BS compared to PLSR models without iPLS selection, where the results for K + , SB, CEC, and BS showed reductions of 30%, 43%, 36% and 29%, respectively. In the validation phase, there was an average increase of 14% in the R 2 values for Ca 2+ , Mg 2+ , K + , P, SB, CECe, CEC, and BS; the R 2 values for Mg 2+ , K + , SB, and BS showed increases of 26%, 19%, 15% and 18%, respectively. In the validation stage, there was an average reduction of 30% of the RMSE with iPLS, with the RMSE for SB, CECe, and CEC showing reductions of 46%, 35%, and 30%, respectively. In the prediction phase, the results for SB, CECe, and BS showed increases of 24%, 19%, and 35% in the values of R 2 , and reductions in the RMSE of 38%, 36%, and 26%, respectively. Another noteworthy result was the prediction RMSE for CEC, with a reduction of 22%.
PLSR models with similar adjustments to the second soil assessment were also found for the third (Table 5) and fourth soil assessments (Table 6), two and three years after the treatment application (2016), respectively. In the fourth soil assessment, the PLSR models of Ca 2+ , K + , SB, CECe, CEC, and BS achieved the highest values of R 2 , ranging from 0.83 to 0.92 in the calibration and validation steps and from 0.70 to 0.90 for the prediction step. This range of R 2 values indicates that these models have a good or very good capacity for quantitative prediction. These values were higher than those reported by Pinheiro et al. [26] and Terra et al. [25], but lower than those found by Rodrigues et al. [28] for Ca 2+ and K + , for all three steps (calibration, validation, and prediction). However, it was not possible to develop good models for S-SO 4 2− , pH, Al 3+ , H+Al, and m% for the calibration, validation, and prediction steps of the third and fourth soil assessments (Tables 5 and 6). The prediction models for CEC and SB for the fourth soil assessment (Table 6) were notable, with R 2 values of 0.89 and 0.90, respectively, indicating a very good or excellent capacity for quantitative prediction.
For the RPD values, the results for the third assessment were similar to those of the second assessment ( Table 5). The PLSR models for the calibration phase were classified as highly reliable for Ca 2+ , K + , P, SB, CECe, and CEC, reliable for Mg 2+ , pH and BS, and unreliable for S-SO 4 2− , Al 3+ , H+Al, and m%. Based on the RPIQ values of the third assessment (Table 5), the PLSR models were classified as poor for Mg 2+ (meaning that only high and low values were distinguishable), regular for Ca 2+ , K + , P, pH, and BS (meaning that quantitative predictions were possible), good for SB and CECe (meaning that these were good quantitative models), and excellent for CEC. There was an improvement in the predictive capacity of the SB model compared with the second evaluation, and this was classified as a good quantitative model. The classification of the PLSR model for CEC remained excellent, based on the RPIQ.
A comparison between PLSR models with iPLS selection (Table 5) and PLSR models without iPLS selection (Table A1) indicates that the R 2 values for the calibration stage increased by 7% for Ca 2+ and Mg 2+ , and 8% for SB and CEC. In terms of the calibration RMSE, there was an average reduction of 14% for Ca 2+ , Mg 2+ , K + , P, SB, CECe, CEC, and BS compared to the PLSR models without iPLS selection; in particular, the SB and CEC models showed reductions of 21% and 38%, respectively. In the validation phase, there was an average increase of 8% in the R 2 values for Ca 2+ , Mg 2+ , K + , P, SB, CECe, CEC, and BS. The R 2 values for Mg 2+ and SB were notable, with increases of 16% and 10%, respectively. Regarding the validation RMSE, there was an average reduction of 18% in the RMSE with iPLS, and in particular, the RMSE values for Ca 2+ , SB, and CEC showed reductions of 21%, 26%, and 43%, respectively. In the prediction phase, the results for K + were notable, and showed an increase of 16% in R 2 and a reduction in RMSE of 29%. However, in the prediction phase, the PLSR models with iPLS developed in the third evaluation were not improved for most features, compared to the PLSR models without iPLS.
The RPD values found in the fourth evaluation (Table 6) indicate that the PLSR models of the calibration phase can be classified as highly reliable for Ca 2+ , K + , P, SB, CECe, CEC, and BS, reliable for Mg 2+ and pH, and unreliable for S-SO 4 2− , Al 3+ , H+Al, and m. The classification for the BS model was elevated to highly reliable, demonstrating the evolution of the BS parameters throughout the evaluations.
Based on the RPIQ values for the fourth evaluation (Table 6), the PLSR models were classified as poor for Mg 2+ , pH, Al 3+ , H+Al (meaning that only high and low values are distinguishable), regular for K + (meaning that quantitative predictions are possible), good for Ca 2+ and P, and excellent for SB, CECe, CEC, and BS. The predictive capacity of the SB, CECe, and BS models improved compared to the other assessments, and these were classified as excellent quantitative models.
The PLSR model for CEC remained excellent based on the RPIQ values over the last three evaluations, a result that can be explained by the association of CEC with mineralogy and soil granulometry [2]. Numerous studies have been performed using Vis-NIR-SWIR reflectance spectroscopy to predict the physical characteristics of soil, such as granulometry [50,67]. The relationship of Vis-NIR-SWIR spectroscopy with the characteristics and physics of soil is already well established; however, works that include exchangeable ions, SB, CECe, CEC, and BS are still scarce, and mainly involve variable selection.
Since both the CEC and BS are necessary to calculate liming recommendations using the BS method [5], the results for R 2 , RPD, and RPIQ show that it is possible to calculate liming recommendations based on the values estimated by the CEC and BS PLSR models.
Once again, an improvement was noted in most PLSR models with iPLS selection (Table 6) compared to models without iPLS selection (Table A1). After iPLS selection, the R 2 values at the calibration stage were on average 4% higher, and in particular, the results for Ca 2+ and BS showed increases of 7% and 9%, respectively. There was an average reduction of 18% in the calibration RMSE for models with iPLS for Ca 2+ , K + , SB, CECe, CEC, and BS compared to the PLSR models without iPLS selection; notable results were seen for Ca 2+ and BS, which showed reductions of 39% and 27%, respectively. In the validation phase, there was an average increase of 4% in the R 2 values for Ca 2+ , Mg 2+ , K + , P, SB, CECe, CEC, and BS, and the values of R 2 for Mg 2+ and BS increased by 7% and 8%, respectively. Again, in the validation phase, there was an average reduction of 13% in the RMSE for models with iPLS, where the values for Ca 2+ , K + , and BS showed reductions of 31%, 14%, and 17%, respectively. The results for K+ and BS were notable in the prediction phase, as these showed increases of 13% and 12% in R 2 and reductions of 11% and 48% in RMSE, respectively. The prediction results for the RMSE for CEC were also significant, with a reduction of 50%. For Ca 2+ and Mg 2+, there were increases in RMSE of 15% and 4%, respectively, in the models with iPLS. These results indicate a greater reduction in the RMSE in the models with iPLS, particularly for the CEC and BS, which are characteristics used for liming recommendation in the BS method [5].
In most cases, the reduction in the RMSE values obtained by the PLSR models with iPLS (in the four evaluations) improved the predictive capacity of these models, since the RMSE measures the differences between the estimated values and the real values, and hence quantifies the accuracy of the model. Reductions in the RMSE values therefore gave rise to better accuracy for these models.
Furthermore, for all the characteristics studied for the four soil samples (Tables 3-6), the systematic error values (BIAS) in the calibration and validation phases were close to zero; these results demonstrate the very low bias of the model in regard to the estimated characteristics [68]. Willmott [69] states that PLSR models are considered to be good when the BIAS approaches zero, similarly to the results in the prediction phase for characteristics with R 2 above 0.6 (Tables 3-6).
Thus, there was a gradual evolution of the statistical parameters from the first to the fourth assessment; this may have been because as the treatments (limestone and phosphogypsum) reacted with the soil over the years, greater amplitudes in the chemical characteristics of the soil were generated, allowing us to find better models. Figure 7 shows the results of a correlation analysis between the estimated liming recommendation using the characteristics obtained by the PLSR models, and the calculated liming recommendation using the characteristics obtained by conventional laboratory techniques. The first soil assessment showed the lowest correlation (r = 0.86), and the correlation values increased up to the fourth soil assessment, for which the correlation coefficient was 0.95. These results therefore show that over the four soil assessments (5,12,24, and 36 months after the application of treatments), the liming recommendations based on the chemical characteristics obtained by the PLSR models became progressively more correlated with the liming recommendations calculated based on the chemical characteristics obtained by conventional laboratory techniques. This was probably due to the limestone reaction time. Several studies of Brazilian soils have demonstrated that three years is sufficient for all of the applied limestone to undergo reaction [6,70,71]. Thus, for our very clayey Oxisol, liming recommendations based on multivariate statistical modeling of soil reflectance should preferably be made three years after the last liming, as in this case, the prediction models for CEC and BS had values of R 2 for prediction of 0.90 (Table 6).

Spectral-Based Method for Liming Recommendation
reaction time. Several studies of Brazilian soils have demonstrated that three years is suf ficient for all of the applied limestone to undergo reaction [6,70,71]. Thus, for our ver clayey Oxisol, liming recommendations based on multivariate statistical modeling of soi reflectance should preferably be made three years after the last liming, as in this case, th prediction models for CEC and BS had values of R 2 for prediction of 0.90 (Table 6).

Figure 7.
Correlations between the calculated liming recommendation (LR) (based on chemical char acteristics obtained by conventional laboratory techniques) and the liming recommendation (LR estimated using PLSR models (based on the chemical characteristics obtained by the PLSR models for four soil assessments (1 st = first assessment; 2 nd = second assessment; 3 rd = third assessment; 4 th fourth assessment).

Conclusions
Our models were shown to discriminate between the layers sampled from the soi profile with high accuracy, although models discriminating among treatments could no be developed.
The use of iPLS selection gave rise to an increase in R 2 and a reduction in the RMSE values for most PLSR models developed for all four soil assessments, which improved th accuracy of these models.
We were able to estimate the soil chemical characteristics by selecting wavelength using the iPLS method and PLSR modeling, and, in particular, for Ca 2+ , SB, CTCe, CTC and BS. The statistical parameters of these models gradually improved until the fourt assessment. The values of R 2 ranged from 0.83 to 0.92 in the calibration and validation steps, and from 0.84 to 0.90 in the prediction step (in the fourth assessment).
The PLSR models of CEC and BS made it possible to calculate a liming recommenda tion based on the chemical characteristics estimated by the PLSR models. This estimat was positively correlated (r > 0.86) with the liming recommendation calculated based o Figure 7. Correlations between the calculated liming recommendation (LR) (based on chemical characteristics obtained by conventional laboratory techniques) and the liming recommendation (LR) estimated using PLSR models (based on the chemical characteristics obtained by the PLSR models) for four soil assessments (1 st = first assessment; 2 nd = second assessment; 3 rd = third assessment; 4 th = fourth assessment).

Conclusions
Our models were shown to discriminate between the layers sampled from the soil profile with high accuracy, although models discriminating among treatments could not be developed.
The use of iPLS selection gave rise to an increase in R 2 and a reduction in the RMSE values for most PLSR models developed for all four soil assessments, which improved the accuracy of these models.
We were able to estimate the soil chemical characteristics by selecting wavelengths using the iPLS method and PLSR modeling, and, in particular, for Ca 2+ , SB, CTCe, CTC, and BS. The statistical parameters of these models gradually improved until the fourth assessment. The values of R 2 ranged from 0.83 to 0.92 in the calibration and validation steps, and from 0.84 to 0.90 in the prediction step (in the fourth assessment).
The PLSR models of CEC and BS made it possible to calculate a liming recommendation based on the chemical characteristics estimated by the PLSR models. This estimate was positively correlated (r > 0.86) with the liming recommendation calculated based on the chemical characteristics obtained by conventional laboratory techniques. The fourth soil assessment yielded the best correlation coefficient (r = 0.95).
Our results demonstrate that the use of Vis-NIR-SWIR spectroscopy in conjunction with the selection of iPLS variables technique has strong potential in terms of predicting soil chemical characteristics, and can complement conventional analyses of these characteristics. This approach has several advantages, as it has a low cost, does not use toxic reagents, and can be used when many soil samples are required, for example in precision agriculture. In addition, the possibility of using a smaller number of wavelengths to estimate soil chemical  Ca 2+ , Mg 2+ , and Al 3+ contents were extracted using KCl 1.0 mol L −1 ; K + and P contents were extracted using Mehlich-1; pH CaCl 2 (0.01 mol L −1 ) at a soil:solution rate of 1:2.5; H+Al was extracted using the SMP method; CECe: effective cation exchange capacity; CEC: cation exchange capacity at pH 7; m: aluminum saturation; BS: base saturation; R 2 : coefficient of determination; RMSE: root mean square error.