Analysis of the Biennial Productivity of Arabica Coffee with Google Earth Engine in the Northeast Region of São Paulo, Brazil

Abstract

Numerous challenges are associated with the classification of satellite images of coffee plantations. The spectral similarity with other types of land use, variations in altitude, topography, production system (shaded and sun), and the change in spectral signature throughout the phenological cycle are examples that affect the process. This research investigates the influence of biennial Arabica coffee productivity on the accuracy of Landsat-8 image classification. The Google Earth Engine (GEE) platform and the Random Forest algorithm were used to process the annual and biennial mosaics of the Média Mogiana Region, São Paulo (Brazil), from 2017 to 2023. The parameters evaluated were the general hits of the seven classes of land use and coffee errors of commission and omission. It was found that the seasonality of the plant and its development phases were fundamental in the quality of coffee classification. The use of biennial mosaics, with Landsat-8 images, Brightness, Greenness, Wetness, SRTM data (elevation, aspect, slope), and LST data (Land Surface Temperature) also contributed to improving the process, generating a classification accuracy of 88.8% and reducing coffee omission errors to 22%.

1. Introduction

Coffee is a perennial crop whose characteristics can restrict the accuracy of satellite image classification [1]. Its cultivation occurs between the tropics, a region with a high concentration of clouds in the rainy season and significant topographic variations [2]. The spectral signature of the coffee tree is very similar to other crops and native vegetation [3,4,5]. Coffee production can have areas combined with another type of crop [6,7], the size of planting areas can vary [8], the maturity stages of the coffee plantations change significantly [9], and the lighting condition age can change, as they are grown in sun or shade [10].

The characteristics of the images also have a direct influence on the mapping results. The use of medium resolutions, such as Landsat, is widespread [6,7,11,12], as the availability of images and their revisit frequency make it possible to create a longer time series, which is fundamental in agricultural monitoring.

The use of multi-temporal mosaics is recurrent in agricultural mapping. When integrated with meteorological indicators, texture types, and topographic information, the performance of classifiers tends to increase [1,6,7,8,9,10,11,13]. The temporal mosaics of satellite images are important in evaluating different phases of crop development. Arabica coffee, in addition to the phenological cycle marked by different stages of seasonal development of the plant, is characterized by a biennial production cycle, with alternation between more productive and less productive years.

Biennial cycles are also important in other agricultural crops, particularly perennials, due to cyclical variations in productivity [14]. The analysis of biennial periods can contribute to optimizing harvest management and improving crop yield prediction models [9,11,15]. The classification of biennial mosaics can improve the accuracy of satellite image classification in several agricultural monitoring systems [16].

Investigating the influence of this biennial period showed potential to improve the image classification process and refine coffee production estimates [17]. Other agricultural research using the GEE platform has corroborated the importance of data integration and indicated that plant phenology could contribute to improving the quality of classifications of wheat and corn [18,19], rice [20], and grape [21]. Some pioneering research on coffee GEE classification was carried out by [4,8]. In both, the authors used temporal mosaics, temperature data, topography data, and addressed seasonal influences on plant phenology. Although they opted for different methodologies, images, resolutions, and study areas, the results indicated an improvement in the accuracy of coffee classification.

It is considered, however, that knowledge about the classification of satellite images of biennial coffee production can be expanded. Thus, the research analyzes the effects of bienniality in the supervised classification of Landsat 8 images of Arabica coffee in the Media Mogiana region, São Paulo, Brazil, between 2017 and 2023. The biennial and annual mosaics were evaluated for the accuracy of the classification of land use (water, urban area, vegetation, coffee, agricultural area, pasture and eucalyptus) and the errors of omission and commission of coffee, on the GEE platform with Random Forest. This case study aims to answer the following question: what is the relevance of bienniality in the supervised classification of Arabica coffee areas?

2. Arabica Coffee

According to the International Coffee Organization, world coffee production in 2022 was 171,268,000 60 kg bags. Of this total, 57.5% are Arabica coffees, and 42.5% Robusta coffees [22]. In Brazil, Arabica coffee crops are predominant under the sun, compared to shaded coffee. Arabica coffees are sold at considerably better prices than Robusta types due to the superior quality of the drink [23]. Total Brazilian coffee production in 2020 was 65.5 million bags and, in 2021, 60.4 million bags. This reduction can be considered a consequence of the biennial production of Arabica coffee and climatic events in producing areas of the country [22].

