The proposal of precision agriculture (PA) aimed at localized soil management involves agronomic, technological and, consequently, financial changes, which hinders adoption of the PA techniques, especially by small producers (Srinivasan, 2006). PA allows studying and analysing the spatial variability of crop yield and the chemical and physical attributes of the area under study. The study of this spatial variability of spatially georeferenced continuous variables can be done using geostatistical techniques (Cressie, 2015), which is crucial for the feasibility of differentiated and localized management practices, and allowing to measure the degree of spatial dependence between the sample elements in a given region and to describe the spatial dependence structure of the georeferenced variable throughout the area, thus generating thematic maps (Cressie, 2015).

An important part of the PA is to obtain data for the elaboration of maps of attributes that influence crop productivity (Bernardi et al., 2014). There are commercially available equipment and sensors that can be highlighted for obtaining these data, as well as for collecting soil samples (Gonçalves et al., 2020).

Agricultural sensors are instruments that can be attached to implements or agricultural machinery, allowing, for example, to analyze soil attributes or characteristics in real time. Scherer et al. (2018) proposed the development of a sensor to determine the electrical conductivity of the soil at low cost. Prudente et al. (2021) compared the NDVI (Normalized Difference Vegetation Index) spectro-temporal profiles obtained by active (GreenSeeker 505 Handheld) and passive (FieldSpec4 Standard-Res model) proximal sensors to monitor soybeans and beans. By means of multispectral sensors onboard the AT120 Remotely Piloted Aircraft Systems, Facco & Pegoraro (2019) captured aerial images, and thus built digital surface models, orthoimages, contour lines, generated vegetation indices and performed the detection of planting failures. Thinking about agility and quality, Resende et al. (2020) used a remotely piloted aircraft equipped with an RGB camera and a MAPPIR 3 camera to capture images of a corn crop, to estimate the leaf area index of a plot infested by Spodoptera frugiperda. Furthermore, the wireless sensor network can be used in underground and terrestrial environments to detect soil and climate conditions, as well as detect pests and insects (Bayrakdar, 2020a,b).

Another important procedure in PA is soil sampling, which represents a critical step in the assessment of soil fertility, and may be responsible for 98% of the errors made in the inappropriate recommendation of fertilizers (Mendes et al., 2001). Still, the choice and influence of the location of the sampled points, as well as the distance that separates them, are essential for the success of the sampling design (Carvalho & Nicollela, 2002).

In this context, it is necessary to determine a soil sampling scheme that is as efficient as possible, with the smallest sample size to minimize operating costs and maximize the quality of the results of the spatial variability analysis (Marchant & Lark, 2010).

For the choice of this sampling scheme, there are traditional spatial sampling designs, such as simple random sampling and the systematic design, which are classified as design-based, assuming that the values of the georeferenced variable are considered fixed. Traditional spatial sampling designs are advantageous in terms of the simplicity of choice of sampling points in the area under study, either because of the occurrence of different spacing between pairs of sampling points (in simple random sampling) or because of coverage uniformity of the area under study (in systematic sampling) (Haining, 2015).

However, the spatial heterogeneity and spatial autocorrelation perform a particular impact on the variability analysis accuracy, in the context of the sampling design. The spatial distribution of an attribute is considered heterogeneous when the attribute varies in the study area, thus requiring division of the area into homogeneous subareas. On the other hand, spatial autocorrelation indicates that the georeferenced samples nearby have similar values.

Characteristics of spatial variability, such as anisotropy, can guide a choice of the most intense number of sampling points in the direction of greater spatial variability (Delmelle, 2009). Or the inclusion of more sampling points with a smaller distance radius in the systematic sampling design, to minimize the nugget effect problem, as proposed by Diggle & Lophaven (2006), in the “lattice plus close pairs” and “lattice plus infill” designs.

