INTRODUCTION
Phenological information supports crop productivity and crop management (^{Sakamoto et al., 2005}; ^{Couto Junior et al., 2013}). Vegetation indices (VIs) are sensitive to phenological changes and have been used to correlate with agricultural productivity (^{Bolton & Friedl, 2013}; ^{Kogan et al., 2013}; ^{Fu et al., 2014}), to estimate attributes such as Foliar Area Index that can be related to agricultural yield (^{Rembold et al., 2013}; ^{Taugourdeau et al., 2014}; ^{Jiang et al., 2014}; ^{Li et al., 2017}; ^{Liaqat et al., 2017}) or to be incorporated into modeling (^{Padilla et al., 2012}; ^{Meroni et al., 2013}; ^{Kowalik et al., 2014}). The incorporation of remote sensing data to the modeling improves the estimation of the agricultural yield, since a multispectral evaluation of the cultures state within a certain area can be obtained favoring the study of the relations of the plant with the environment (^{Delécolle et al., 1992}; ^{Rudorff & Batista, 1990}).
The spectral behavior of an agricultural crop is related to the way interactions occur between electromagnetic radiation and vegetation characteristics such as vegetation structure, phenological stage, vegetation density, spatial orientation and background soil effects (^{Galvão et al., 2005}; ^{Shimabukuro & Ponzoni, 2012}). In areas of perennial crops such as coffee, factors such as soil, systematic use of agricultural implements, internal and interrow shading, and crop seasonality increase the complexity of the study of coffee spectral characteristics (^{Epiphanio et al., 1994}).
The variation in lighting and target geometries, as well as the nonLambertian characteristic of the surface and topography also influence the spectral response of the target (^{Mattar et al., 2014}; ^{Moreira & Valeriano, 2014}). According to ^{Ediriweera et al. (2013)}, the topography influences the radiance measured by the sensors and reduces the accuracy of the estimates derived from the vegetation cover. According to the same authors, the topographic correction on the reflectance values leads to a better prediction of biophysical parameters on vegetation.
As reinforced by ^{Rezende et al. (2014)} and ^{Taugourdeau et al. (2014)}, because it is a complex analysis target, the research on indirect methods of measuring coffee agronomic parameters is scarce. In this context, the objective of this study was to evaluate the relationship between yield of coffee crops and vegetation indexes with and without topographic correction derived from the OLI /Landsat8 sensor for the 2013/2014 and 2014/2015 crops.
MATERIAL AND METHODS
Study area
This study was carried out in the coffee plantations of Fazenda Conquista from Ipanema Coffees group, located in the municipality of Alfenas in the state of Minas Gerais, between 21° 14'57”and 21° 20' 26” South latitude and 45° 54'00” and 45° 59'25” West latitude. The climate of the region is Cwb type with dry winter and mild summer where the average temperature of the hottest month is below 22° C (^{Alvares et al., 2013}).
In relation to the slope classes, 64.4% of the study area is located in undulating relief (8  20 °), 33.5% in smooth undulating relief (3  8 °), 1.8% in flat ground (0  3 °) and 0.3% in heavily undulating (2045 °). Regarding the slope orientation classes, 41.6% of the study area is located in the West facing slopes; 27.9% to the North; 22.7% to the South; and 7.8% to the East. The radial slope orientation frequency profile and slope frequency histogram obtained from the TOPODATA project (^{Valeriano & Rossetti, 2011}), are shown in Figure 1.
Remote sensing data
Surface reflectance values obtained from OLI / Landsat8 (Operational Land Imager) images from bands 4 (0.64  0.67 μm), 5 (0.85  0.88 μm) and 6 (1.57  1.65 μm) from orbit / point 219/75 for the crop years 2013/2014 and 2014/2015 years of low and high productivity, respectively. The dates of acquisition of the images used are shown in Table 1. The phenological stages described by ^{Camargo & Camargo (2001)} in which the authors divide the coffee cycle of the Arabica species into the Southeast of Brazil into: (I) stage of dormancy of floral buds / beginning of flowering: occurs in the months of July  August  September; (II) flowering stage / beginning of grain formation: occurs between the months of October  November  December; (III) stage of grain formation / beginning of maturation: occurs in the months of JanuaryFebruaryMarch.
Crop  Phase  

