Maize Crop Detection through Geo-Object-Oriented Analysis Using Orbital Multi-Sensors on the Google Earth Engine Platform

Abstract

Mato Grosso state is the biggest maize producer in Brazil, with the predominance of cultivation concentrated in the second harvest. Due to the need to obtain more accurate and efficient data, agricultural intelligence is adapting and embracing new technologies such as the use of satellites for remote sensing and geographic information systems. In this respect, this study aimed to map the second harvest maize cultivation areas at Canarana-MT in the crop year 2019/2020 by using geographic object-based image analysis (GEOBIA) with different spatial, spectral, and temporal resolutions. MSI/Sentinel-2, OLI/Landsat-8, MODIS-Terra and MODIS-Aqua, and PlanetScope imagery were used in this assessment. The maize crops mapping was based on cartographic basis from IBGE (Brazilian Institute of Geography and Statistics) and the Google Earth Engine (GEE), and the following steps of image filtering (gray-level co-occurrence matrix—GLCM), vegetation indices calculation, segmentation by simple non-iterative clustering (SNIC), principal component (PC) analysis, and classification by random forest (RF) algorithm, followed finally by confusion matrix analysis, kappa, overall accuracy (OA), and validation statistics. From these methods, satisfactory results were found; with OA from 86.41% to 88.65% and kappa from 81.26% and 84.61% among the imagery systems considered, the GEOBIA technique combined with the SNIC and GLCM spectral and texture feature discriminations and the RF classifier presented a mapping of the corn crop of the study area that demonstrates an improved and aided the performance of automated multispectral image classification processes.

1. Introduction

Mato Grosso state is the biggest maize producer in Brazil, according to Companhia Nacional de Abastecimento (Conab) in its 12° Brazilian grains crop survey (2021/22), presenting an area of near 6.55 million hectares (ha) and 41.62 million tons (t) of grains [1]. It is also worth noting that despite the large domestic consumption, most of the production is destined for exports, which have broken records in recent years [2]. At the state level, the predominance of cultivation occurs on the second harvest, which corresponds to approximately 99% of the maize crop total area [3]. In view of the large production, the crop area’s estimate plays an important role on national demand supply as well as ensuring that transportation and storage capacity are not compromised [4].

Brazilian official crops estimates are based on subjective surveys by Conab and Instituto Brasileiro Geografia e Estatística (IBGE). Most of the agricultural crops data are based on surveys conducted by technical agents, which rely on cultivated areas, production, and economic data from agricultural producers, agricultural inputs sellers, and other related interviews, which are poorly reliable data for such a survey [5].

Field data surveys occur at low frequency and are being gradually reduced by the scarcity of financial and human resources. Moreover, the large territorial extension of the state of Mato Grosso makes the surveys costly and time-consuming. In this way, the use of remote sensing techniques and geographic information systems (GIS) can be applied to avoid the inconveniences related to the search for agricultural production data [6]. In comparison to other productive sectors, agricultural activities face uncertainties and thus demand frequent and large-scale monitoring [7,8].

Remote sensing and geoprocessing are the most use techniques for land use data generation over time, as they allow the evaluation of changes in the landscape [9]. Remote sensing imagery provides important historical series for crop dynamics identification, and spatial resolution allows precision agriculture-level interventions [10] Geographic object-based image analysis (GEOBIA) has emerged as a powerful methodology for image analysis and classification, proving effective in accurately identifying and classifying land cover types, mapping and monitoring deforestation, monitoring vegetation health and growth, predicting crop yields, and classifying remote sensing imagery [11]. By utilizing spectral, spatial, textural and topological features, GEOBIA enables comprehensive image analysis, providing valuable information about the characteristics of and changes in natural and agricultural landscapes [12].

These methodologies offer strong integration with GIS and use advanced machine learning techniques for image classification. By combining the power of random forest (RF) with the insights derived from GEOBIA, which uses a set of decision trees where each tree is trained on a subset of data and features, providing robust and accurate predictions, researchers and practitioners have achieved significant advances in the analysis of remotely sensed imagery, enabling comprehensive understanding and informed decision making [13,14].