The Arabica coffee tree is a perennial woody shrub with a complex phenological cycle [24]. Coffea arabica has a phenological development divided into two vegetative and four reproductive phases: (1) vegetation and formation of flower buds (September to March); (2) induction and maturation of floral buds (April to August); (3) flowering (September to November); (4) fruit granulation (January to March); (5) maturation (April to June); and (6) rest and senescence of the branches (July to August) [25].

This phenological characteristic of Arabica coffee generates a variation in productivity. The years in which the plant is in the first vegetative phase, in which there is a concentration of leaves, are called “negative biennially”, that is, the least productive phase. The years in which the plant is in the second vegetative phase, in which flowers and fruits are produced, are called “positive biennially”, in which there is more expressive productivity [26].

Table 1. Variation in biennial productivity of Arabica coffee in São Paulo.

Cycles	Year	Biennial Period	Bags Harvested (Millions)	Climate Observations
1st cycle	Annual 2017–2018	Positive	6.2	Good weather conditions
	Annual 2018–2019	Negative	4.3	Standard weather conditions
	Biennial 2017–2019	-	10.5	-
2nd cycle	Annual 2019–2020	Positive	6.2	Good weather conditions
	Annual 2020–2021	Negative	4.0	Bad weather conditions (droughts and frosts)
	Biennial 2019–2021	-	10.2	-
3rd cycle	Annual 2021–2022	Positive	3.9	Bad weather conditions
	Annual 2022–2023	Negative	5.0	Recovery from the previous harvest
	Biennial 2021–2023	-	8.9	-

Source: [26].

shows the variation in the biennial productivity of Arabica coffee in the state of São Paulo, from 2017 to 2023, the production of bags of coffee, and the climatic observations of each cycle studied.

Coffee production also varies in relation to the planted area. During the evaluated period, fluctuations in coffee planting in Brazil can be observed in Figure 1, and in the state of São Paulo, in Figure 2.

Exploring the distinct phenological and biennial phases of the plant can be significant in improving the quality of mapping coffee areas [25]. In high-production phases, some plants tend to lose more leaves than in low-production phases. As a result, the difference in coffee productivity can be related to leaf biomass and monitored through time series of satellite images. The production volume, when associated with climatic conditions, can provide additional information about the situation of the coffee plant in situ [17].

3. Materials and Methods

3.1. Study Area

The case study was carried out in an area formed by 45 municipalities in the northeast of the state of São Paulo (Brazil) and covers approximately 13,500 km² (Figure 3). The region is called Média Mogiana because of the Mogiana railway, built at the end of the 19th century, a period of great expansion of coffee culture in the state.

The region studied has average monthly temperatures that vary from 18 °C to 24 °C (64.4 °F and 75.2 °F) and an average annual precipitation of 1378 mm, distributed as follows: rainy season, from October to March, and dry season, from April to September. Concerning geomorphological aspects, the Atlantic Plateau predominates in the region, with altitudes varying between 800 and 900 m. In the eastern region, there are higher areas (around 1300 m) and, in the Peripheral Depression, to the north, altitudes are more modest and vary between 500 and 700 m. These characteristics are significant, as coffee productivity depends on altitude and climate throughout its development cycle. Water deficit, for example, is very harmful during the vegetative and fruiting periods (October and May), but it should occur moderately during harvest, from June to September [25].

3.2. Methodology

The research methodology was designed based on the integration of the following parameters: phenology, seasonality, bienniality, resolution and time series of images, mosaic composition and processing platform.

The incorporation of the development phases of the coffee plant into the image classification process determined the temporal interval of the research. The annual cycles begin in October, when the vegetation and floral buds are formed, coinciding with the start of the rainy season in the study area. This period extends until September of the following year, when flowering and the dry season end. In the following month, in October, the second annual period begins, characterized by the development of the fruits, which continues until September, when harvesting occurs. At the end of this cycle, the plant enters the resting and senescence phase. Thus, the two annual periods were combined to form a complete biennial cycle. And the period from 2017 to 2023 was defined to encompass the last three complete biennial production cycles of Arabica coffee in the study area.

The selection of images and input data were defined in preliminary tests. The Landsat-8 satellite was chosen due to the temporal availability and consistency of the data provided. A total of 429 scenes with less than 80% cloud cover were used, distributed as follows: 70 scenes for 2017–2018, 70 for 2018–2019, 70 for 2019–2020, 72 for 2020–2021, 76 for 2021–2022 and 71 for 2022–2023.

