SOIL PROPERTIES MAPPING USING PROXIMAL AND REMOTE SENSING AS COVARIATE

ABSTRACT

Obtaining knowledge about the distribution of spatial variability of soil properties is crucial to the proper site-specific management. One way to improve the quality of soil mapping is by using auxiliary information (covariate). The objective of this study was to test whether remote and proximal sensing data can assist in soil properties mapping through geostatistical prediction. We worked in an experimental area cultivated with sugarcane located in Sao Paulo State, Brazil, and selected five soil properties: organic matter, CEC, base saturation, K and P availability. Two covariates often used to express soil variation were chosen, one obtained by remote sensing (SWIR2 band) and the other by proximal sensing (apparent soil electrical conductivity – ECa). These covariates were individually and together used in geostatistical interpolation method (kriging with external drift). We found that ECa is a more promising covariate than SWIR2 band from orbital imaging. Such proximal sensing can identify the soil short-range spatial variability. However, when the soil property variability is well explained by the sampling procedure, multivariate geostatistical methods may not improve the mapping accuracy.

KEYWORDS
satellite imagery; apparent soil electrical conductivity; kriging with external drift; site-specific management

INTRODUCTION

Precision agriculture (PA) identifies and treats the spatial and temporal variability of crops. To satisfactorily characterize such spatial variation of the soil, dense sampling grids are needed (Nanni et al., 2011). However, this type of sampling is costly and time-consuming, which often makes economically feasible sampling not sufficiently detailed to allow the proper identification of spatial variability of soil properties. Thus, the economic return of site-specific management is compromised.

Covariates can be used to assist in spatial predictions, without increasing the costs with soil sampling and laboratory analysis. Such auxiliary data are usually obtained at a greater spatial resolution and is easier and cheaper to obtain compared to manual samplings. The sole requirement is that the covariates must present empirical spatial correlation with the properties of the soil to be mapped. Among these possible covariates are the traditional remote sensing products, i.e. satellite imagery (Mirzaee et al., 2016; Sayão et al., 2018), and the most recent proximal sensing products. The first has the advantage of easy acquisition of data; the last have been developed for PA usage, which are sensors conducted near or in contact with the crop or soil (Adamchuk et al., 2004). The proximal sensing approach can avoid problems related to cloud cover that might appear in remote sensing data. One of the most used proximal sensing techniques to quantify the variability of agricultural soils is the apparent electrical conductivity (ECa), since it presents a direct relation with texture, water quantity and organic matter variations within the soil (Kühn et al., 2009).

Although the literature reports the potential of both remote and proximal sensing in assisting predictions of some soil properties, previous studies usually focus only on a single soil property that is highly correlated with the specific covariate under study (García-Tomillo et al., 2017; Sanches et al., 2018; Sayão et al., 2018), often conveying the idea of a potential that is not always widely applicable in other situations. Goovaerts (2000) argued that multivariate geostatistical methods (MGM) provide better predictions than ordinary kriging (OK) when working with correlations above 0.75. However, such high correlation values are difficult to find in PA data, due to randomness on soil spatial variability caused by anthropic interference. Nevertheless, Goovaerts & Kerry (2010), working with PA data and Kriging with External Drift (KED), argued that a correlation between 0.25 and 0.5 may be sufficient to improve prediction accuracy. Nevertheless, Sanches et al. (2018) improved chemical soil properties mapping through KED predictions working with lower relationships between soil properties and ECa data, ranging from 0.1 to 0.3. However, it remains unclear what is the minimum relationship level required between soil properties and covariates to significant prediction gains when compared to univariate interpolation methods, such as ordinary kriging. Thus, more research is needed to clarify how data relationship influences the quality of soil predictions through MGM.

Thus, our objective was to investigate whether remote sensing and proximal sensing could improve the quality of prediction of soil properties. Therefore, we selected five soil properties with different relationship levels with the covariates to evaluate the performance of spatial prediction using KED algorithm. Moreover, we used Bivariate Moran’s index (LISA) to express the relationship between soil sampling and covariates.

