The Field data from an avocado farm of 1.1 ha of avocado (Persea americana) var Hass planted at high density, the distance between trees is 2 m and between rows 5 m was collected (Figure 1). The study area is located in the city of Barranca, a province located 190 km northwest of the city of Lima, Perú, in the coastal desert, in the foothills of the western slope of the central Andes (Lira-Camargo et al. 2020). Given this location, an almost total absence of precipitation was perceived, as well as a high level of atmospheric humidity and persistent cloudiness.

Figure 1

The average annual temperature at the study site is 17.5 to 19 °C, with an annual summer maximum of about 29 °C. Temperatures during the summer months, from December to April, range from 29 to 30 °C during the day and 21 to 22 °C at night (Ruiz et al. 2019).

The avocado trees have a height of 5.20 m, which were collected by the same expert. A total of 202 observations in which the measurement variable is the quantity of the crystalline mite Oligonychus sp. per leaf were collected with a Qfield (2019), which allows working on QGIS projects for field work equipment (QField 2019). The projection system with the European Petroleum Survey Group (EPSG) code 32718 - WGS 84/UTM zone 18S was used. Data manipulation and analysis was performed using R statistical software and QGIS software.

For the estimation of the risk model, a geostatistical approach was applied; in this, the random variable Z(x) with spatial index varies continuously across the spatial region D (Cressie 1992). The behavior of the frequency distribution and the possibility of a transformation were evaluated, as geostatistical analysis and predictions are sensitive to outliers (Giraldo 2002).

Applied to agronomic sciences, geostatistics considers each sample value Z(x) it is associated with a position and uses this same dependence to make inferences about the distribution of the data, in this approximation a model of spatial variation that contains at least three components: first a general structure, which can be defined as a trend; a second, superimposed structure, related to spatial correlation and gradual variation; and finally a third component consisting of random variation caused by sampling errors or spatial variations at different scales sampling errors or spatial variations at scales smaller than the sample network (Cressie 1992; Li and Heap 2008).

In other words, it should be evaluated whether the process generating the data is stationary by mean or formally expressed as follows E(Z(x_i)) = μ (Giraldo 2002), to test this situation, a first-order polynomial surface was estimated and the significance of the model and the adjusted R² were evaluated, since this indicator essentially imposes a penalty for including additional predictors that do not contribute much to explaining the variation observed in the variable response (Diez et al. 2020).

Secondly, it is necessary to test whether the process is stationary in variance , i.e., whether the data is isotropic or anisotropic (Giraldo 2002; Plant 2018), this be analyzed by defining the behavior of spatial dependence in cardinal directions (Bivand et al. 2008). Thirdly, the spatial dependence must be modeled using the experimental semivariogram, where is the semivariance for all samples located in space, separated by the distance interval h; N_(h) is the total number of pairs of samples separated by a distance interval h; Z_(x) is the value of the sample at a location x; Z_(x+h) is the value of the sample at interval distance h from x . Variogram modelling and estimation is important for structural analysis and spatial interpolation (Oliver and Webster 2015; Li and Heap 2008).

For the adjustment of the theoretical semivariogram, the method proposed by Bivand et al. (2008) was used, which utilizes non-linear regression to fit the coefficients. For this, a weighted sum of square errors with the value according to the parametric model is minimized. A theoretical Spherical semivariogram model was adjusted, following Equation (1).

The parameter C₀ represents the nugget value induced by the spatial error when the distance is smaller than the lag distance. Whereas the variance of the process is denoted by C₁ and a combination between C₀ and C₁, where C₀ is referred to as sill. The parameter C₂ is the range in which it indicates a distance from the origin to the point of sill achieved (Moonchai and Chutsagulprom 2020).

The parameters previously analyzed, anisotropy and spatial dependence, allowed us to determine the type of geostatistical predictor interpolator most suitable to be used, since the three large families of geostatistical interpolators differ in three basic conditions: 1) Simple Kriging: the process is stationary by mean and the population mean is known; 2) Ordinary Kriging: the process is stationary by mean and the population mean is unknown; 3) Universal and Residual Kriging: the process is not stationary by mean and the population mean is unknown (Bivand et al. 2008; Li and Heap 2008; Cressie 1992; Oliver and Webster 2015).

