Geostatistical analysis of Arabic coffee yield in two crop seasons

Análise geoestatística da produtividade do cafeeiro arábica em duas safras R E S U M O Para tornar a atividade cafeeira competitiva é utilizada, por alguns produtores, a cafeicultura de precisão. Desta forma é possível a criação de mapas temáticos que direcionam as práticas de manejo para as regiões em que há limitações para o desenvolvimento das plantas. Objetivou-se, neste trabalho, identificar a dependência espacial da produtividade de plantas de uma lavoura de café, em 2012 e 2013. A área experimental está localizada em um Latossolo Vermelho distrófico argiloso, em Três Pontas, Minas Gerais. Foram georreferenciados 100 pontos para a coleta dos dados da produtividade, por meio de derriça manual em pano. Também foi avaliada a diferença da produtividade entre as safras. Os dados foram processados por meio de análise geoestatística. Foi possível identificar e caracterizar a dependência espacial das duas safras, bem como a criação dos mapas de isolinhas. Observaram-se diferenças entre os mapas de 2012 e 2013 devido ao ciclo fenológico bienal do café, o que pode ser confirmado pelo mapa da diferença entre as safras. Recomenda-se um manejo da lavoura que considere a variabilidade espacial da produtividade, para maior retorno econômico.


Introduction
Coffee cultivation is one of the agricultural activities of greatest importance in the Brazilian agribusiness.According to a survey conducted by CONAB (2016), Minas Gerais is expected to produce 28.5 million processed sacks in 2016, leading the national production, due to the increase of area and yield gain resulting from the biennial cycle of the crop, in relation to 2015.Technological innovations in all coffee production stages have been developed and evaluated by researchers and rural producers to make coffee cultivation more competitive in the market, through low production cost, achievement of optimal yields and reduced environmental impacts.For that, precision coffee cultivation stands out, which allows to identify in the plantation the sites that require specific managements.
Despite being available in grain plantations, the yield sensors that allow to quantify the harvest are not common in crops like coffee.Leal (2002) identified the spatial variability of coffee yield using an automatic grain-weighing system and a GPS receiver, in a coffee harvester.The use of yield measuring systems in coffee harvesters also allowed to identify the yield variability, as reported by Sartori et al. (2002).Despite the use of these technologies, coffee yield maps are currently created through georeferenced sampling and manual harvest, allowing a more detailed knowledge on the plantation.Thus, it is possible to act in a localized way, in regions with distinct characteristics of yield.
Since the georeferencing of the plantation is performed considering the entire area as homogeneous, the hypothesis of this study is that there is no structure of spatial dependence of the yield of plants in a coffee plantation.This study aimed to identify and characterize the spatial variability of coffee yield, in 2012 and 2013, and visually compare the thematic maps.