The emergence of several cloud computing platforms, which store images captured by a range of satellite sensors as well as geospatial analysis and geoprocessing tools, has expanded access to free images to supports a wide range of remote sensing research. The Google Earth Engine stands out in this scenario [15]. The platform has imagery data from the Sentinel, Landsat, and Terra/Aqua satellites and provides conditions for the development of geospatial algorithms involving large data sets [16].

In view of the economic importance of maize for Canarana and the Mato Grosso state, information regarding the extension of cultivated areas, which is widely disclosed by the official estimates of area and agricultural production made by IBGE and Conab, supports political and economic planning. However, there are currently differences between these agencies’ publications that were evidenced in the 2019/20 harvest. In that crop year, Conab estimated an area of 5414.4 ha [3], while IBGE estimated 5337.3 thousand ha for second-harvest maize cultivation in Mato Grosso [17].

This area difference is one of the factors that contributes to the difference in maize production of about 930,000 tons in the state, which is a significant amount for agribusiness, and this has raised uncertainty in the sector about the numbers released by official agencies. Therefore, this study presents a geoprocessing- and remote-sensing-based approach to mapping maize crops at Canarana in the state of Mato Grosso.

2. Materials and Methods

The mapping of agricultural land with remote sensing tools endorsed by geoprocessing and based on well-stablished, cost-effective, time-efficient, and accurate machine learning was carried out as determined by the definition of approaches of each of the following steps. Here, the maize crops mapping at Canarana in the state of Mato Grosso considered the steps in the following workflow (Figure 1).

2.1. Study Area

The study area comprises the municipality of Canarana, in northeastern Mato Grosso, over the geographic coordinates 12°36′17″ to 13°47′12″ S and 51°22′32″ to 53°06′12″ W (Figure 2). Canarana has an area of 10,855.181 km² and estimated population of 21.842 inhabitants. The average altitude of the locality is 390 m, and it is characterized by a humid tropical climate (Köppen–Geiger climate classification: Aw), with mean temperature of 25 °C. It exhibits two well-defined seasons, namely dry (from May to September) and rainy (from October to April) [18,19], with an average annual rainfall for the 2019/2020 crop year ranging around 1650.54 mm to 1866.10 mm [20]. The predominant soil class is dystrophic red-yellow latosol [21].

Canarana was involved in the Superintendence for the Development of Amazon (SUDAM), originated in 1972 from a colonizing cooperative company (Cooperativa Colonizadora 31 de Março Ltd.a—COOPERCOL) appliance. The purpose of the project was to attract to the regional rural entrepreneurs (large, medium, and small) as well as the multinational and family producers. Currently, the main economic activities of the municipality are cattle ranching, agriculture (rice, maize, sesame, and soy), and agroindustry. In view of agricultural aptness and phytogeography, the Amazon–Cerrado biomes transition can be seen in the Canarana area, where 34.51% is attributable to the Amazon biome and 65.49% to Cerrado [18,19,22,23].

2.2. Pre-Processing

At first, the IBGE cartographic basis for Canarana was followed, based upon the methodological approach, and then Google Earth Engine (GEE) uploading was carried out. Due to processing limitations in the GEE platform resources, the acquisition, segmentation, and machine learning steps for processing geospatial data were performed in the Google Collaboratory platform (Colab) and can be accessed at https://colab.research.google.com/ (accessed on 12 November 2022). The Google Colab platform is based on the Phyton language in the Jupyter Notebook, allowing free access at runtime to scripts. The PlanetScope NICFI imagery was added to GEE in addition to other multispectral imagery present in Google Earth Engine, namely MODIS Terra/Aqua, OLI/Landsat-8 (Operational Land Imager), and MSI/Sentinel-2 (MultiSpectral Instrument). Multiple imagery-based stages in this methodological approach were utilized for maize crops detection.