MATERIAL AND METHODS

Study site and soil data

This study used a 116 ha field (Figure 1) cultivated with sugarcane in the southeast of Brazil (21°30’29” S, 48°09’04” W), with a predominance of Oxisols. The climate of the region is Cwa – humid subtropical with dry winter and hot summer. Sugarcane is a semi-perennial crop, which means that with only one planting, the same plants (ratoons) may be generally harvested for three to seven cropping seasons, depending on soil quality, management practices, weather, and plant genotype. Generally, fertilizers are applied at the beginning of the crop cycle, immediately after the previous ratoon harvest.

FIGURE 1
Location of the study site and representation of the sampling points, indicating the soil samples.

We performed grid soil sampling in November 2011, after the final ratoon was harvested (September 2011), but before any kind of soil management practice was done for field reform. One hundred and sixteen soil samples were collected within the area, resulting in a density of approximately one sample per hectare. The soil was sampled in the 0-0.4 m depth using an auger. Each sample consisted of six subsamples collected within a 5 m radius to ensure proper representation of the sampling point since soil chemical variations may occur at very short distances depending on the position of the samples to the crop rows (1.5 m spacing between rows).

The following soil properties were chosen for this study: soil organic matter (OM, g kg⁻¹), cation exchange capacity (CEC, mmol_cdm⁻³), base saturation (V, %), available potassium (K, mmol_c dm⁻³) and available phosphorus (P, mg dm⁻³). Such choice was due to our research group focus on the improvement of variable-rate fertilization maps required for precision agriculture approach; OM is related to soil texture and can be used in some regions of Brazil to guide nitrogen prescription, whereas CEC and V may be used in liming prescription, while available P and K are generally used in the prescription of phosphate- and potassium-based fertilizers. Soil analyses were conducted in a commercial laboratory, using the standard analytical methods for São Paulo State, Brazil (Camargo et al., 2009). The phosphorus (P) availability was estimated via the extraction method with ion exchange resin, and potassium (K) availability was estimated via ammonia acetate 1N pH 7.0 (Camargo et al., 2009). OM was estimated via extraction (Camargo et al., 2009), while CEC was obtained through a compulsive exchange in BaCl2 (Camargo et al., 2009). V was estimated through the ratio between the sum of bases (K, Ca and Mg) and CEC in a 7 pH (Camargo et al., 2009).

Covariates

We employed two covariates that may provide information about soil variability, one from remote sensing and other from proximal sensing.

Apparent soil electrical conductivity (ECa) was the employed proximal sensing technique. This data was acquired in January 2012 using the tractor-driven Veris-3100 Soil Electrical Conductivity Sensor (Veris®, Salina, KS, USA). Data were collected before soil preparation for sugarcane planting in 30 m passes (Figure 2A). The ratoons were mechanically destroyed, and, subsequently, the soil was subsoiled and harrowed 15 days before ECa readings were collected. Since the objective was to estimate soil properties in the shallow soil depth (0-0.4 m), only the most superficial measuring depth of the sensor was used (up to 0.3 m deep). Discrepant data were removed to avoid negative interference in interpolations. For such, the histogram was created with equal classes and analyzed to search for discrepant values; when these were not inserted in a region with similar values (extreme values), they were removed (∼5% of the data were excluded as outliers). Further, we did a dataset reduction of ECa. Hence, we eliminated two thirds of the entire number of samples to improve computation performance for KED predictions. We followed the elimination using the direction of the sensor pass throughout the field, selecting three samples and excluding two of them. The data were interpolated to 10 m resolution using ordinary kriging, since KED demands covariate values in each node of the intended surface.