I  II  III  
DF BF(1)  F BGF(2)  GF  BM(3)  
13/14  Low productivity  07/31/2013  11/20/2013  02/08/2014 
14/15  High productivity  03/08/2014  07/11/2014  01/10/2015 
^{(1)}Dormancy of floral buds / beginning of flowering;
^{(2)}Flowering / beginning of grain formation;
^{(3)}Grain formation / beginning of maturation.
Vegetation index
For this study, three vegetation indices were used: the Normalized Difference Vegetation Index (NDVI), the SoilAdjusted Vegetation Index (SAVI) and the Normalized Difference Water Index (NDWI) (Table 2). NDVI is based on the contrast between the presence of photosynthetic pigments (red reflectance) and the internal scattering of the radiation in the leaf (near infrared reflectance) (^{Rouse et al., 1973}). ^{Huete (1988)} proposed the SAVI that seeks to minimize the influences of soil brightness, incorporating the L factor to NDVI. According to ^{Huete (1988)}, the value of L varies from 0 to 1 according to the density of the vegetation being necessary analysis for different biomes and agricultural situations. In this study was adopted L = 1. NDWI is sensitive to changes in the net water content of vegetation canopies expressing the relation between the reflectance values in the shortwave infrared spectral regions affected by the water absorption, and near infrared (^{Gao, 1996}).
Vegetation Index  Equation *  Authors 

Normalized Difference Vegetation Index 

Rouse et al. (1973) 
SoilAdjusted Vegetation Index 

Huete (1988) 
Normalized Difference Water Index 