The process of validation and quality control of the interpolation is similar to that of traditional statistical models, evaluation of the behavior of the residual semivariogram, homoscedasticity, normality of residuals (Bivand et al. 2008; Cressie 1992; Gujarati and Porter 2010). Also, Root-mean-square deviation (RMSE) and Mean Error (ME) were estimated as a validation process, this is a standardized methodology used by Andrades-Grassi et al. (2020, 2021). In addition, it was estimated the z-score, computed as: with the cross validation prediction for xi, and σ_[i](x_i) the corresponding kriging standard error (Bivand et al. 2008).

Finally, an indicator function of risk was estimated by Oligonychus sp., for this the kriging predictor method selected according to the validation process was used. To convert the prediction into a risk function, the Indicator kriging was used, in this type of kriging we are interested in the event: at a given point x the value Z(x) exceeds the level Z₀. The event can be represented by a binary function, the indicator function, valued 1 if the event is true, and 0 if it is false, whose expected value is the probability of the event Z(x), exceeds Z₀ (Chilès and Delfiner 2009). In this case it must be assumed that the process of the regionalized variable is stationary in which the following transformation is defined: (Chilès and Delfiner 2009; Giraldo 2002). Values which are much greater than a given cut-off, will receive the same indicator value as those values which are only slightly greater than that cut-off (Giraldo 2002).

Thus, indicator transformation of data is an effective way of limiting the effect of very high values. Simple or Ordinary kriging of a set of indicator-transformed values will provide a resultant value between 0 and 1 for each point estimate. This is in effect an estimate of the proportion of the values in the neighborhood which are greater than the indicator or threshold value (Chilès and Delfiner 2009). In this case, the threshold was set between 50-100 motile forms of the Hass variety of Oligonychus sp. per leaf which is equivalent to 40% infestation, since leaf drop occurs when the percentage of damaged area is equal to or above 7.5-10 %, which corresponds to spider populations exceeding 100-500 individuals per leaf (Hoddle et al. 2002). As can be seen, this is a wide range, so the risk function cutoff was set at five thresholds (50, 60, 70, 80 and 100 count/leaf).

As it was evidenced, up to this moment exclusively univariate geostatistical models were evaluated, therefore, it is not known if there are mappable geomorphometric causal factors that have an influence on the growth of Oligonychus sp. mite, the Digital Elevation Model (DEM) of ALOS Parsal of 12.5 m was used (Racoviteanu et al. 2007), as this was the free product with the best spatial resolution.

It is known that the spectral response of vegetation changes according to whether it is affected or vigorous. The typical behavior of vigorous vegetation shows reduced reflectivity in the visible spectrum, but high reflectivity in the near infrared, gradually decreasing towards the mid-infrared.

These spectral characteristics are primarily related to the action of photosynthetic pigments and water stored in the leaves. Specifically, the low reflectivity in the visible portion of the spectrum is due to the absorbing effect of the leaf pigments. As for the high reflectivity in the near infrared, it is due to the internal cellular structure of the leaf. In particular, the spongy layer of the mesophyll, with its internal air cavities, plays a leading role by scattering and dispersing most of the incident radiation in this band of the spectrum (Harris 1987). Therefore, the healthy leaf offers high reflectivity in the near infrared (between 0.7 µm and 1.3 µm), in clear contrast to the low reflectivity it offers in the visible spectrum, especially in the red range (Crippen 1990).

There are different vegetation indices that could indicate if there is a potential problem in a crop and use broad information from the electromagnetic spectrum (Crippen 1990). The Normalized Difference Vegetation Index (NDVI) is defined as the parameter calculated from reflectance values at different wavelengths and is particularly sensitive to vegetation cover. The use of this index is based on the particular radiometric behavior of vegetation, which is characterized by the contrast between the red band (between 0.6 μm and 0.7 μm), which is largely absorbed by the leaves, and the near infrared (between 0.7 μm and 1.1 μm), which is mostly reflected. The NDVI is estimated as follows NDVI = (IRC - R) / (IRC + R), where IRC is the reflectivity of the near infrared band and R is the reflectivity of the red band.