FIGURE 2
Apparent soil electrical conductivity data (spatial density from ∼272 samples ha⁻¹) (a) and bare soil reflectance in the SWIR2 band of a composition image of a historical series between 1984 and 2020 years (30 m spatial resolution) (b).

The remote sensing covariate is a shortwave infrared band 2 (SWIR2 - Figure 2B). The employed image was composed by a time-series image composition to guarantee a high-quality bare soil scene. Scenes from 1984 to 2020 (n = 432) were accessed from the Landsat Collection 2 – Tier 1 (Landsat 5, 7, and 8, comprising the respective spectral bands of 2.08-2.35, 2.09-2.35, and 2.10-2.30 μm), which are atmospherically corrected to surface reflectance by Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm, and geometrically corrected using ground control points and elevation data provided by a Digital Elevation Model (DEM). Then, following the procedure recommended by Demattê et al. (2020), we removed clouds and shadows using the quality assessment band (pixel_qa) and masked vegetation and straw pixel values to get pure bare soil pixels, resulting on an average of 42 reflectance value for each 30m pixel. For spatial predictions via KED, the time-series was reduced to a single SWIR2 image by applying the median over the entire spectral dataset, which was further resampled to 10 m resolution using bilinear method. This band was chosen due to its known correlation with the content of organic matter, clay and iron oxides in the soil (Sayão et al., 2018). Thus, it has the potential to indicate soil variation, although there is a lack of explored application in precision agriculture.

Data relationship

Bivariate Moran index (LISA) (Chen, 2015) was applied to evaluate the spatial relationship between soil properties and covariates (Figure 3). A Global LISA was applied to identify spatial correlation between variables. Yet, we also employed a Local LISA (Anselin, 1995) to identify the correlations for each observation (sample). Local LISA offers a notion of locals that present significant correlation in addition to the presence of high values of soil properties near to high values of covariates (high-high class); high values of soil properties near to low values of covariates (high-low class); low values of soil properties near to high values of covariates (low-high class); and low values of soil properties near to low values of covariates (low-low class).

FIGURE 3
Local and Global bivariate Moran index between soil properties and covariates. OM - organic matter; CEC – cation exchange capacity; V – base saturation; K – potassium; P – phosphorus; ECa – soil apparent electrical conductivity; SWIR2 – shortwave infrared band. ^ns – not significant correlation (p> 0.05), ^*p-value < 0.05, ^**p-value < 0.01, ^***p-value < 0.001.

Geostatistical interpolation

In this study, we used kriging with external drift (KED) as the multivariate geostatistical method (MGM), and in KED, the trend is estimated with a regression method (Hengl, 2009). Moreover, KED and Regression Kriging has mathematical similarities, delivering the same results when the stochastic component is modeled with a linear function (Hengl, 2009). We compared such MGM predictions to a univariate geostatistical method (UGM) to evaluate the prediction performance improvement when the covariates were used. Thus, Ordinary kriging (OK) or Universal Kriging (UK) were used as UGM, depending on whether the data shows a trend according to the coordinates.

Thus, three scenarios were built according to the covariates used: 1) “KED (ECa)”, where the drift was modeled as a function of ECa data; 2) “KED (SWIR2)”, where the drift was modeled as a function of SWIR band; and 3) “KED (Eca + SWIR2)”, where the drift was a combination of both covariates.

We employed the Restricted Maximum Likelihood (REML) to model the theoretical variogram. For such, two covariance models were tested: spherical and exponential. The one with the best performance measured by the leave-one-out cross-validation was chosen. We used the gstat package (Pebesma, 2004) from R for this geostatistical modeling.

Evaluation of prediction performance

The interpolation performances of the three KED scenarios and the reference method (UGM) were evaluated through leave-one-out cross-validation (LOOCV). Thus, the following prediction quality statistics were computed using observed and predicted values: mean error (ME), root-mean-squared error (RMSE) and Nash-Sutcliffe efficiency (NSE - Krause et al., 2005).