The spatial sampling designs that consider the sample as a realization of a probability process are classified as model-based design or geostatistical design (Diggle & Lophaven, 2006). The spatial sampling design that considers the spatial correlation makes it possible to reduce the sample size by restricting redundant spatial information, in relation to traditional sampling (Haining, 2015; Dal’Canton et al., 2021).

Several methodologies for spatial re-planning of the sample, with sample size reduction, are developed in the literature, such as calculation of the effective sample size (Dal’Canton et al., 2021) or the optimized sample re-planning (Maltauro et al., 2019; 2021). Moreover, prior knowledge of the spatial autocorrelation of the variable can be used in the spatial stratified sampling, which provides a lesser loss of precision of the estimate than spatial simple random and systematic sampling. The spatial stratified sampling can use relevant previous historical data in the same area, such as prior knowledge. Associated with the multivariate clustering analysis, this information allows generating zones, which will represent the stratified space (Wang et al., 2013).

In the PA context, the identification of homogeneous Management Zones (MZs) within the cultivation areas becomes a more technically and economically viable strategy for the implementation of localized soil management, allowing a long-term that all area treatments to be carried out uniformly within each zone (Molin, 2006; Srinivasan, 2006; Bernardi et al., 2014). Thus, to produce more stable MZs, it is recommended to use variables that do not vary significantly over time, such as topographic data and physical data (Aikes et al., 2021).

On the other hand, the chemical attributes are important and useful to generate the zones for recommendations of variable rate fertilizer applications (Aikes et al., 2021). These zones are called application zones (AZs) (Molin, 2006). Therefore, the difference between MZs and AZs is related to the variables available (stables or not) and with the intention of generating the zones (long-term use or only for a future application of fertilizers). Moreover, the spatial statistics and the multivariate cluster analysis are methods commonly used by MZ and AZ.

These MZs or AZs make up a database that can be used to direct, in the future, the reduction of the sample size in areas more homogeneous (Rodrigues Jr. et al., 2011). In this context, there are several ways to select samples, which depend on the purpose and the resources available for the survey (Haining, 2015). When we do not have information about the spatial variability, the sampling points can be selected using traditional sampling criteria (Benedetti et al., 2015), or multi-phase adaptive sampling schemes (Marchant & Lark, 2010).

Determination of the sample configurations with reduced size by the use of AZs as a spatial stratified sampling provide the following advantages: AZs generate more homogeneous regions from the point of view of spatial similarity, allowing the researcher to choose a smaller number of sampling points in each AZ; the methodology considers the spatial information redundant; there is a design-based and model-based combination; use of less computational time than in optimized sampling; and greater simplicity of execution.

Given the above, the objective of the paper was to assess the acquisition of reduced sample configurations contained in application zones generated by spatial multivariate clustering. The sampling protocol proposed in this work is intended to guide the definition of new intervention cycles in the study area, with a reduced sample size, but maintained the refined localized management in this agricultural area.

Regarding all the soil’s chemical attributes for the 2013-2014 and 2014-2015 harvest years, there was a wide variation in the coefficient of variation (CV), from 22.41 to 60.18 and from 25.11 to 126.28 respectively, what indicates high dispersion (20 < CV ≤ 30) or heterogeneity of the soil chemical attributes (CV > 30). The 2015-2016 harvest year presented a CV between 10.58 and 40.92, ranging from average (10 ≤ CV ≤ 20) to high (20 < CV ≤ 30) dispersion of the soil chemical attributes (Table 1).

Soil chemical attributes Ca, C, Cu, Fe, H+Al³, Mn, and SB had mean values considered average or high, whereas the mean value of Zn can be classified as low or average, and that of Al can be classified as low regarding the need for the soil (Table 1). The Zn content of the 2014-2015 harvest year presented a moderate linear association of its respective values with the X-axis coordinates, that is, showing a mean linear trend of the deterministic term in relation to the East-West direction (Table 1) (r ≥ 0.30).