The NDVI were generated from a Sentinel-2 image taken on October 10, 2020, at 10 m spatial (Copernicus 2021) resolution, this was preprocessed and transformed to reflectance using the Qgis SPC plug-in (Congedo 2021). It was as close as possible in time to the date of field spatial data collection of Oligonychus sp. Information derived from vegetation indices are not causal risk factors, but rather a consequence of the effect of mite damage. Low values of vegetation indices usually indicate low vigorous vegetation, while high values indicate very vigorous vegetation (Crippen 1990). As this DEM and NDVI do not have a spatial resolution and level of detail that is sufficient for geostatistical analysis, a multi-Gaussian conditional simulation process was performed. One of the main problems of Kriging predictions is that the smoothing of the map of estimated values is smoother than the map of true values (Chilès and Delfiner 2009).

A multi-Gaussian random function of variogram g(h) was simulated at the sites x₁,... x_n in the space. The sequential simulation is performed as follows: 1) simulate a Gaussian value U₁ (mean 0, variance 2) with Y(x₁)=U₁. 2) for each i {2,…n}, with: Y(x_i)=Y(x_i)^SK + σ_SK(x_i)U_i, where Y(x_i)^SK is the simple Kriging of Y(x_i) from the previously simulated values {Y(x₁),…Y(x_i-1)}, σ_SK(x_i) is the Kriging standard deviation U_i is a Gaussian value independent of U₁,… U_i-1. At each stage, the value at one site is simulated and added to the conditioning data, the set of simulated values {Y(x₁),…Y(x_i-1)} has a multi-Gaussian distribution, with mean 0 and variogram g(h), a simulation of the random function is obtained at the sites x₁,…x_n (Giraldo 2002).

This procedure, although complex, allows evaluating the behavior of spatial autocorrelation and forces the exploration of the NDVI and DEM covariate data. Upon completion of this procedure, the DEM issued the Topographic Wetness Index (TWI), this indicator allows the identification of potential wetland concentration sites or water accumulation zones and the Topographic Position Index (TPI) (Sörensen et al. 2006). The conditional simulation process has already been used by Andrades-Grassi et al. (2021), for the improvement of continuous random variables of low spatial quality.

The Kriging predictions equations can be simply extended to obtain multivariable prediction equations. The general idea is that multiple variables may be cross correlated, meaning that they exhibit not only autocorrelation but that the partial variability of variable A is correlated with variable B, and can therefore be used for its prediction, and vice versa (Bivand et al. 2008). This statistical technique requires the availability of information of secondary variables at all points being estimated, co-Kriging is proposed to use non-exhaustive secondary information and to explicitly account for the spatial cross correlation between the primary and secondary variables (Goovaerts 1997; Li and Heap 2008).

Co-Kriging requires that both target and co-variable be, individually, spatially autocorrelated, and in addition that they be spatially cross-correlated. To perform co-Kriging, one must first model the spatial structure of a co-variable and its covariance with the target variable, this is called a co-regionalization, and it is an extension of the theory of a single regionalized variable (Bivand et al. 2008).

The easiest way to ensure this is to fit a linear model of co-regionalization: all models (direct and cross) have the same shape and range but may have different partial sills and nuggets. The Fit uses an iterative procedure of two steps: first, each variogram model is fitted to a direct or cross variogram; next each of the partial sill coefficient matrices is approached by least squares (Bivand et al. 2008). As with the model obtained with univariate indicator kriging the RMSE and Mean Error were estimated. It should be noted that co-Kriging was evaluated by introducing a single covariate (NDVI, TWI and TPI). Finally, the best co-Kriging prediction was transformed into an indicator function, using the previously described indicator Kriging.