RESULTS AND DISCUSSION

We found an improvement in KED performance for mapping CEC, V, and K when we used ECa as covariate (Table 1). The improvement in prediction performance can be explained by the fact that ECa helped to explain data variance, which is corroborated by the reduction of the sill in the semivariograms (sill = c0 + c) mainly for CEC and V (Figure 4). The prediction of these two soil properties presented the highest improvement in performance compared to UGM (Table 1 – NSE). As K presented less reduction in the sill value resulted in a lower improvement in performance to UGM. On the other hand, KED did not bring any additional gain in OM and P prediction compared to UGM (Table 1); yet there was almost no reduction in the sill of KED (ECa). Thus, we understand that improvement of performance occurs when ECa helps to explain data variance and to reduce the sill, because it means the covariate can capture the deterministic part of soil properties spatial variability.

TABLE 1
Semivariogram parameters (c - partial sill; c_o - nugget; a – range; and theoretical model) and statistics of the prediction performance in the leave-one-out cross validation (LOOCV) of soil properties for the different geostatistical methods (UMG and KED scenarios), using ECa and SWIR2 as covariates.

Soil property	Geostatistical method	Semivariogram parameters				LOOCV
		c	co	a	Model**	EM	RMSE	NSE
OM (g kg−1)	UGM *	5.60	1.08	469.15	sph	0.00	1.76	0.84
	KED(ECa)	5.28	1.09	467.08	sph	0.00	1.75	0.84
	KED(SWIR2)	6.05	0.97	470.52	sph	−0.01	1.77	0.84
	KED(ECa + SWIR2)	5.72	0.99	468.16	sph	-0.01	1.77	0.84
CEC (mmolc dm−3)	UGM *	225.15	8.89	931.87	sph	0.09	6.39	0.85
	KED(ECa)	26.95	0.11	137.70	sph	−0.01	5.01	0.91
	KED(SWIR2)	107.95	7.77	428.32	sph	0.09	6.44	0.85
	KED(ECa + SWIR2)	27.11	0.18	138.50	sph	-0.01	5.03	0.91
V (%)	UGM	148.53	27.48	899.35	sph	0.07	7.19	0.52
	KED(ECa)	28.05	30.35	456.73	sph	0.00	6.57	0.60
	KED(SWIR2)	142.56	28.49	899.37	sph	0.08	7.21	0.51
	KED(ECa + SWIR2)	24.34	31.32	432.78	sph	-0.01	6.64	0.59
K (mmolc dm−3)	UGM *	0.30	0.03	205.10	exp	0.00	0.48	0.62
	KED(ECa)	0.14	0.09	385.44	sph	0.00	0.47	0.64
	KED(SWIR2)	0.21	0.05	319.90	sph	0.00	0.49	0.61
	KED(ECa + SWIR2)	0.15	0.09	385.37	sph	0.00	0.48	0.63
P (mg dm−3)	UGM	0.03	0.04	271.56	sph	0.04	5.71	0.02
	KED(ECa)	0.03	0.04	282.31	sph	0.04	5.74	0.01
	KED(SWIR2)	0.03	0.04	277.95	sph	0.04	5.75	0.01
	KED(ECa + SWIR2)	0.04	0.04	290.47	sph	0.04	5.80	−0.01

• OM - organic matter; CEC – cation exchange capacity; V – base saturation, K – available potassium; P – available phosphorus; ECa – soil apparent electrical conductivity; SWIR2 – shortwave infrared band. EM – mean error, RMSE – root-mean squared error, NSE – amount of variance explained.
• ^*
  when trend was modeled by coordinates, we used universal kriging as UGM method.
• ^**
  Semivariogram model: Sph – Spherical; Exp – Exponential.