Regarding the maize phenological cycle, images from 1 April 2022 to 31 May 2022 were considered since this period represents the stage of greatest vegetative vigor typical of the second-crop maize culture, and these images utilized cloud filtering of up to 35%. In order to obtain a single image to represent the collected period, the median was used.

To reduce the effects of terrain irregularity, before calculating the subsequent spectral indices and image classification, it was necessary to perform topographic correction on GEE [24,25]. This correction was based on a semi-empirical method that takes into account the topography of the area and the solar angle (zenith and azimuth), called the correction method Solar Canopy Sensor+C (SCSc) [25,26,27]. The SCS+C model is based on the canopy, allowing changes in the direction of illumination to be considered during the processing of light correction from inclined to horizontal surfaces.

The digital elevation model (DEM) information was taken from the Shuttle Radar Topography Mission (SRTM); the SRTM V3 (SRTM Plus) product was provided by NASA JPL, with a resolution of 1 arc second (approximately 30 m) [28,29]. The SRTM metadata were used to obtain information such as the satellite angle at solar zenith for topographic correction effectiveness.

Mapping maize crops was dependent on the vegetation indices (VI). Therefore, the NDVI (Equation (1)) was considered since this VI creates an image displaying the green tone (relative biomass), illustrating the contrast of the chlorophyll pigment absorption characteristics in the red band and the high reflectivity of plant materials in the near-infrared band (NIR) [30].

The enhanced vegetation index (EVI) (Equation (2)) was also considered. Despite being similar to Normalized Difference Vegetation Index (NDVI), the EVI is improved in regards to the soil and atmospheric effects for vegetation mapping, and it also does not saturate dense green vegetation areas [30]. Further, the perpendicular vegetation index (PVI) (Equation (3)) allows to nullify the background reflectance of soil a at the crop emergence stage; an important part of this reflectance registered by the sensor refers to exposed soil. Finally, the perpendicular crop enhancement index (PCEI) (Equation (4)) was used to determine the minimum and maximum crop development period [30].

$$NDVI=\frac{{\rho }_{NIR}-{\rho }_{RED}}{{\rho }_{NIR}+{\rho }_{RED}}$$

(1)

$$EVI=g\times \frac{{\rho }_{NIR}-{\rho }_{RED}}{{\rho }_{NIR}+\left(c_1\times {\rho }_{RED}\right)-\left(c_2\times {\rho }_{BLUE}\right)+1}$$

(2)

$$PVI=\frac{{\rho }_{NIR}-(a\times {\rho }_{RED})-b}{\sqrt{1+a^2}}$$

(3)

$$PCEI=g\times \frac{\left(MaxPVI+S\right)-\left(MinPVI+S\right)}{\left(MaxPVI+S\right)+\left(MinPVI+S\right)}$$

(4)

where:

ρNIR—Reflectance in the near-infrared spectral range;

ρRED—Reflectance in the red spectral range;

ρBLUE—Reflectance in the blue spectral range;

g—Gain factor (10²);

c1—Atmospheric effects correction coefficient for red (6.0);

c2—Atmospheric effects correction coefficient for blue (7.5);

a—Soil line slope (1.17);

b—Soil line intercept (3.37);

MaxPVI—Maximum PVI value observed during the period of maximum maize crop development;

MinPVI—Minimum PVI value observed in the pre-planting and/or emergence period;

S—Enhancement coefficient. (the value of 10² is assigned for the enhancement of the amount of energy deposited in the active layer of the cell due to the reduction of the reflection intensity.)

2.3. Segmentation and Classification

The segmentation and classification process were initially based on textural feature extraction via gray level co-occurrence (GLCM). The GLCM statistical approach relies on the texture features from the distribution of observed intensity combinations at specified positions relative to each other in the same image (Equations (5)–(9)) [31,32]. This algorithm requires an eight-bit gray-level image generated by linearly combining near-infrared, red, and green bands, which leads to eighteen different textural indices [25,33]. From this, three texture models were selected: inverse difference moment (IDM), sum of entropy (SENT), and dissimilarity (DISS) [34].