Although images from the Sentinel-1 and Sentinel-2 satellites are frequently used for agricultural purposes [14,27,28,29], they were not used in this research because we intend to develop a methodology that may produce longer time series of data and Landsat images are available with the same resolution since 1985. The advantages of Sentinel-2 over Landsat-8 include a higher revisit frequency (approximately 5 days) and better spatial resolution (10 m) [27,28]. However, its availability of images was limited in the study area in 2017 and early 2018, which would prevent the evaluation of the first biennial cycle (2017–2018).

The Sentinel-1 SAR sensor is very useful in mapping coffee in tropical zones, as it provides data regardless of atmospheric conditions [1,4,8]. SAR data also allow the analysis of vegetation structure and contributes to differentiating coffee production systems in the shade of surrounding vegetation [28]. The use of SAR texture metrics allowed for a more effective distinction between coffee trees grown in shaded and full sun systems [8]. However, the integration of these images still presents challenges. In the case of Sentinel-2, it is necessary to consider the differences in the width of the bands, which can generate alternations in the solar and viewing angles, causing variations in reflectance [28]. In Sentinel-1, the effects of backscatter signal displacement and relief shadows can limit the effectiveness of data from areas with very irregular topography [1]. The integration of multi-sensor data (Sentinel-2, Sentinel-1 and Landsat-8), on the other hand, demonstrates a promising path for future research, as it can increase the efficiency of agricultural mapping [8,14,28].

Numerous methodological references [1,2,3,4,6,8,10,11] recommend the assessment of coffee bienniality using multitemporal mosaics. After evaluating the images, spectral bands, indices, and additional data, the input data were tested and validated. The adopted mosaics represent two distinct sets of information. The best results were obtained with Mosaic 1, which was selected as the reference for the composition of the final mosaics. Its characteristics are:

• Mosaic 1: bands B2, B3, B4, B5, B6, B7, B8, B9, B10, B11, Brightness, Greenness, Wetness, LST, elevation, aspect, and slope;
• Mosaic 2: Brightness, Greenness, Wetness, LST, elevation, aspect, and slope.

The main references of this research are related to research on supervised image classifications in which: (I) combined the spectral analysis of Landsat images and the topographic model to map coffee [6]; (II) predicted the distribution of hybrid zones of shaded coffee trees with high-resolution images [10]; (III) mapped Coffea arabica with Landsat-8 images and the Random Forest algorithm [3]; (IV) integrated spectral mixture analysis and data mining for coffee classification [11]; (V) reviewed remote sensing methods for mapping coffee production [1]. The references that used cloud processing were: (VI) classified shaded coffee with GEE [4]; (VII) developed a meta-analysis and systematic review of GEE for geo-big data applications [2]; (VIII) integrated optical and Sentinel radar data to map smallholder coffee production [8].

Agricultural crop monitoring is based on the classification of time series of images and the composition of mosaics, whose most frequent metrics are averages or medians [2]. High spectral heterogeneity can become a significant problem in unitary satellite images; therefore, the association or fusion of temporal images into mosaics by calculating the median tends to minimize this difficulty [1].

The selected classifier was Random Forest (RF), which creates classification through statistical decision trees [30]. Each tree is generated with random ranges or indices, and the final class of each pixel is assigned by most of the results. The RF showed good accuracy when differentiating vegetation areas and other cropping systems from coffee planting areas [3,27,31,32,33].

The Kauth–Thomas linear transformation [34], or Tasseled Cap method, was used to integrate spectral data from some bands with the scene characteristics. This method helps to highlight certain aspects of the image, generating three basic components: Brightness, which represents the amount of light reflected by the surface; Greenness, which corresponds to the reflection of vegetated areas at the expense of non-vegetated ones; and Wetness, which generates reflection of the humidity on the surface [1,34]. The accuracy in classifying coffee areas has been improved by merging these components [4,35].

Land Surface Temperature (LST) data, extracted from Landsat-8, were also incorporated into the mosaics. These surface temperature data are often used in precision agriculture studies. LST captures thermal radiation emitted by the Earth’s surface and helps detect subtle differences in leaf and canopy density, morphology, biomass, species composition, and canopy water status in coffee plantations [4,6,7,9]. The GEE methodology and algorithms [15] were used to create the LST data of the study area (bands 4, 5, 7, and 8).

Data from the Shuttle Radar Topography Mission (SRTM) was also used. Topographic information and relief assessment are important in coffee mapping, particularly slope, and elevation, as they influence the accuracy of classified areas [1,8].