FIGURE 4
Semivariograms for all five soil properties (principal variables) according the geostatistical method: UGM - Ordinary Kriging; KED(ECa) – Residual Semivariograms when Kriging with External Drift was used with apparent soil electrical conductivity as auxiliary variable; KED(SWIR2) - Residual Semivariograms when Kriging with external Drift was used with the Short Infrared band SWIR2 as auxiliary variable; KED(ECa+SWIR2) - Residual Semivariograms when Kriging with external Drift was used with both covariates. OM – soil organic matter; CEC - cation exchange capacity; V - base saturation; P – available phosphorus, K - available potassium. The “circle” symbol is the experimental semivariogram and the solid lines are the theoretical semivariograms adapted by residual maximum likelihood (REML). NOTE: KED_ECa+swir2 are overlapping the KED_ECa.

The use of soil ECa showed a positive and complementary effect on the spatial prediction process. This happens because the insertion of covariates makes the residuals spatially stationary and normally distributed (Wu et al., 2019). Thus, ECa as a covariate provided information gains, especially when there was a high correlation with the soil property (i.g. Global Moran Index of 0.76 to CEC and 0.57 to V variables – Figure 3). This may be because this sensing technique is influenced not only by the most superficial soil layer– like SWIR2 –, but its signal represents the variation of soil up to 0.3 m deep, closer to the soil sampling depth. In addition, ECa measures a greater volume of soil which can assist in identifying the short-range spatial variability of the soil. This result in greater spatial correlations with some soil properties when compared to SWIR2 (Figure 3), thus leading to better prediction performance. Moreover, the use of ECa seems more promising given the widespread adoption of the no-tillage farming system, where the traditional soil management processes – i.e., horizontal preparation of the soil by plows and harrows – are decreasing. Thus, there will be fewer opportunities to obtain orbital imagery in the entirely bare soil situation, hindering auxiliary data acquisition on periods close to soil sampling.

Unlike ECa, remote sensing data (SWIR2 band) did not show good suitability as covariate in predicting the soil properties analyzed in this study (Table 1). In general, predictions using SWIR2 showed equivalent or even inferior performance to the control method used (UGM). However, Sorenson et al., (2021), report that the SWIR2 band from bare soil was important for CTC prediction, leading to the inference that CTC responds to relative changes in clay content. Moreover, KED scenario where we used both remote and proximal data did not yield better predictions than where ECa were solely used in KED (Table 1). Such an unsatisfactory result of using remote sensing data may have occurred due to some reasons: 1) this technique only measures the reflectance of the soil surface, which can impair the relationships with soil samples collected deeper; 2) there may be small scattered portions of straw or weeds on the soil surface that were not filtered by Demattê et al. (2020) technique employed, which influenced the reflectance of the pixels; 3) we used a mosaic of images of 36 years, and since it is an agricultural field, the combination of images could impair the relationship of SWIR2 with soil properties. The main advantage of using orbital imagery for this purpose is the lower acquisition cost, dispensing additional fieldwork, once there is a significate amount of freely-available databases. However, imagery acquisition of bare soil has some difficulties because of the presence of clouds and the adoption of no-tillage systems as well as integrated crop-livestock systems, which are more soil-use-intensive. This makes it even more challenging to access bare soil data through remote sensing.

In the case of the gathering of soil ECa, there are several commercial sensors available since this technology has been more incorporated in precision agriculture practice. Because this measure gives valuable information about soil variability, it has been employed in research for soil sampling optimization (Sanches et al., 2020) and delimitation of management zones (Betzek et al., 2019). Furthermore, the variations of ECa over the area is relatively stable over time, once mineralogy and texture do not usually change over the crops’ seasons. Thus, we should expect better results for usage ECa as a covariate for mapping soil properties that are stable over time, like soil texture and CEC. On the other hand, variables that are easily changed by management, like P and K, tend to have a lower gain when using methods combined with ECa data, corroborating the results we found.