$$mean=\frac{1}{\mathrm{♯}{P}_{\mathrm{O}\mathrm{b}\mathrm{j}}}\sum _{(x,y)\in P_{\mathrm{O}\mathrm{b}\mathrm{j}}}{C}_k(x,y)$$

(5)

$$Standard deviation=\sqrt{\frac{1}{\mathrm{♯}{P}_{\mathrm{O}\mathrm{b}\mathrm{j}}}\sum _{\left(x,y\right)\in P_{\mathrm{O}\mathrm{b}\mathrm{j}}}\ast \left(c_k\left(x,y\right)-\frac{1}{\mathrm{♯}{P}_{\mathrm{O}\mathrm{b}\mathrm{j}}}\sum _{\left(x,y\right)\in P_{\mathrm{O}\mathrm{b}\mathrm{j}}}{c}_k\left(x,y\right)\right)}$$

(6)

$$GLCM homogeneity=\sum _{i,j=0}^{N-1}\frac{P_{i,j}}{1+{(i-j)}^2}$$

(7)

$$GLCM dissimilarity=\sum _{i,j=0}^{N-1}{P}_{ij}|i-j|$$

(8)

$$GLCM entropy=\sum _{i,j=0}^{N-1}{P}_{ij}(-\mathrm{l}\mathrm{n}{P}_{ij})$$

(9)

where:

$P_{\mathrm{O}\mathrm{b}\mathrm{j}}$—{(x, y):(x, y) $\in P_{\mathrm{O}\mathrm{b}\mathrm{j}}$} Set of pixels of an image object;

$\mathrm{♯}{P}_{\mathrm{O}\mathrm{b}\mathrm{j}}$—Total number of pixels contained in the $P_{\mathrm{O}\mathrm{b}\mathrm{j}}$;

c_k(x, y)—Pixel value of the image layer (x, y), where (x, y) are pixel coordinates;

i—Row number of the co-occurrence matrix;

j—Column number of the co-occurrence matrix;

$P_{i,j}$—The normalized value in the cell i, j: $P_{i,j}=(V_{i,j}/{\sum }_{i,j=0}^{N-1}{V}_{i,j})$;

$V_{i,j}$—The value in cell i, j of the co-occurrence matrix;

N—The number of rows or columns of the co-occurrence matrix.

After the procedure, the principal component (PC) analysis was carried out in view of data dimensionality reduction to maximize the amount of original information in the smallest number of principal components. Here, the data rely on 16 vectors among bands, vegetation indices, and texture features [25,35].

This process significantly reduces the computational burden for feature extraction by transforming a set of correlated variables (original bands) into distinct uncorrelated variables (principal components) that contain a maximum of primary information, significantly speeding up the maximum likelihood classification process [36]. Moreover, these calculations have been widely applied in remote sensing to classify land use and land cover and their changes [35,37].

Next, the images went through the segmentation step, qualified by the application of the simple non-iterative clustering (SNIC) algorithm [38]. In this algorithm, similar pixels are grouped into image objects, possessing spectral and textural information that will be employed in the classification step. SNIC is an improved version of the simple linear iterative clustering (SLIC) segmentation algorithm that benefits from a non-iterative procedure and imposes the connectivity rule from the initial stage [34].

At first, SNIC initializes the centroids of pixels in the regular image grid, and the pixel distance in dimensional color space and spatial coordinates are used to determine the dependence of each pixel on the centroid. Finally, the integrated spatial and color distances result in efficient, compact, and nearly uniform polygons, providing the identification of objects (clusters) according to the input parameters and generating a multi-band raster including clusters and additional layers containing average values of the input features [25,33,39,40]. The main parameters of the SNIC algorithm are image, size, compression, connectivity, neighborhood size, and seeds [39,40].