The definition of land-use classes considered the similarity in spectral responses between coffee trees, vegetation, eucalyptus, and other agricultural areas, mainly perennial crops [7,36,37]. The annual and biennial mosaics were classified into seven land-use classes: water, urban area, vegetation, coffee, agricultural area, pasture, and eucalyptus.

Classification accuracy demonstrates credibility in the methodology and reliability of the final maps. Factors such as the number of points collected, the differentiation of classes, and the sampling strategy directly influence the results [1]. To validate the classes and evaluate the coffee classification performance, the accuracy method based on [38] was used. This method generates an error matrix comparing the reference samples with the mapped classes. The results provide an overview of the overall mapping accuracy and the specific commission and omission values for each mapped class.

3.3. Procedures and Data

Figure 4 shows the pre-processing and mosaic composition steps.

Landsat-8 images (USGS Landsat 8 Collection 2 Tier 1 TOA Reflectance) from 2017 to 2023 were used, with bands B2 (blue) to B11 (TIRS-2). Calculating the median to create the mosaics reduced the image pixel values. In addition to the reflectance bands, LST data and data derived from SRTM were used.

The mosaics were grouped into three cycles, with biennial and annual periods, according to the phenological phases of plant development (positive and negative), as shown in

Table 2. Cycles, biennials, years, and periods (2017 to 2023).

Cycles	Year	Period	Biennial
1st cycle	Annual 2017–2018	1 October 2017 a 30 September 2018	Positive
	Annual 2018–2019	1 October 2018 a 30 September 2019	Negative
	Biennial 2017–2019	1 October 2017 a 30 September 2019	-
2nd cycle	Annual 2019–2020	1 October 2019 a 30 September 2020	Positive
	Annual 2020–2021	1 October 2020 a 30 September 2021	Negative
	Biennial 2019–2021	1 October 2019 a 30 September 2021	-
3rd cycle	Annual 2021–2022	1 October 2021 a 30 September 2022	Positive
	Annual 2022–2023	1 October 2022 a 30 September 2023	Negative
	Biennial 2021–2023	1 October 2021 a 30 September 2023	-

. The following filters were used: dates (annual and biennial), study area, and 80% cloud mask. The corrections adopted were radiometric and scale corrections and pixel quality control (QA_PIXEL and QA_RADSAT bands).

In agricultural landscapes, where the composition of vegetation cover resembles plantations, the sampling strategy and sample collection determine the quality of image classification [1]. Although the incorporation of spectral information is relevant, this measure alone does not guarantee the effectiveness of the classifier. Well-structured sample collection is essential, since the lack of representative samples for each class can compromise the quality of classification results [27].

The samples were collected during fieldwork, starting in 2021, and aimed to differentiate coffee-growing areas from other agricultural crops, both perennial and seasonal. In relation to past periods, samples were created through the visual interpretation of images and aerial photographs [27]. High-resolution images and historical data from the Mogiana Region from the Google Earth Pro platform were used to complement the research samples.

The division of classes in the land use and land cover map was designed to facilitate the discrimination of coffee. It was necessary to separate eucalyptus from agricultural areas and natural vegetation. Samples from each class were spatially distributed across the study area and were individually interpreted. The validation of the results depends on the heterogeneity of the mapped classes, as well as on the sampling strategies and sample sizes [1]. The number of samples should vary according to the study area [1]; however, it is recommended to stratify samples by class to limit distortions caused by less prevalent land cover classes [38].

The samples collected from the seven classes were divided as follows: 70% of the collected polygons were used to draw random training points and 30% to draw random validation points, thus ensuring independence between both. This pattern of training and validation points was maintained for the RF classifier in all sets of mosaics in all years and used in all classifications.

From the collected samples (658 polygons), training and validation points were generated for each class.

Table 3. Land-use classes.

Class	Definition
Agricultural area	Areas with diverse crops, such as perennials (orange trees and vines) and seasonal crops (corn and sugar cane), as well as areas destined for agriculture, were collected (109 samples).
Eucalyptus	Areas of eucalyptus plantations (separate class for distinction from the coffee class) (70 samples).
Urban area	Urban areas, clusters of rural houses, and other buildings (75 samples).
Vegetation	Areas of preserved and/or native vegetation (121 samples).
Coffee	Coffee plantation areas (108 samples).
Pasture	Areas of pastures without agricultural uses during the analyzed period (105 samples).
Water	Water concentration areas (rivers and lakes) (70 samples).

shows the classes and Figure 5 shows examples of samples collected from each class (approximate visualization scale of 1:50,000).