The spatial correlation between soil properties and covariates is an essential source of information to evaluate the possibility of working with multivariate geostatistical methods. In previous studies, Sanches et al. (2018) and Jurado-Expósito et al. (2019) argued that statistically significant relationships between the primary and covariates improve the predictions of soil properties. However, in our study, we did not observe improvements in OM estimation, which despite having presented high and statistically significant correlations (Figure 3) with the covariates did not provide improvements in prediction accuracy (Table 1-RMSE). This probably happened because soil OM was sufficiently explained by the sample density, corroborated by a good variogram fitting (Figure 4). However, for the K and V properties that presented intermediate but statistically significant correlations (Figure 3) and a lower spatial relationship than OM, the results were divergent. V predictions using covariates showed a more significant increase in prediction performance than for K (table 1 - NSE). These properties (K and V) showed trend behavior in the semivariogram, which means that other factors are responsible for the variation were found for the analyzed property. The modeling with covariates helps in this sense, modeling trend and more clearly displaying the effect of micro variation that is, variations that occur at distances smaller than the sampled. By recognizing this effect of trend, we have smaller errors by adding covariates in predicting soil attributes. Therefore, when the sample density used is sufficient to capture the variability of properties, the addition of covariates, although presenting high correlations, does not contribute to the modeling of variability and, consequently, in improving predictions.

Although P was not correlated with the auxiliary variables used (Figure 3), all analyses were performed to compare the results. A global LISA of 0 indicates that the relationship of P with the covariate occurs randomly. As expected, we did not find gains in the prediction quality when modeling the variability of this property with covariates. Therefore, the KED scenarios presented results similar to UGM, all with low quality in predictions (NSE < 0.02 - Table 1). Furthermore, as the NSE value was close to zero (Table 1) for P, this indicates that the best way to manage phosphate fertilizers would be with uniform application throughout the area.

The spatial distribution of CEC and V (maps on Figure 5) showed similarities to the ECa distribution. However, in the other soil properties, where the KED (ECa) did not improve (OM, P) the predictions, the spatial distribution of the ECa did not clearly appear in the maps. Otherwise, we observed the representation of some high ECa value spots in the spatial distribution of K (Figure 2 compared to Figure 5). Thus, as indicated by local LISA (Figure 3), the sites with high K contents coincide with the sites with higher ECa. Thus, when we insert covariates that even marginally explains the soil variability, we may obtain a better spatial distribution characterization, which may lead to better fertilizer prescriptions in variable-rates.

FIGURE 5
Spatial distribution of soil properties for the different geostatistical methods: UGM - Ordinary Kriging; KED(ECa) – Kriging with External Drift using apparent soil electrical conductivity as covariate; KED(SWIR2) - Kriging with external Drift with the Short Infrared band (SWIR2) as covariate; KED(ECa+SWIR2) - Kriging with external Drift with both covariates. OM - organic matter; CEC- cation exchange capacity; V - base saturation; P – available phosphorus; K – available potassium.

CONCLUSIONS

Our goal was to investigate whether remote sensing and proximal sensing can improve the quality of prediction of soil properties through KED. Therefore, we selected five soil properties to evaluate the performance of spatial prediction when adding different levels of correlation between soil property and covariates. Thus, we concluded that ECa is a more promising covariate than SWIR2 band from orbital imaging. Moreover, ECa data collection is more flexible, not depending on bare soil image collection. In addition, ECa data is densely collected, which can assist in identifying the short-range spatial variability of the soil. However, when the soil property variability is well explained by the sampling procedure, multivariate geostatistical methods (MGM - like KED) may not improve the mapping accuracy, even when the spatial relationship between the variables is high.

ACKNOWLEDGEMENTS

We would like to thank professor José Paulo Molin, from University of São Paulo, for the financial support for data collection. This study was possible thanks to Ph.D scholarship granted by Coordination for the Improvement of Higher Education Personnel – Brazil (CAPES) – Financing Code 001, to the first author, and the master scholarship granted by São Paulo Research Foundation (FAPESP), process n° 2018/25473-2.