The classification based on the random forest (RF) classifier was applied due to its great accuracy classification compared to other classifiers and also considering the dataset traits and the applied methods [25,33,39,40,41,42] in addition to providing a reduced probability of explanatory variables to the training data and fitting them perfectly, regardless of the large number of decision trees, where each element uses a random subset of the training data and a limited number of randomly selected predictor variables [43].

2.4. Segmentation and Classification

We collected 2200 random field points of maize crops and other different land-use and land-cover types. For this, 74 points of water, 8 of cotton, 10 of agricultural expansion area, 9 of recovering area, 6 of urban area, 239 of Cerrado, 7 of crotalaria, 4 of beans, 347 of forest, 247 of sesame, 84 of millet, 725 of second-crop corn, 159 of pasture, 149 of degraded pasture, 14 of fallow soil, 4 of exposed soil, and 104 of sorghum for the RF training and validation were determined (Figure 3). The Locus Map application was used to collect these points with an average precision of 5 m. For each class of samples, we used 70% of the points in the training phase and 30% in the validation phase; these were selected randomly [25,33,40].

After processing and classifying the images and obtaining the thematic maps of the corn areas, we proceeded to the analysis of the numerical confusion matrix, which is responsible for determining the method’s accuracy by comparing the percentage of objects classified by class with the real class verified in the field, indicating a posteriori the correct evaluation and errors among the strata studied. The confusion matrix provides the classification’s overall accuracy (OA) (Equation (10)); the producer accuracy (PA) (Equation (11)) or omission error that, in turn, indicates the probability that the result classified in the image actually represents that category in reality; and the user accuracy (UA) (Equation (12)) or commission error, which indicates the percentage of correctness of a polygon or a true pixel (reference) to have been correctly classified [25,33,39,40,44,45]. In addition to these analyses, kappa coefficient metrics (Equation (13)) were applied to assess the reliability and accuracy of the classified data [46].

$$\mathrm{O}\mathrm{A} (\mathrm{\%})=\frac{\sum _{i=1}^{n}Pii}{N}\times 100,$$

(10)

$$\mathrm{P}\mathrm{A}\left(\mathrm{\%}\right)=\frac{P_{ii}}{P_{+i}}\times 100$$

(11)

$$\mathrm{U}\mathrm{A}\left(\mathrm{\%}\right)=\frac{P_{ii}}{P_{i+}}\times 100,$$

(12)

$$\mathrm{k}\mathrm{a}\mathrm{p}\mathrm{p}\mathrm{a}\left(\mathrm{\%}\right)=\frac{N\times \sum _{i=1}^{n}Pii-\sum _{i=1}^n(P_{i+}\times P_{+i}) }{N^2-\sum _{i=1}^n(P_{i+}\times P_{+i}) }\times 100$$

(13)

where:

n—Total number of columns in the confusion matrix, i.e., the total number of categories;

P_ii—Number of correct classifications of the top crop type sample in row i and column i of the confusion matrix;

P_i+—Total number of samples of the crop type in row i;

P_+i—Total number of samples of the culture type in column i;

N—Total amount of samples used for verification.

3. Results and Discussion

Through the digital elevation model, two land-use classes were derived with a 12% slope threshold based on the mechanization factor as one of the attributes to be observed in the agricultural aptitude and considering that some machines, mainly harvesters, that are available on the market are adapted for slopes of up to 12% [47,48,49]. The area amount of 10,465.36 km² with up to 12% slope was obtained, corresponding to 96.55% of the total area of the municipality.