For all the soil’s chemical attributes of the 2013-2014 harvest year, an estimated value for the spatial dependence radius (range) from 157.70 m to 707.86 m was observed (Table 1). For the 2014-2015 harvest year, the practical range varied from 128.61 m to 453.07 m. The 2015-2016 harvest year presented greater variation in the practical range, between 149.73 m to 855.10 m. Regarding spatial dependence intensity, and according to the criteria used, it is observed that, for all the harvest years, the soil chemical attributes presented moderate or strong spatial dependence (Table 1).

Regardless of considering the same variables year to year to generate zones, Schenatto et al. (2016), working with different soil attributes, also found different amounts of the zones for an agricultural area considering different harvest years. The sampling points obtained by the reduced sampling configurations (points in red in Fig. 3) within each AZ (Fig. 3), showed that sampling points selected throughout the study area, seeking greater concentration of points in regions where there is greater variability, as well as a reduction in sample density in more uniform locations (Rodrigues Jr. et al., 2011).

The similarity in descriptive statistics when comparing the reduced sample configurations with the initial sample configurations was also found by Maltauro et al. (2019); and Dal’Canton et al. (2021) working with the sample reduction in the same agricultural area with grain cultivation, considering the chemical attributes of the soil.

The results regarding spatial dependence are explained by the bilateral relationship of sample size in the spatial scale (Kerry et al., 2010). A reduction in sample size can generate a reduction in spatial dependence intensity, as evidenced by these authors, who presented examples with sample size reduction that generated an increase in the nugget effect. These authors also mention the importance of defining an adequate sample size, to remove the effect of the “noise” that may have been generated exclusively by the choice of sample size.

In general, most of the soil’s chemical attributes showed moderate spatial dependence; this fact makes the thematic maps more accurate than those generated considering weak spatial dependence (Cambardella et al., 1994). In addition, a high nugget effect indicates low spatial dependence and leads to estimates around the sample mean (Alencar et al., 2019), as the nugget effect is associated with sampling or analysis errors, which indicates that two observations close to each other have very different values. Regarding the range, Dal’Canton et al. (2021) concluded that the larger the spatial dependence radius, the smaller the sample size, estimated by the effective size of the univariate sample and the greater the area of homogeneity between sampling points.

The best estimates for the accuracy indices were observed comparing the initial sample configurations with those reduced by 25% (Table 2). This result was already expected, for containing more sampling points than the 50% sample reduction. However, it was possible to observe that the spatial variability pattern was maintained in most classes of the thematic map of the soil chemical attributes, even considering the largest sample reduction. This same trend was observed in Maltauro et al. (2019) and by Dal’Canton et al. (2021), even the authors working with different methods for sample reduction in an area with grain cultivation.

For all the harvest years, the clustering methods were efficient for defining the AZs and, for the 2013-2014 and 2015-2016, the best number of clusters for all the clustering methods was kc=2. For the 2014-2015 harvest year and for most of the clustering methods, the ideal number of clusters was kc=3. Considering the evaluation criteria, K-means was the best clustering method. Therefore, from a practical point of view, it is concluded that the AZs allow for localized application of inputs in the agricultural area.

Among the configurations analysed, the sample configuration proportionally reduced by 25%, when compared to the initial sample configuration, presented the best estimates for the values of the spatial dependence radius and the highest values for the accuracy indices; while the worst estimates for the accuracy indices were shown when comparing the initial sample configuration and those randomly and systematically reduced by 50%.

Overall, the AZs allowed dividing the agricultural area into more homogeneous sub-regions, as well as to select a smaller number of points within each AZ. Thus, the definition of AZs is efficient in obtaining a reduced sample configuration and in defining future soil samplings in the study area, aiming to minimize, in the long term, the spatial variability of the soil chemical attributes soil in this agricultural area, also allowing the producer to reduce costs when carrying out the soil analysis. This methodology used presents a design-based and model-based combination, greater simplicity in execution and brief computational times.

As future works, we seek to apply the methodology developed in spatio-temporal data.