REFERENCES

• Adamchuk VI, Hummel JW, Morgan MT, Upadhyaya SK (2004) On-the-go soil sensors for precision agriculture. Computers and Electronics in Agriculture 44: 71–91. DOI: http://dx.doi.org/10.1016/j.compag.2004.03.002
  » http://dx.doi.org/10.1016/j.compag.2004.03.002
• Anselin L (1995) Local Indicators of Spatial Association - LISA. Geographical analysis 27: 93–115. DOI: http://dx.doi.org/10.1111/j.1538-4632.1995.tb00338.x
  » http://dx.doi.org/10.1111/j.1538-4632.1995.tb00338.x
• Betzek NM, Souza EG de, Bazzi CL, Schenatto K, Gavioli A, Magalhães PSG (2019) Computational routines for the automatic selection of the best parameters used by interpolation methods to create thematic maps. Computers and Electronics in Agriculture 157: 49–62. DOI: http://dx.doi.org/10.1016/j.compag.2018.12.004
  » http://dx.doi.org/10.1016/j.compag.2018.12.004
• Camargo OA, Moniz AC, Jorge JA, Valadares JMAS (2009) Methods of chemical, mineralogical and physical analysis of soils of the Agronomic Institute of Campinas. Campinas, Instituto Agronômico, 2 ed. 94p. (Boletim Técnico 106)
• Chen Y (2015) A New Methodology of Spatial Cross-Correlation Analysis. Plos ONE 10: 1–20. DOI: http://dx.doi.org/10.1371/journal.pone.0126158
  » http://dx.doi.org/10.1371/journal.pone.0126158
• Demattê JAM, Safanelli JL, Poppiel RR, Rizzo R, Silvero NEQ, Mendes WS, Bonfatti BR, Dotto AC, Salazar DFU, Mello FAO, Paiva AFS, Souza AB, Santos NV dos, Nascimento CM, Mello DC, Bellinaso H, Gonzaga Neto L, Amorim MTA, Rezende MEB de, Vieira JS, Queiroz LG de, Gallo BC, Sayão VM, Lisboa CJS (2020) Bare earth’s surface spectra as a proxy for soil resource monitoring. Scientific Reports 10: 4461. DOI: http://dx.doi.org/10.1038/s41598-020-61408-1
  » http://dx.doi.org/10.1038/s41598-020-61408-1
• García-Tomillo A, Mirás-Avalos JM, Dafonte-Dafonte J, Paz-González A (2017) Estimating soil organic matter using interpolation methods with an electromagnetic induction sensor and topographic parameters: a case study in a humid region. Precision Agriculture 18: 882–897. DOI: http://dx.doi.org/10.1007/s11119-016-9481-6
  » http://dx.doi.org/10.1007/s11119-016-9481-6
• Goovaerts P (2000) Geostatistical approaches for incorporating elevation into the spatial interpolation of rainfall. Journal of Hydrology 228: 113-129. DOI: http://dx.doi.org/10.1016/S0022-1694(00)00144-X
  » http://dx.doi.org/10.1016/S0022-1694(00)00144-X
• Goovaerts P, Kerry R (2010) Using ancillary data to improve prediction of soil and crop attributes in precision agriculture. In: Oliver M. Geostatistical Applications for Precision Agriculture. Dordrecht, Springer. DOI: http://dx.doi.org/10.1007/978-90-481-9133-8
  » http://dx.doi.org/10.1007/978-90-481-9133-8
• Hengl T (2009) A practical guide to the geostatistical mapping. Amsterdam, 266p.
• Jurado-Expósito M, Castro AI de, Torres-Sánchez J, Jiménez-Brenes FM, López-Granados F (2019) Papaver rhoeas L. mapping with cokriging using UAV imagery. Precision Agriculture 20: 1045-1067. DOI: http://dx.doi.org/10.1007/s11119-019-09635-z
  » http://dx.doi.org/10.1007/s11119-019-09635-z