Material and Methods
The experiment was carried out at the Brejão Farm, in a dystrophic Red Latosol -LVd, with clayey texture, located in the municipality of Três Pontas,southern Minas Gerais,Brazil (21º 25' 58" S;45° 24' 51" W;914.7 m).The experimental area is 22 ha cultivated with coffee (Coffea arabica L.) cv.'Topázio' , and the plantation was installed in 2005, at spacing of 3.8 m between rows and 0.8 m between plants.
Except the crop seasons of 2007/2008 and 2008/2009, which received different fertilization based on the spatial variability, as described by Ferraz et al. (2011), the other seasons were managed in the conventional way.In the conventional management, the plantation is fertilized in three periods along the year: February, with 250 kg ha -1 ; October, with 500 kg ha -1 of 20-05-20 and December, with 350 kg ha -1 of 25-00-25.Soil correction is made with the application of 1.0 t ha -1 of dolomitic limestone, with relative neutralizing value (RNV) of 80%, in August.Foliar application of micronutrients is performed in August, October, December and January, at the doses of 0.6, 0.8, 0.8 and 0.6 L ha -1 .Weeds were controlled through the application of 2.0 L ha -1 of herbicide in September, and in December and January mechanical mowings are performed.Diseases are controlled through the application of 0.6 L ha -1 of systemic herbicide in August and January, and 0.8 L ha -1 in October and December.
The experimental area was georeferenced using the geodetic GPS receiver Topcon FC-100 (Topcon Positioning Systems Inc, Livermore, California, USA) aided by a receiver positioned on a fixed base, for the post-processed correction of the data.100 sampling points (coffee plants) were georeferenced.Coffee yield was evaluated in the crop seasons of 2012 and 2013, as well as the difference between seasons.Coffee yield (L plant -1 ) was determined through manual threshing on a canvas and later wagging, which required a container graduated in milliliter (mL).Yield was determined through the mean of observations of four plants formed by the plant relative to the point referenced with the GPS receiver, one plant beside it in the same row, one plant in the row below and another in the row above the referenced plant.The difference between the yield of the seasons was determined by subtracting the data of 2012 from those of 2013.
The results were analyzed through descriptive statistics, Pearson's linear correlation and geostatistics.The descriptive statistics determined the minimum value, maximum value, mean, median, standard deviation, coefficient of variation, coefficient of skewness and coefficient of kurtosis, besides the Kolmogorov-Smirnov normality test.The correlation between the yields of 2012 and 2013 was determined.Geostatistical analysis was performed through the identification of spatial variability, using fits of semivariograms by the classical estimator (Vieira, 2000).
The semivariograms were fitted using the Ordinary Least Squares (OLS) method and the spherical model.The selection of the spherical model is justified because it is the one that best adapted to the fits in studies in the field of geostatistics related to the coffee crop (Alves et al., 2009;Ferraz et al., 2012a, b, c;Silva et al., 2008).From the fit of a mathematical model to the calculated values of γ (h), the following coefficients of the theoretical model for the semivariogram are calculated: nugget effect (C 0 ), sill (C 0 + C 1 ) and range (a).In addition, the spatial dependence degree (SDD) of the variables was analyzed according to the classification proposed by Cambardella et al. (1994).The validation (Isaaks & Srivastava, 1989) was used to determine the mean error, standard deviation of the mean error, reduced mean error and standard deviation of the reduced mean error.The validation was performed to verify if the fits in the semivariograms met the requirements and if they were satisfactory.
After the fit of the semivariograms to identify the spatial variability, the data were interpolated through ordinary kriging, allowing to create thematic maps for each variable.Maps of the standard errors of the prediction were also created.
The geostatistical analysis and maps were made using the free software environment R (R Development Core Team, 2014), through the package geoR (Ribeiro Júnior & Diggle, 2001).The maps were created in Universal Transverse Mercator (UTM) coordinates in the zone 23K.

Results and Discussion
The descriptive statistics and the result of the data correlation are presented in Table 1.Regarding the results of the Kolmogorov-Smirnov test, normality was observed only for the yield in 2013 and for the difference between seasons.However, data normality is not required for the geostatistical analysis (Cressie, 1991).For 2012 and 2013, the values of mean and median are close, indicating symmetric distributions, confirmed by the coefficients of skewness close to zero.In addition, the coefficient of skewness of the variables is between -1 and 1; thus, it is not necessary to transform the data before calculating the semivariograms (Kerry & Oliver, 2007).Furthermore, in 2012 the yield varied from 0.04 to 17.46 L plant -1 , with a coefficient of variation of 71.58% (Table 1).In 2013, the yield varied from 0.07 to 21.57L plant -1 , with a coefficient of variation equal to 51.65%.These results show the amplitude of values that these variables may have.There was a high correlation between the seasons of 2012 and 2013, equal to -0.6875, significant at 0.01 probability level.
To understand how the coffee yield behaves spatially, a geostatistical analysis was applied and its results are presented in Table 2.
It was possible to identify and quantify the spatial dependence of all variables, because the absolute value of the difference between two samples increased with the increment in the distance of separation between them, until reaching a value at which there was no more spatial influence, causing the stabilization of the semivariogram.
The nugget effect is an important parameter of the semivariogram and indicates the unexplained variability or measurement errors, considering the utilized sampling distance (Carrasco, 2010).None of the seasons showed nugget effect equal to zero, a condition only observed for the yield difference between the years.Since it was not possible to quantify the individual contribution of these errors, the nugget effect can be expressed as percentage of the sill, thus facilitating the comparison of the spatial dependence degree (SDD) of the studied variables (Trangmar et al., 1985).All variables showed strong SDD, according to the classification proposed by Cambardella et al. (1994).The values of range relative to the semivariograms have a considerable importance in the determination of the limit of the spatial dependence and can also be an indication of the interval between the mapping units, important to optimize future samplings in precision coffee cultivation.The range (a) varied between 155.4 and 204.5 m, for 2013 and 2012, respectively.For the difference between the seasons, the range was equal to 170.4 m.
Interpolation was performed through kriging and allowed the creation of thematic maps of the spatial distribution of coffee yield (Figures 1A and 1C), for the years 2012 and 2013, and of the difference between the years (Figure 1E), in L plant -1 , as well as their respective maps of standard deviations of the prediction (Figures 1B,1D and 1F).
According to Figures 1A and 1C, there was large variation of yield in both maps and regions with lower yields are represented by the dark red color, while higher yields are in regions with dark green color.Regarding the maps of the standard error of the prediction for 2012 and 2013 (Figures 1B,1D and 1F), lower standard errors are close to the sampling points.Yield maps can be used as first step in the identification of limiting and optimal sites, compared with maps of variables of the soil-plant-atmosphere system.In addition, they serve as a tool for the management of harvest, mechanized or manual, and prediction of yield, as described in Ferraz et al. (2012c), besides allowing to plan the post-harvest steps.
In case the mean value (Table 1) is considered to represent the plantation, for both maps, there are sites with very different values of yield.For 2012, the mean value was 7.29 L plant -1 , so that through the generated thematic map, there are large regions on the left and right side of the area that show yield close to zero.The same occurs for 2013, in which the mean value was 9.93 L plant -1 , and there are large regions in the map with values close to 0 (center) and 20 L plant -1 (left side).The yield maps demonstrate the potential problems that may appear in the management of the plantation, when only the mean values are considered for decision-making.It was observed that, even if the management operations (fertilization, control of diseases, among others) are performed in the entire area, this will not result in uniformity of yield in the plantation.Yield variation is influenced by factors such as: biennuality (Camargo & Camargo, 2001); soil fertility and leaf nutrition (Wadt & Dias, 2012;Scalco et al., 2014); occurrence of pests, diseases and weeds (Fialho et al., 2011;Carvalho et al., 2012;Lopes et al., 2012), among others.
In the coffee crop, a significant factor that interferes with the variation of its production, characteristic of its physiological nature, is the biennial alternation, with high and low yields; thus, the plant needs to vegetate in one year to produce well in the next year (Rena & Maestri, 1985).A visual comparison demonstrates the occurrence of biennuality of yield, because regions that in 2012 had the highest yields showed the lowest   values in 2013.Plants that in 2012 produced a lot (regions with dark red color in Figure 1A) used their reserves for fruiting, negatively influencing the growth of branches and, consequently, reducing the yield in 2013 (regions with dark red color).This can be confirmed by the difference between the yields of 2013 and 2012 (Figure 1E), because this map shows that the areas with greater difference, positive or negative, coincide with the areas with higher or lower yield in 2013, respectively.
According to Matiello et al. (2008), the coffee plant undergoes a stage of low metabolism, after high yields.For regions with dark green color in 2012, nutrients must be replaced to the plants using soil amendments and fertilizers, besides foliar fertilization.Hence, plants will find a favorable environment for recovery, promoting stability in the production.However, as the data evaluated in the present study, variables of soil fertility and leaf nutrition may also exhibit spatial dependence (Silva et al., 2010;Silva & Lima, 2015), thus requiring knowledge on them to perform a localized and efficient management in the coffee plantation.

Figure 1 .
Figure 1.Spatial distribution of coffee yield (L plant -1 ): in 2012 (A), in 2013 (C), difference between years (E), and maps of the standard deviation of the prediction for 2012 (B), 2013 (D) and difference between years (F)

Table 1 .
Descriptive statistics of coffee yield (L plant -1 ), for 2012 and 2013 Corr.2013 -Pearson's correlation with the yield of 2013; SD -Standard deviation; CV -Coefficient of variation; Cs -Coefficient of skewness; Ck -Coefficient of kurtosis; KS -Kolmogorov-Smirnov normality test; ++ Correlation significant at 0.01 probability level; NE -Not evaluated; *Significant; ns -Not significant C 0 -Nugget effect; C 1 -Contribution; C 0 +C 1 -Sill; a -Range (m); SDD -Spatial dependence degree; ME -Mean error; S ME -Standard deviation of the mean error; RME -Reduced mean error; S RME -Standard deviation of the reduced mean error

Table 2 .
Parameters estimated for the semivariograms fitted to the spherical model and using the Ordinary Least Squares method for coffee yield (L plant -1 ), for 2012 and 2013