For improving supervised classification results, PC analysis was carried out, and as input, the reflectance data of 16 bands based on red (R), green (G), blue (B), and near-infrared (NIR) bands of the four sensors as well as the terrain slope and the vegetation indices were utilized. Here, the first three components (PC1, PC2, and PC3) were defined, which can express the dataset variance (

Table 1. Principal component variation in each sensor.

PC	OLI/Landsat-8	MODIS	PlanetScope	MSI/Sentinel-2
PC01	97.02%	94.52%	98.22%	95.28%
PC02	2.65%	4.31%	1.52%	4.51%
PC03	0.29%	1.08%	0.18%	0.15%
PC04	0.03%	0.04%	0.04%	0.03%
PC05	0.01%	0.03%	0.03%	0.02%
PC06	0.00%	0.02%	0.01%	0.01%
PC07	0.00%	0.00%	0.00%	0.00%
PC08	0.00%	0.00%	0.00%	0.00%
PC09	0.00%	0.00%	0.00%	0.00%
PC10	0.00%	0.00%	0.00%	0.00%
PC11	0.00%	0.00%	0.00%	0.00%
PC12	0.00%	−0.00%	0.00%	0.00%
PC13	0.00%	−0.00%	0.00%	0.00%
PC14	−0.00%	−0.00%	0.00%	−0.00%
PC15	−0.00%	−0.00%	0.00%	−0.00%
PC16	−0.00%	−0.00%	−0.00%	−0.00%

).

Our results corroborate other findings [33,35,50] that obtained more than 90% of the information of the original bands in the first three principal components, which can be expressed in false-color composition as shown in Figure 4.

Hereinafter, the segmentation based on first three principal components was carried out. The combination of texture features with vegetation indices, the original bands, and the PC analysis improved the segmentation results [33].

In the segmentation, one of the factors that may have influenced the results was the decrease of the spatial resolution of the MSI/Sentinel-2 and Planet NICFI sensor images to 30 m, making it possible to evaluate the spectral, textural, contextual characteristics, and hierarchical features of all multispectral images. Also, the limitation of the processing capacity of the GEE platform in the free account was a constraint since the platform is of global use, and therefore free access should be limited [33,51,52,53].

In a previous step of the final classification of the algorithm, the accuracy test with different numbers of decision trees in random forest showed different results for the four sensors used, which were analyzed a range from 10 (ten) to 200 (two hundred) trees as the ideal combination in terms of accuracy and computational costs [33].

In view of highest accuracy of random forest classification, the optimal quantities of decision trees for imagery sensor were 120, 110, 150, and 180 for OLI/Landsat-8, MODIS, Planet NICFI, and MSI/Sentinel-2, respectively (Figure 5). This process concomitantly evaluated the importance of each image band for land-use and land-cover classification.

From the classified data (Figure 6 and Figure 7), the confusion matrices were generated with OA from 86.41% to 88.65% (Figure 8, Figure 9, Figure 10 and Figure 11;

Table 2. Overall accuracy (OA), kappa coefficient, and classified maize crop area (ha) for each sensor image classification.

	Landsat-8	MODIS	Planet	Sentinel-2
Overall accuracy	88.65%	86.83%	86.79%	86.41%
Kappa coefficient	84.61%	82.01%	82.06%	81.26%
Second-harvest maize area (ha)	450,766.60	424,715.59	329,557.85	432,422.91

). The Landsat imagery had the highest accuracy, with OA ranging from 88.65% to 84.61% and PAs between 65.16% and 91.53%. Yet, the lower accuracy relied on the “other land uses” class, and the higher corresponded to the second-harvest maize class, where UA for this class was 91.98%. [53,54,55].