• Krause P, Boyle DP, Bäse F (2005) Comparison of different efficiency criteria for hydrological model assessment. Advances in Geosciences 5: 89–97. DOI: http://dx.doi.org/10.5194/adgeo-5-89-2005
  » http://dx.doi.org/10.5194/adgeo-5-89-2005
• Kühn J, Brenning A, Wehrhan M, Koszinski S, Sommer M (2009) Interpretation of electrical conductivity patterns by soil properties and geological maps for precision agriculture. Precision Agriculture 10: 490–507. DOI: http://dx.doi.org/10.1007/s11119-008-9103-z
  » http://dx.doi.org/10.1007/s11119-008-9103-z
• Sanches GM, Magalhães PSG, Luciano ACS, Camargo LA, Franco HCJ (2020) Comprehensive assessment of spatial soil variability related to topographic parameters in sugarcane fields. Geoderma 362: 114012. DOI: http://dx.doi.org/10.1016/j.geoderma.2019.114012
  » http://dx.doi.org/10.1016/j.geoderma.2019.114012
• Sanches GM, Magalhães PSG, Remacre AZ, Franco HCJ, (2018) Potential of apparent soil electrical conductivity to describe the soil pH and improve lime application in a clayey soil. Soil and Tillage Research 175: 217–225. DOI: http://dx.doi.org/10.1016/j.still.2017.09.010
  » http://dx.doi.org/10.1016/j.still.2017.09.010
• Mirzaee S, Ghorbani-Dashtaki S, Mohammadi J, Asadi H, Asadzadeh F (2016) Spatial variability of soil organic matter using remote sensing data. Catena 145:118–127. DOI: http://dx.doi.org/10.1016/j.catena.2016.05.023
  » http://dx.doi.org/10.1016/j.catena.2016.05.023
• Nanni MR, Povh FP, Demattê JAM, Oliveira RB, Chicati ML, Cezar E (2011) Optimum size in grid soil sampling for variable rate application in site-specific management. Scientia Agricola 68: 386–392. DOI: http://dx.doi.org/10.1590/S0103-90162011000300017
  » http://dx.doi.org/10.1590/S0103-90162011000300017
• Pebesma EJ (2004) Multivariable geostatistics in S: the gstat package. Computers & Geosciences 30: 683–691. DOI: http://dx.doi.org/10.1016/j.cageo.2004.03.012
  » http://dx.doi.org/10.1016/j.cageo.2004.03.012
• Sayão VM, Demattê JAM, Bedin LG, Nanni MR, Rizzo R (2018) Satellite land surface temperature and reflectance related with soil attributes. Geoderma 325: 125–140. DOI: http://dx.doi.org/10.1016/j.geoderma.2018.03.026
  » http://dx.doi.org/10.1016/j.geoderma.2018.03.026
• Sorenson PT, Shirtliffe SJ, Bedard-Haughn AK (2021) Predictive soil mapping using historic bare soil composite imagery and legacy soil survey data. Geoderma 401: 115316. DOI: http://dx.doi.org/10.1016/j.geoderma.2021.115316
  » http://dx.doi.org/10.1016/j.geoderma.2021.115316
• Wu Z, Wang B, Huang J, An Z, Jiang P, Chen Y, Liu Y (2019) Estimating soil organic carbon density in plains using landscape metric-based regression kriging model. Soil and Tillage Research 195: 104381. DOI: http://dx.doi.org/10.1016/j.still.2019.104381
  » http://dx.doi.org/10.1016/j.still.2019.104381

Edited by

Area Editor: Teresa Cristina Tarlé Pissarra

Publication Dates

• Publication in this collection
  17 Dec 2021
• Date of issue
  Nov-Dec 2021

History

• Received
  05 July 2021
• Accepted
  29 Oct 2021