Gao (1996) 
^{*}ρ_{nir} = near infrared reflectance; ρ_{red} = reflectance in red; ρ_{swir} = shortwave infrared reflectance; and constant L = 1.
Correction of surface reflectance according to soil topography
In flat soil the influence of the topography on the reflectance is zero. However, in areas in which topography presents variations in slope and slope orientations, the influence of the topographic effect on the spectral data can be significant (^{Moreira, 2014}). The effect of topography can be minimized by taking into account the angle of solar incidence on the surface of the terrain (cosine factor), determined from the azimuthal and solar zenith angles, the slope (zenith angle from normal to the surface) and the hillslopes orientation (azimuthal angle on the normal surface). The calculation of the cosine factor (cos i), according to ^{Sellers (1965)} can be obtained according to [eq. (1)]:
On that,
θ_{s} = solar zenith angle, [degree];
θ_{t} = slope, [degree];
φ_{s} = solar azimuth angle, [degree]; and
φ_{t} = hillslope orientation, [degree].
The cos i can present values between 0 and 1. Values near 0 correspond to dimly lit areas (θ_{s} ~ 90° and elevation angle ~ 0°) and values close to 1 correspond to the illuminated areas (θ_{s} ~ 0° and elevation angle ~ 90°) (^{Lima & Ribeiro, 2014}). According to the results found by ^{Moreira & Valeriano (2014)} it was tested for this study the SCS + C method of ^{Soenen et al. (2005)} (Equations. 2 to 4).
On that:
And,
On that,
ρcorr_{λ} = corrected band reflectance λ;
ρ_{λ} = original band reflectance λ;
b_{λ} = linear coefficient;
m_{λ} = coefficient of the line, and
c_{λ} = constant determined individually for each studied field.
The nonLambertian correction method presented in [eq. (2)] considers the effects of the irradiance variation in each spectral band through constants obtained by the regression method between the reflectance and the angle of solar incidence on the surface. The SCS + C method is based on the C correction of ^{Teillet et al. (1982)} with addition of the SunCanopySensor method (SCS). By including the slope of the terrestrial surface, the SCS + C method considers the influence of the diffuse irradiance on the canopy spectral response. In the ^{Teillet et al. (1982)} method, the constant C_{λ} is used to model diffuse irradiance and exerts a moderating influence on cosine correction increasing the denominator and reducing overcorrection of dimly lit pixels (^{Soenen et al., 2005}).
We used the solar angles (θs and φs) available in the metadata image and the soil angles (θt and φt) that were determined from the digital elevation model (DEM) from the TOPODATA project. After the determination of the cosine factor, the reflectance was corrected for each pixel and spectral band of the image. Since the coffee areas show variations in planting, spacing and plant age, the determination of the cosine factor and the topographic correction were done individually by the coffee field, so that the crop variations on the reflectance values were taken in consideration.
Determination of average values on vegetation index
From the thematic map with the spatial distribution of the coffee plantations given by the technicians from the Conquista farm, it was generated a mask of the coffee plantations with pixels size of 30 m, compatible with the images from the Landsat satellite. On this mask were selected areas with 100% of coffee occupation to obtain reflectance values, thus excluding edge pixels with presence of other targets.
The reflectance values were used to determine the NDVI, SAVI and NDWI index and their corresponding values corrected by topography (NDVI_{corr}, SAVI_{corr} and NDWI_{corr}). Subsequently, we obtained the average values of this vegetation index per plot.
Relationship between coffee productivity and VIs with and without topographical correction
On the determination of the average values of the vegetation index mentioned in item 2.5 per plot, only pure cloudfree pixels were used, totaling 1,185 pixels. It was determinate the relationship between productivity data (t ha^{1}) and the mean values of VIs with and without topographic correction per plot. The analysis was conducted by the simple regression method obtaining the coefficient of determination (R^{2}) which was submitted to the analysis of significance by means of the pvalue of the Ftest for regression analysis at a significance level of 5%.
The quantitative analysis of the validity of the correction method used in the study was made by determining the standard deviation of the data with and without topographic correction, once the topographic effect was removed, the corrected image variance was lower (^{Hantson & Chuvieco, 2011}), besides the calculation of the coefficient of determination between the reflectance values (with and without topographic correction) and the cosine factor, since this relation decreases after the topographic correction (^{Lima & Ribeiro, 2014}).
RESULTS AND DISCUSSION
Analysis of the relationship between VIs with and without correction of topographic effect
Figure 2 shows the frequency histograms on mean values of cos i and zenith angles for the considered period in this study. The mean cosine factor varied between 0.49 and 0.94 with a minimum of 0.33 and a maximum of 0.97. Low values of cos i indicate that some areas receive poor lighting while areas of cos i close to 1 indicate welllit ones. There were no values of cosine factor lower than or equal to zero, so all analyzed plots received direct solar illumination at the moment of the satellite crossing in accordance with the results by ^{Moreira (2014)}. High mean values of cos i observed in periods of low solar zenith (SZ) (lower than 30 °), indicate welllit areas and little shade as stated by ^{Lima & Ribeiro (2014)}.
From the analysis on the values of determination coefficient (R^{2}) (Figure 3) it is verified that, in general, the NDWI index presented higher values of R^{2} with the cosine factor, except for the dates of July 31, 2013 and November 20, 2013 in which the sensitivity of NDVI and SAVI index in relation to the topography varied between plots. On July 31, 2013 SAVI presented the highest values of R^{2} for plots 1, 4, 5 and 8, while NDVI showed higher correlations with cosine factor for plots 2, 3 and 7. On November 20, 2013 SAVI was more sensitive to topography, that is, higher correlation with cos i, except for plot 5 and 8. However, in the year of high productivity (2014/2015), SAVI presented lower correlation with cosine factor, with NDVI being more sensitive to topographic effects.
In relation to the topographic correction (Figure 3) there was decreasing in R^{2} values among the vegetation index and cosine factor (R^{2} values close and / or equal zero) indicating a decrease in the dependence of the index on topographic effects. Similar results were observed by ^{Moreira et al. (2016)}.
Figure 4 shows the standard deviations of vegetation indices before and after the topographical corrections on spectral bands. In general, there was a reduction in the standard deviation values of the NDWI. It was not observed on corrected image, values of standard deviation higher than the original image. For NDVI and SAVI the values of standard deviation of some corrected data were higher than those obtained from the original images indicating inefficiency of the method by ^{Soenen et al. (2005)} to remove the variance from the original image. Similar results were reported by ^{Hantson & Chuvieco (2011)} and ^{Moreira & Valeriano (2014)} in the study of individual bands.
Analysis of the relationship between VIs and coffee productivity
Table 3 shows the values of determination coefficients between three vegetation indexes with and without topographic correction and coffee productivity.
Crop  Date  Phase(1)  NDVI MEDIUM  NDVI_{MEDIUM  Corrected}  

