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 (20-45 °). 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.

FIGURE 1 Radial profile of slope orientation and slope frequency histogram on the study area. 

Remote sensing data

Surface reflectance values obtained from OLI / Landsat-8 (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 January-February-March.

Vegetation index

For this study, three vegetation indices were used: the Normalized Difference Vegetation Index (NDVI), the Soil-Adjusted 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 short-wave infrared spectral regions affected by the water absorption, and near infrared (^{Gao, 1996}).

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 non-Lambertian 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 Sun-Canopy-Sensor 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 cloud-free 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²) which was submitted to the analysis of significance by means of the p-value of the F-test 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}).

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 well-lit 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 well-lit areas and little shade as stated by ^{Lima & Ribeiro (2014)}.

FIGURE 2 Histogram of frequency values of cos i medium and zenith angle (°) in the study area for the images obtained in the year of low productivity: (a) July 31, 2013, (b) November 20, 2013, (c) February 8, 2014; and in the year of high productivity: (d) August 3, 2014, (e) July 7, 2014 and (f) January 10, 2015. (SZ = solar Zenith). 

From the analysis on the values of determination coefficient (R²) (Figure 3) it is verified that, in general, the NDWI index presented higher values of R² 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² 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² values among the vegetation index and cosine factor (R² 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 3 Coefficient of determination (R²) between the cosine factor and the uncorrected VIs, and Coefficient of determination (R² CORRECTED) between the cosine factor and VIs with correction by the SCS + C method for obtained images in the year of low productivity: (a) July 31, 2013, (b) November 20,2013 and (c) February 8, 2014; and for the images obtained in the high productivity year: (d) August 3, 2014, (e) November 7, 2014 and (f) January 10,2015. 

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.

FIGURE 4 Standard deviation of the VIs without corrections (NDVI, SAVI and NDWI) and VIs corrected (NDVI_corr, SAVI_corr and NDWI_corr) by the SCS + C method for obtained images in the year of low productivity: (a) July 31, 2013, (b) November 20, 2013, (c) February 8, 2014; and for the images obtained in the year of high productivity: (d) August 3, 2014, (e) November 7, 2014 and (f) January 10,2015. 

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.

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² 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 two-year crop with determination coefficient (R²) 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.