The main recommendations for assessing accuracy and estimating area [38] were carried out. The accuracy of the process was analyzed in two ways: (I) a percentage of total map hits (all classes), and (II) a percentage of coffee omission and commission errors. In addition to only considering successes, the research evaluated errors to highlight what was mapped beyond and below expectations. The observation and variation of errors, mainly those of omission, were important to define the methodological parameters since the reduction in the percentage of this error indicates a better ability to distinguish between land uses that are confused with coffee, such as eucalyptus, agricultural areas, and vegetation. Figure 6 reveals all stages of classification and evaluation of results.

The GEE’s JavaScript code is available at Github (https://github.com/mariaceciliamanoel/Mapping-Coffee-areas/tree/V1.0.0, accessed on 25 September 2024).

4. Results

The classification results for Arabica coffee are presented in

Table 4. General, annual and biennial hits, per cycle (%).

Cycles	1st Cycle			2nd Cycle			3rd Cycle
Year	Annual 2017–2018	Annual 2018–2019	Biennial 2017–2019	Annual 2019–2020	Annual 2020–2021	Biennial 2019–2021	Annual 2021–2022	Annual 2022–2023	Biennial 2021–2023
Mosaic 1	85.3	85.4	88.8	83.8	81	83.7	80.8	81	81
Mosaic 2	79.6	81.3	85.5	77.3	80.3	78.1	74.8	74.5	74.2

,

Table 5. Coffee commission errors per cycle (%).

Cycles	1st Cycle			2nd Cycle			3rd Cycle
Year	Annual 2017–2018	Annual 2018–2019	Biennial 2017–2019	Annual 2019–2020	Annual 2020–2021	Biennial 2019–2021	Annual 2021–2022	Annual 2022–2023	Biennial 2021–2023
Mosaic 1	35.1	34.7	31.8	33.5	38.6	31.7	34.6	38.8	39.9
Mosaic 2	45.3	42.7	43.2	41.5	51.2	46.2	52.5	50.8	56.5

and

Table 6. Coffee omission errors per cycle (%).

Cycles	1st Cycle			2nd Cycle			3rd Cycle
Year	Annual 2017–2018	Annual 2018–2019	Biennial 2017–2019	Annual 2019–2020	Annual 2020–2021	Biennial 2019–2021	Annual 2021–2022	Annual 2022–2023	Biennial 2021–2023
Mosaic 1	25.2	23.8	22	30.1	32.2	25.8	43.7	32.2	45.3
Mosaic 2	50.9	35.1	36.7	26.4	38.5	25.7	46.4	35.4	52.4

. Table 4 shows the overall accuracy of the two mosaics, according to the annual and biennial periods of each cycle. Table 5 and Table 6 contain the commission and omission errors for coffee during the annual and biennial periods of each cycle.

The results of the supervised classification of Arabica coffee are summarized in

Table 7. Results of biennial and annual mosaics (%).

Cycles	Year	General Hits (All Classes)	Commission Errors (Coffee)	Omission Errors (Coffee)
1st cycle	Annual 2017–2018	85.3	35.1	25.2
	Annual 2018–2019	85.4	34.7	23.8
	Biennial 2017–2019	88.8	31.8	22.0
2nd cycle	Annual 2019–2020	83.8	33.5	30.1
	Annual 2020–2021	81.0	38.6	32.2
	Biennial 2019–2021	83.7	31.7	25.8
3rd cycle	Annual 2021–2022	80.8	34.6	43.7
	Annual 2022–2023	81.0	38.8	32.2
	Biennial 2021–2023	81.0	39.9	45.3

. The data show the general accuracy of the classification of the seven land-use classes and the errors of commission and omission obtained in each cycle in the biennial and annual periods.

Figure 7 shows the general accuracy of the classification of the seven types of land use (%) of the three annual and biennial cycles (2017–2019; 2019–2021; 2021–2023), starting in October and ending in September. Figure 8 shows the errors (%) of omission and commission for coffee from the same period.

Maps of all classes and coffee from the annual and biannual periods are shown in Figure 9 and Figure 10. The results reveal discrepancies in the mappings, especially in recent years. Figure 11 highlights the best results observed in the 1st biennial cycle (2017–2019), both for the classification of the seven land-use classes (left) and for Arabica coffee (right).

In Figure 11, a more relevant concentration of coffee plantations can be seen in the eastern portion of the study area, where the relief is higher and has steeper slopes. To the south of the region, spectral confusions were more evident with vegetation, and, to the north, the errors were associated with other types of agricultural use since these municipalities (São João da Boa Vista, Casa Branca, and Mococa) have a greater diversity of crops (perennial and temporary).

5. Discussion

The results confirm the importance of bienniality in coffee classification. The time interval allowed for the comparison of results in three complete biennial periods and showed coherence with the plant’s development cycle. This period can be extended in future research, with observations of cycles before 2017 and after 2023. In this way, the understanding of bienniality and data validation should become even more important.

The composition of the mosaics and the choice of input data were decisive for the quality of the classification of coffee areas. The integration of additional data [1,3,4,6,8], as well as the adoption of strategies to map the different stages and characteristics of coffee, such as age [3], shaded cultivation [4], different sizes of producing areas and topography can increase the quality of the classification [6,7]. In this context, one of the contributions of this research was to highlight the use of individual spectral bands, different from previous works, which emphasized the RGB (B4, B3 and B2), NIR (B5) and SWIR (B6 and B7) bands [3,6,7].

Although the discussion about the division of mosaics into seasonal periods is a consolidated methodological approach [3,4,5,6,7,8,10,11,27], another contribution of this research was the composition of mosaics every two years, when there is complete development of the plant cycle. The results demonstrate that the use of the biennial cycle, instead of the annual cycle, increased the accuracy of classifying coffee planting areas.

5.1. Data Evaluation

The Landsat-8 bands contributed to the multitemporal composition of the mosaics. The availability of data was very important to analyze the last three coffee production cycles, even if the resolution and periodicity of the images are lower than those of Sentinel. In addition to the individual bands, it was found that the integrated use of KT (Brightness, Green, Humidity) data, still little discussed in coffee mapping [4,15], presented great potential, as they helped to discriminate coffee from different seasonal and perennial crops observed at the study site.

SRTM data, whose use is consolidated and is usually used in land-use mapping [1], contributed to coffee mapping, helping to distinguish coffee trees from higher areas and with more significant slopes. Integrating SRTM data into multi-temporal mosaics reduces omission and commission errors in the coffee class [4].

The LST, integrated with the Landsat-8 bands, KT data, and SRTM, contributed to the classification process. This combination improved the discrimination of coffee areas, as it detected subtle changes in the plant, especially in annual periods, when the phenological variation of coffee trees was more significant.

The Normalized Difference Vegetation Index (NDVI) is frequently used to evaluate the spectral response of various crops [3,6,8,11] and presented better results compared to the EVI, NDWI and SAVI indices [1]. The NDVI provided by the Copernicus Global Land Service was considered a relevant source of information to estimate and predict coffee productivity. However, the phenological variables were not evaluated by its authors [39]. Another study used PlanetScope orbital images and productivity data, which resulted in a direct correlation between NDVI and productivity in three coffee harvests. However, the best NDVI performance occurred in specific periods, after harvest, and during the plants’ dormancy phase (July and August) [27].

When evaluating the research results, it was found that the NDVI was not very useful for accurately classifying coffee in annual and biennial evaluations. Despite improving the overall accuracy of mosaics containing all classes by up to 1%, the NDVI increased errors specifically for the coffee class. Commission and omission errors for coffee increased by an average of 1.0% and 1.5%, respectively. As a result, NDVI was excluded from the composition of the mosaics due to the asymmetry between the general hits and the errors it produced.

The increase in NDVI errors may have been influenced by factors unrelated to vegetation, such as the presence of exposed soil or water [27]. Furthermore, it is believed that the presence of different types of trees, such as eucalyptus, native vegetation, and coffee, which show considerable variation in leaf density, may have reduced the effectiveness of this index [29]. One possible approach to mitigate these issues is to combine it with the EVI (Enhanced Vegetation Index) [8,40] and SAVI (Soil Adjusted Vegetation Index) [8,27], which can improve accuracy in areas with exposed soil and dense vegetation. The use of NDVI can be enhanced by utilizing the Red-Edge band (B11 and B12) of Sentinel-2 instead of the NIR band, generating the NDVIre (Normalized Difference Vegetation Index Red-edge) or NDRE (Normalized Difference Red-Edge) [27]. This index composition improved performance in mapping and distinguishing six different crops at more detailed spatial scales [41].

The use of CHIRPS data can contribute to agricultural mapping, especially due to its spatial scale of 5 km and the long time series (since 1981) [42]. However, CHIRPS data depend on meteorological information and, therefore, on the density and distribution of meteorological stations. Another aspect that can limit the accuracy of CHIRPS data is the morphological variation of the relief (irregular topography). It is also necessary to consider the resolution of the images in the integration process with CHIRPS data, as the number of errors can increase when differentiating between coffee grown in the shade and native vegetation [5]. For agricultural applications, integrating CHIRPS data with other climate sources can increase the quality of classifications [5,42].

The use of RF produced good classification results. However, there is a need to consider the distortions and errors generated by the sensitivity of this classifier to spatial autocorrelation in the training samples. This sensitivity highlights the challenges of classifying coffee areas with high precision [1,3,4,8].

There are inherent limitations to the spatial resolution of Landsat-8 images, but, on the other hand, its temporal resolution can be considered advantageous in evaluations of longer periods. Sentinel-2 images are often incorporated into supervised classifications to minimize this Landsat spatial resolution issue [5,8,43]. This incorporation should be evaluated in the future, as studies consider that Sentinel-2’s Red-Edge bands (B5, B6 and B7) can be used to detect small changes in leaf chlorophyll and improve the classification of heterogeneous landscapes [28]. The SWIR bands, despite having similar lengths on Landsat-8 (B6 and B7) and Sentinel-2 (B11 and B12), can be very useful for detecting leaf water content [28]. SWIR data contributed to distinguishing crops and pastures, improving the classification of agricultural areas [28].

The occurrence of clouds, especially in summer, can generate spectral noise in areas not treated by atmospheric corrections and cloud masks. One of the most adopted alternatives to reduce this problem is the combination of data from SAR sensors, both with Landsat-8 and Sentinel-2 images [44,45,46]. By having longer wavelengths, the Sentinel-1 SAR can penetrate clouds and, when combined with other optical sensors, can fill data gaps arising from cloud cover and maintain temporal resolution. Even so, there are limitations associated with steep topography because of the displacement of the sensor’s backscatter signal [1].

Discussions about the phenology of crop types, such as grains, cereals, legumes, and perennial plantings, are common, as they assist in agricultural monitoring, ecological studies, and the assessment of climate change [47]. In the specific case of coffee, the discussions have addressed the concept of seasonality [4,8,11], that is, adapting mapping to the different growth phases of the plant and associating them with distinct climatic variations (dry, wet, warm periods). However, the inclusion of the concept of bienniality has been less developed. Remote sensing techniques have been used to evaluate the influence of biennial production cycles on coffee productivity [15,25]. Nevertheless, the integration of phenology (plant development), seasonality (initial and final months of the phenological phases), and the concept of bienniality (a production cycle that is completed every two years) had not yet been applied to map Arabica coffee areas.

Considering that bienniality is a recurring phenomenon in different types of plantations, especially perennial ones, it is believed that future research could explore its use in other agricultural crops. There are investigations into the influence of bienniality on the development of perennial crops, such as raspberries [16], and on cocoa mapping, aiming to improve the accuracy of classifications [48,49].

It is possible to reproduce the methodology of this research in other Arabica-coffee-growing areas, particularly in countries such as Colombia, Ethiopia and Costa Rica [1]. In Brazil, the methodology can be used throughout the southeast region, as coffee plantations have similar characteristics (phenology, seasonality and bienniality). To classify Robusta-type coffees, it is suggested to use annual mosaics instead of biennial ones. Conilon coffee, a Robusta variety, which predominates in the state of Espírito Santo, is not characterized by biennial production [24].

It is believed that the trend of future research will be to further expand the integration of input data, aiming to distinguish with greater precision the spectral characteristics of Arabica coffee in relation to other types of vegetation and nearby plantations. To this end, it is considered that the use of Sentinel-2 images and texture data from the SAR sensor, from Sentinel 1, can improve the distinction of coffee trees, with the use of mosaics after 2018, especially in small areas of production, in addition to incorporating climate data and other indices that can increase the accuracy of mapping.

5.2. Map Evaluation

When observing the entire period, it was found that the results in the 1st cycle (2017–2019) were the best, especially in terms of general successes. Regarding coffee errors, the 1st cycle showed a decrease in omission errors (22%). The results of the 3rd cycle (2021–2023), on the other hand, were less efficient in the general classification, especially in the annual period 2021–2022 (80.8%). Coffee classification errors from the 3rd cycle also stood out negatively, in the annual and biennial periods. The smallest errors were 32.2% omission and 34.6% commission. This fact can be related to circumstances exogenous to processing, such as production losses, in the harvest and planted areas, caused by climate variations.

The climatic interferences of the 3rd cycle [24] and observed in the reduction in productivity, reveal how the mapping is susceptible to external variables, mainly rain, and temperature, which alter the seasonal pattern of the plant. In the biennial period of the 1st cycle, in which climatic conditions remained within the expected standard, the overall classification accuracy percentage was the best of the entire evaluated period (88.8%).

Unfavorable weather conditions changed the areas planted with coffee during the research period. The occurrence of frosts and droughts, which caused the reduction of coffee plantations, and the severe pruning carried out by farmers, which aimed to save the plants or improve the productivity of the next harvests, influenced these transformations. As a result, the classification process was affected, and errors of commission and omission increased. This increase may also be related to the vegetative loss of the coffee plant, especially in the biennial assessment from 2021 to 2023, whose omission errors reached 45.3%, the highest recorded in the period.

These findings corroborate the considerations of [4] on the relationship between seasonal cycles and climate issues. For the authors, tree vegetation, such as coffee, presents sensitive changes in spectral responses in dry and humid periods. Therefore, the climatic variations may have contributed significantly to the increase in classification errors.

When comparing the biennial and annual periods with the general classification hits (Figure 9 and Figure 10), it was not possible to observe trends or patterns: in the 1st cycle, the highest hits were concentrated in the biennial period 2017–2019 (88.8%); in the 2nd cycle, both the annual period 2019–2020 and the biennial period 2019–2021 presented very similar results (83.8% and 83.7%); in the 3rd cycle, both the annual period 2021–2022 and the biennial period 2021–2023 had an overall accuracy rate of 81%.

However, when comparing the results of the biennial and annual coffee periods, the biennial percentages were more expressive, with fewer errors of commission and omission. Regarding commission errors, there was a reduction of approximately 3% from the annual period to the biennial period in the 1st cycle. In the 2nd cycle, there was a reduction of almost 7% in commission error from the annual period 2020–2021 to the biennial period 2019–2021. Commission errors also reduced from annual to biannual periods: in the 1st cycle, around 3% and, in the 2nd cycle, around 6%.

In the 3rd cycle, there were more significant errors compared to the other cycles, and they did not exhibit similar characteristics. The lowest commission error was recorded in the 2021–2022 annual period (34.6%), while the lowest omission error was observed in the 2022–2023 annual period (32.2%).

The average number of correct answers for the general classification was 82.9% for annual cycles and 84.5% for biennial cycles. In coffee accuracy errors, the average commission error for annual periods was 35.9%, and, for biennial cycles, it was 34.5%. In the omission errors, the average for annual cycles was 31.2% and 31% for biennial cycles.

Coffee plantations worldwide vary in the surrounding land cover and depending on the type of cultivation itself, whether in the sun or shaded [1]. They also vary depending on the types of management, such as irrigation, the existence of areas with windbreaks, the size of the property, and spacing. As these variables influence the classification [36], it was observed that Arabica coffee crops in full sun were better discriminated against when compared to crops close to bordering areas with another class of land-use cover.

In the periods evaluated, errors were also motivated by the diversity of perennial plants (vegetation, agriculture, and eucalyptus). Plant height and age affect processing accuracy [3,8]. Young and recently planted coffee trees were mistakenly assigned as other land-use classes, such as pasture, agriculture, and urban areas. Evaluating and distinguishing the age of the plant can contribute to reducing classification errors since the spectral signature of the coffee tree changes according to its age, as well as changing according to its phenological phase [1,8,15].

Furthermore, the significance of biennial mosaics in classifying Arabica coffee can be indirectly inferred from the findings of other studies. Previous studies in Central America and Vietnam have shown that models for predicting coffee productivity have been improved by including flowering date [50,51,52], as monitoring coffee phenological development and predicting coffee productivity is a complex activity, given the plant’s biannual cultivation cycle [53].

6. Conclusions

The effects of biennial cycles on the supervised classification of Arabica coffee in the southeast region of the State of São Paulo, Brazil, from 2017 to 2023 were examined. The results indicate that the classification of seasonal multi-temporal mosaics improved the level of accuracy of the process when considering the phases of phenological aspects of plant development (from October to September).

The classification accuracy revealed an increase when complete biennial cycles were processed. The most relevant biennial classifications showed 88.8% accuracy in land-use classes and 78% for coffee, both in the 1st cycle (2017–2019). The mosaics with bands B2 to B11, from Landsat-8, associated with KT (Brightness, Greenness, Wetness) information, SRTM (elevation, aspect, slope), and LST data provided the best arrangement for the general accuracy of the biennial classifications, as well as annual. The classification carried out on a shared access platform showed the great potential for analyzing temporal mosaics of coffee.