R^{2}  pValue  R^{2}  pValue  
Low  07/31/13  I  0.62  0.02  0.71  0.01 
11/20/13  II  0.58  0.03  0.51  0.05  
02/08/14  III  0.46^{n.s.}  0.06  0.25^{n.s.}  0.21  
High  08/03/14  I  0.89  0.00  0.57  0.03 
11/07/14  II  0.73  0.01  0.81  0.00  
01/10/15  III  0.22^{n.s.}  0.24  0.28^{n.s.}  0.17  
Crop  Date  Phase  SAVI_{MEDIUM}  SAVI_{MEDIUM  Corrected}  
R^{2}  pValue  R^{2}  pValue  
Low  07/31/13  I  0.60  0.02  0.47^{n.s.}  0.06 
11/20/13  II  0.44^{n.s.}  0.07  0.37^{n.s.}  0.11  
02/08/14  III  0.55  0.03  0.43^{n.s.}  0.07  
High  08/03/14  I  0.85  0.00  0.78  0.00 
07/11/14  II  0.88  0.00  0,89  0.00  
01/10/15  III  0.10^{n.s.}  0.44  0.16^{n.s.}  0.33  
Crop  Date  Phase  NDWI_{MEDIUM}  NDWI_{medium  Corrected}  
R^{2}  pValue  R^{2}  pValue  
Low  07/31/13  I  0.70  0.01  0.65  0.01 
11/20/13  II  0.46^{n.s.}  0.06  0.35^{n.s.}  0.12  
02/08/14  III  0.48^{n.s.}  0.06  0.32^{n.s.}  0.14  
High  08/03/14  I  0.89  0.00  0.78  0.00 
11/07/14  II  0.64  0.02  0.81  0.00  
01/10/15  III  0.16^{n.s.}  0.33  0.23^{n.s.}  0,23 
^{(1)}I: Dormancy of floral buds / beginning of flowering; II: Flowering / beginning of grains formation; III: Grains formation / beginning of maturation.
Pvalue for the F test of the regression analysis for significance of 5%; n.s.(not significant)
The data in Table 3 show that agricultural productivity is correlated with the NDVI values in phases I and II, for both low and high production years. In the low production year (2013/14) the determination coefficient values were 0.62; 058; 0.46 (not significant) for phase I, II and III. In the year of high production (2014/15) the determination coefficient was 0.89; 0.73; 0.22 (not significant), for phases I, II and III, respectively, which explained the coffee productivity in 89%, 73% and 22% respectively. The determination coefficient values decreased from phase I to III. The results corroborate with values found by ^{Bernardes et al. (2012)} who found higher correlation between coffee productivity and NDVI in the period of August / September. This fact may be associated to an increase in leaf area index that cause an increase in internal shading and, consequently, a change in crop reflectance.
The determination coefficient of NDVI values after topographic correction increased from 0.62 to 0.71 in phase I in the year of low productivity and from 0.73 to 0.81 in phase II in the year of high productivity.
In general, the SAVI index presented a lower relation with productivity data when compared to the other index. As discussed for the NDVI, in the year of low production the determination coefficient values were lower compared to the year of high coffee production. The determination coefficient between SAVI and productivity was 0.60 and 0.85 in phase I in the crop years 2013/2014 and 2014/2015 respectively. In phase II, the R^{2} was 0.44 (not significant) in the year of low productivity and of 0.88 in the year of high productivity. In the period of grain formation and maturation (phase III), the SAVI presented determination coefficient of 0.55 in the year of low productivity and 0.10 (not significant) in the year of high productivity.
The NDWI presented higher correlation with the productivity data in the dormancy period (phase I) for the twoyear crop with determination coefficient (R^{2}) of 0.70 and 0.89 for the crop years 2013/2014 and 2014/2015, respectively. In addition to the relationship with lower canopy internal shading (^{Epiphanio et al., 1994}), coffee cultivation is highly dependent on water availability which is the main factor affecting productivity (^{Picini et al., 1999}). Thus, index that use spectral information of the medium infrared (energy absorbed by water) may present better results. After the topographic correction was verified correlation with productivity only for phases II and III in the year of high productivity with determination coefficient increasing from 0.64 to 0.81 in phase II and from 0.16 (not significant) to 0.23 (not significant) in phase III.
CONCLUSIONS
The objective of this study was to evaluate the relationship between yield of coffee crops and vegetation index with and without topographic correction derived from the Landsat8 / OLI satellite for 2013/2014 and 2014/2015 crops, and according to the results of the research it can be conclude that the SAVI index presented a lower relation with the productivity data when compared to the other two index. In addition, the three indexes showed higher values of determination coefficients in the years of high productivity, coinciding with the lower index of leaf area, besides greater contribution on the bottom surface.
The analyses of the correlation of the three vegetation indexes with the productivity for coffee trees conclude that the best correlations among the VIs and productivity were obtained for the stages of dormancy and flowering. Among the analyzed indexes the NDWI presented higher determination coefficient (R^{2}) at the dormancy stage for two analyzed years (R^{2} between 0.69 and 0.89).
In two years of study the NDVI obtained the best correlations with productivity in the dormancy and flowering stages. This index explained the productivity variation between 62% and 89% of the data observed in the dormancy stage and from 58% to 73% in the flowering stage.
The topographic correction method reduced the correlation between vegetation index and cosine factor to zero in most of the analyzed dates. However, inefficiency of the method was observed in decreasing the standard deviation of NDVI and SAVI data.
Topographic correction increased the determination coefficient between productivity and the three indexes in phase II in the year of high productivity and in phase I in the year of low productivity for the NDVI index.
It is emphasized that additional studies are needed to evaluate the efficiency of other methods of topographical correction in coffee areas and its implication on the correlation between vegetation indexes and agricultural productivity.