The classification of the second-harvest corn areas with Planet images deviated from the results of the other sensors, with a total area for this class of 330,000 hectares, indicating that despite the higher spatial resolution, the classification algorithm confused the second-harvest corn crop with the other second-harvest crop (Table 2). Considering the analysis of variance and Tukey test’s for means comparison, there was no significant difference at 5% probability among the results of the classifications performed, implying that the responses of the sensors were statistically similar (

Table 3. Analysis of variance (ANOVA) and means multiple comparison by Tukey’s test at 5% probability for the overall accuracy and kappa coefficient results.

Analysis of Variance
Variable	N	Mean	SD	SE	95% Conf	Interval
MODIS	2	84.42	3.42	2.42	53.72	115.12
Landsat-8	2	86.63	2.86	2.02	60.95	112.31
Planet	2	84.43	3.35	2.37	54.37	114.49
Sentinel-2	2	83.84	3.64	2.57	51.18	116.50
Tukey’s test HSD
Group 1	Group2	Meandiff	P-adj	Lower	Upper	Interval
Landsat-8	MODIS	−2.2098	0.9	−15.7536	11.3339	False
Landsat-8	Planet	−2.2037	0.9	−15.7474	11.34	False
Landsat-8	Sentinel-2	−2.7929	0.8225	−16.3367	10.7508	False
MODIS	Planet	0.0061	0.9	−13.5376	13.5498	False
MODIS	Sentinel-2	−0.5831	0.9	−14.1268	12.9606	False
Planet	Sentinel-2	−0.5892	0.9	−14.1329	12.9545	False

).

Despite the lower spatial resolution in relation to the other sensors, MODIS images showed values very close to those of the OLI/Landsat sensor, obtaining OA and kappa coefficients of 86.83% and 82.01% (Table 2), respectively, with the focus class of the study having CA (consumer accuracy) of 90.85% and PA (producer accuracy) of 89.40%.

The independence of spatial resolution with LULC classification accuracy were seen in the comparison between the OA and kappa coefficient results of MODIS and MSI/Sentinel-2 sensors, indicating that high-spatial-resolution data are not always superior to low-resolution data in identifying land uses when objects have varied and complex attributes [25,33,52].

The PlanetScope imagery, despite its high spatial resolution, reached inferior overall results compared to OLI/Landsat-8 [33]. This result is probably related to the lower temporal and spectral resolution of this dataset and to the training points being collected in a concentrated manner and not spread over the entire area of interest [56]. This combination of factors did not produce enough information to adequately differentiate the classes, suggesting that the spatial resolution of 4.77 m did not provide enough textural information for the effective separation of the classes studied. The results of Planet NICFI images in relation to Sentinel-2 were similar to each other in both the OA of 86.79% and 86.41%, and kappa of 82.06% and 81.26%, respectively. For CA, the results were higher than 87% and for PA higher than 89%.

Classification in a highly heterogeneous and fragmented agricultural region is challenging due to the similarities in reflectance among the second-crop crops as well as the different temporal and spatial resolutions of the images used. Given this situation, the combination of several types of images with spatial and temporal resolution as well as the use of time series images can be seen in past research, in which better results were achieved than those obtained by analyses using separate, once-monthly images [55,57,58].

Based on our studies, the result obtained in this research, especially for Landsat 8 and MODIS Terra images, is very consistent with respect to classification accuracy and suggests the use in further studies of the combination of the aforementioned images with PlanetScope NIFCI and Sentinel-2 images as well as other sensors in order to achieve more refined results than the present study.

The results of this study were obtained due to the open-access remote sensing datasets and the computing power of the GEE platform, which has demonstrated considerable versatility and adaptability due to its integrative capabilities and its efficient platform for scripting in JavaScript and Python. This study was based on the use of open-source technologies with a focus on processes, and these are robust and scalable over large spatial extents; thus, good results were obtained.

4. Conclusions

The evaluations performed showed satisfactory results since the kappa coefficient and OA presented values higher than 80%. Therefore, it is concluded that the GEOBIA methodology, which employed the combination SNIC + GLCM with the random forest classifier, was successful. However, it is important to highlight that there is still room for improvement in the segmentation step in order to make the methodology applicable in large agricultural areas. The use of other techniques with more spectral models as well as using other types of sensors in the geo-object-oriented analysis input could improve our findings. For future work, the use of hyperspectral data to improve identification should be considered.

Furthermore, this study highlighted the overall reliability of the GEOBIA methodology, although its complexity results in higher computational demands. This can affect the execution of the GEE code, especially when using high-spatial-resolution data such as that from Sentinel-2 and PlanetScope.