Citrus orchard mapping in Juybar, Iran: Analysis of NDVI time series and feature fusion of multi-source satellite imageries

Gray Level Co-occurrence Matrix (GLCM) textural features, and spectral features such as vegetation, built-up, bare-soil indices, and the proposed Vegetation Dynamic Index (VDI). A semi-automatic sample selection paradigm is also developed based on a time-series analysis of 12 monthly Normalized Difference Vegetation Indices (NDVIs), Otsu thresholding, multi-level thresholding (MLT), and using two proposed indices called Ever-greenness Index (EGI) and Water-covered or No-vegetation (WCNV) index, and finally human post-revision. The output of the Support Vector Machine (SVM) classification using 60,000 samples and the post-classification operation showed that the classified map has an average overall accuracy (OA) and an average kappa coeffi- cient (KC) equal to 99.7% and 0.992, respectively. The results show that multispectral bands lonely extracted orchards with high accuracy (OA: 99.55%, KC: 0.986), and the rest of the features only made a slight improvement to that. For the year 2019, an area of about 4351 ha is estimated as citrus orchards, which is in accordance with the agriculture department ’ s reports (~4700 ha). The results indicate that from 2016 to 2019, despite a “ citrus to non-citrus ” land-use conversion of about 754 ha, the citrus orchards area was totally expanded by about 17%. Comparing the results with the Google Earth images indicates the precise extraction of orchards with a 10 m spatial resolution. To use the proposed method for extensive cases or areas with other types of evergreen trees, it is recommended to use high-resolution normalized DSMs (nDSMs) and textural features.


Introduction
Nowadays cropland mapping and monitoring, showing the spatial extent/distribution of crops at local, national, and global levels, is of great importance (Mallick et al., 2021;Oliphant et al., 2019;Qiu et al., 2017;Qiu et al., 2018). Accurate and precise information on crops and their behavior, i.e., growth, changes, etc. helps for better agricultural resources management (Breunig et al., 2020;Li et al., 2019). It is clear that these issues are more essential in countries having extensive fertile croplands and a high amount of agricultural production. Today, a variety of the above-mentioned information can be gained using new sciences and technologies like remote sensing (RS) which helps humans in precision agriculture (Sishodia et al., 2020).
Countries that are located near the sea and have fertilized soil play a significant role in the agriculture of the world. Among the agrarian countries, Iran, due to its geographical situation, is the source of a variety of agricultural productions (FAO, 2014). This country is located in the western part of Asia, the Middle East. Although a huge fraction of its provinces has an arid or semi-arid climate, the situation is quite different for the northern parts of Iran such as Mazandaran, Gilan, and Golestan. Having fertilized soil, good weather, high precipitation, relatively high sunny days, 1 etc. leads to the cultivation of a variety of crops all over Mazandaran, as the first-rank province in Iran's agriculture (Mirkarimi et al., 2013;Saadat et al., 2019). In developing countries like Iran, the extent and distribution of agricultural croplands and quality/quantity estimations of cultivated crops have been traditionally managed for many years, but today, the implementation of up-to-date methods and techniques to monitor farms and natural resources is necessary and inevitable. Moreover, the arid and semi-arid climate of the majority portion of Iran has led policymakers to increasingly think about using Geomatics branches such as RS and geospatial information system (GIS) in the field of limited agricultural resource management in the north of Iran (Ahmed et al., 2021;Hu et al., 2021;Jalili et al., 2008). It is noteworthy that the good quality of these farms led to an intense increase in their prices very recently (IRNA, 2020), so more efficient and modern procedures are needed to be utilized to manage them due to their very high prices. In this case, satellite RS as an alternative to the traditional methods can be so effective to help us conduct a comprehensive study on agronomical purposes (Kumar et al., 2008;Shanmugapriya et al., 2019).
A variety of studies have been conducted in the field of crop mapping or crop type classification using aerial and satellite-based data so far. We categorize the literature based on different keywords such as "time-series analysis", "phenology cycle information", "feature fusion", "change detection", "ecological informatics", and "citrus mapping" as follows: • Methods based on time series analysis: A group of studies was conducted based on time series analysis of satellite imageries. For instance, Chen et al. (2018) mapped some types of croplands in Mato Grosso, Brazil using MODIS time series analysis. Hao et al. (2018) improved a neural network-based algorithm for crop mapping (including orchards) using Sentinel-derived Normalized Difference Vegetation Index (NDVI) time series for Hebei province, China. In three distinct but relatively similar studies, Clauss et al. (2018), Bazzi et al. (2019), and Saadat et al. (2019) provided rice crop maps for their case studies using Sentinel-1time series analysis. In another case, Chen et al. (2019aChen et al. ( , 2019b implemented an automatic orchard mapping method using Landsat satellite images time series. Liu et al. (2020) proposed a method for large-scale crop mapping based on multi-source RS images and NDVI time-series analysis in Google Earth Engine (GEE). They concluded that the classification accuracy of the multi-source feature set is noticeably better than the classification outcome derived from single source features. Ofori-Ampofo et al. (2021) mapped crop type from optical and radar time series using a deep learning-based method. • Methods that benefit both phenology information and time series analysis: Zhong et al. (2011) implemented a phenology-based method based on NDVI profile to map crop types in a case study from California. A phenology-based method was proposed by Liu et al. (2018) to map cropping patterns by time-series analysis of MODIS vegetation products. Rad et al. (2019) developed an automatic phenology-based procedure for rice crop detection using the time-series of the optical Sentinel-2 dataset and tested it on some case studies in Iran. Hu et al. (2019) implemented a method based on phenology-based spectral and temporal feature selection to extract crops from satellite time series. Luciani et al. (2019) proposed an automatic procedure for crop mapping, agricultural monitoring, and yield estimation based on phenology information, and implemented that over a case study in Kenya. Momm et al. (2020) developed a procedure to generate a historical crop type based on the time series analysis of Landsat-5 TM-derived NDVIs and crop phenology curves. Pan et al. (2021) developed a phenology algorithm for mapping winter crops using time-series of Sentinel-2 and Landsat images in GEE. Gao and Zhang (2021) summarized the phenology analysis technique for near real-time crop mapping using satellite images. • Papers with multi-source RS data fused scheme: Riedel et al. (2007) integrated optical and Synthetic Aperture Radar (SAR) data to improve land cover mapping in agricultural areas. Ursani et al. (2011) proposed a method based on the fusion of spectral and textural features for tree extraction and mapping other agricultural covers with very-high-resolution (VHR) satellite imageries. Orynbaikyzy et al. (2020) classified crop type using the fusion of different features extracted from Sentinel-1 and Sentinel-2 images. Akbari et al. (2020) implemented a crop mapping scheme on multi-temporal Sentinel-2 images from Tehran, Iran using a fused feature set, metaheuristic feature selection, and an artificial intelligence-based classification. Pott et al. (2021) benefited from the fused Sentinel-1, Sentinel-2, and Shuttle Radar Topographic Mission (SRTM) Digital Elevation Model (DEM) data for crop type classification and mapping in a case study in Brazil. • Literature focused on croplands change detection: Using two SPOT 5 photos obtained in 2003, Li et al. (2014 developed a Support Vector Machine (SVM)-based classification system for orchard identification and change detection. More recently, Sun et al. (2017) used three Landsat images with a 10-year time gap to track dynamic changes in rubber plantations in Xishuangbanna. • Literature from ecological informatics: In the product mapping literature, some studies have focused exclusively on ecological informatics. In this regard, Albayrak et al. (2021) used Random Forest and Decision Tree algorithms, which are among the ensemble learning techniques, and estimated the behavior of migratory bees. Using the results of the modeling study, practical information has been given to the migratory beekeepers with thematic maps for the whole country.
In a similar study, Ayadi et al. (2022) performed a comprehensive survey of resolution methods for constraint networks using static and dynamic algorithms and showed that ecological informatics is the most appropriate constraint satisfaction problem resolution method for a given context. In another study, by applying remote sensing proxies and the breakpoint approach, Royimani et al. (2022) estimated information on the onset of grass aging to understand the forage preparation period and confirmed the value of this method in estimating autumn grass aging and determining its onset. Citrus trees and shrubs (Rutaceae and Aurantioideae) are flowering trees and shrubs found in tropical, subtropical, and mild temperate climates. These plants are among the most productive fruit crops and are adaptable to a variety of climates and soil types, making them valuable economic assets for poverty eradication and agricultural development in many nations. Citrus plantations have become a worldwide phenomenon as a result of rising market demand, particularly in Asia, where over half of the world's citrus fruits are grown (Xiao et al., 2015). For a better understanding of land use management, defining the spatial extent of citrus orchards and recording their temporal changes are essential. Satellite RS offers distinct advantages in terms of detecting cash crop dynamics in a timely and cost-effective manner, and numerous attempts have been made on various scales (Li et al., 2014;Sun et al., 2017). However, representing complicated land surface movements with a few discrete photos simplifies real-world scenarios dramatically. This paradox is especially relevant for citrus orchard dynamics monitoring because the extent of the orchard may change several times due to various impact variables; thus, a small number of photos cannot accurately capture high-order complexity such as acceleration and other non-linearities in time. Managers and scientists must also monitor the progress of citrus orchards at a higher temporal frequency (e.g., annual basis).
Nowadays a great portion of the agricultural products of Iran (as an agrarian country) is devoted to citrus and its derivatives (FAO, 2021). In this regard, mapping citrus orchard extent and distribution on a regular basis using the RS technique is one of the most successful methods for addressing the protection of socioeconomic and agricultural resources in Iran. Due to the prevalence of cultivation of different types of citruses in countries like Iran, Iraq, Lebanon, India, Pakistan, Turkey, etc., the need for new and effective RS-based techniques and algorithms is felt to map and monitor citrus orchards.
In this regard, this paper attempts to propose a robust and applicable RS-based approach for citrus orchard mapping and to the best of our knowledge, no other study has implemented the proposed strategy, neither in regions inside Iran nor outside that. The objective of this study is to provide thematic classified maps representing the spatial distribution of the citrus orchards all over the case study, calculate their area, and detect their changes within a short time. Methodologically, this is conducted using the time series analysis of a group of VIs derived from multi-temporal Sentinel-2 satellite imageries and intrinsic characteristics of citrus. The characteristic is based on the phenology cycle of citrus, i.e., the fact that citrus trees have a relatively high and not variable NDVI value in a period of time. This approach reveals the phenological behavior of vegetations during a period of time. As a result, due to the distinct phenological properties of citrus trees (having evergreen leaves), time-series analysis is used in this research to extract citrus orchards. In summary, our method benefits from the following titles simultaneously: (1) time series analysis of VIs based on the phonology of crops, (2) multi-modal imageries feature fusion, and (3) image classification based on machine learning (ML) algorithms.
The time-series data of a whole year relates the time (in the monthly interval) to NDVI values in such a way that there is no overlap between the values obtained, i.e., each month in the horizontal axis is related to one and only one NDVI value. This approach reveals the phenological behavior of vegetations during a period of time. As a result, due to the distinct phenological properties of citrus trees (having evergreen leaves), time-series analysis provides good features for automatic training feature collection. Furthermore, using multimodal RS data helps to extract citrus orchards as precisely as possible.
The main core of the proposed method is entirely developed in the GEE environment, exploiting its full potential to test different combinations of variables and obtain data collection to improve each image classification step. Moreover, this approach has made it possible to exploit the strengths of supervised classification and implement a new method for the collection of training samples to facilitate the classification and subsequently achieve the highest levels of accuracy. The novelty of this research relies upon the proposed workflow for analyzing different images and variable composites (i.e., monthly image collections), automatic and semi-automatic sample selection, and proposing new spectral indices to improve citrus orchard mapping classification accuracy. Ultimately, the in-depth investigation of these data and their combination enabled us to optimize training and validation points selection, which is a delicate step in the classification process.
We will contribute to exploring whether or not the proposed scheme is able to extract citrus orchards. To do so, by using the classification confusion matrix, visual assessment, and also comparing the results with the published reports from the agricultural organizations, we will figure out how precisely and correctly it works.

Study area description
Juybar as one of the coastal counties of Mazandaran has all the factors that are needed for the cultivation of different crops. In this region, besides rice cultivation, plenty of croplands are devoted to citrus orchards 2 and the people's livelihood intensely relies on them. The citrus products include different types such as orange, tangerine, lemons, and marmalade orange. It is one of the major producers of citrus in Iran and the produced citrus is exported to other provinces and also to other countries (JKMAZ, 2016). The study area of this paper (Juybar) is located in one of the northernmost parts of Mazandaran province in the UTM zone of 39 N with an area of about 285.5 km 2 (Fig. 1a). Its annual precipitation is about 782 mm , and its annual temperature is about 15 • C (MAZMET, 2020). This plain county study site is approximately at the same level as the high seas and based on Google Earth (GE) images, it is seen that it has no jungle and mountain coverages. The geographic coordinates of its center and also its northernmost, southernmost, easternmost, and westernmost points are reported in Table 1.
Considering GE's images and field survey, it is figured out that the citrus orchards are from different categories in terms of size and area. In other words, there are both small and large citrus orchards (respectively about one-tenth of a ha and tens of ha) that belong to farmers and feudalists, respectively. Considering the GE view and SRTM DEM of the study area ( Fig. 1b and c), and as Juybar is located in the coastal and plain areas, citrus is mostly cultivated in flat rural areas (with no rough topography), out of the cities.

Datasets
In this study, three types of satellite images are used, which provide us with data with different modalities. The dataset includes Sentinel-1 SAR images, Sentinel-2 MS images, and ALOS Digital Surface Model (DSM). The data are described as follows. Furthermore, Table 2 reports the characteristics of each dataset, and Fig. 2 provides a representation of them.
• Sentinel-2 MS images: Out of 16 spectral bands of Sentinel-2, 12 top of atmosphere (TOA) reflectance bands including visible (RGB); red edge-1, 2, 3, 4; near-infrared (NIR); and short-wave infrared (SWIR)-1, 2 bands are utilized in this study. The images were pre-processed in terms of radiometric, atmospheric, and geometric corrections, so they are not needed to be pre-processed by the users anymore. Sentinel-2 provides middle-spatial-resolution (10 m to 60 m) images with rich spectral quality, contributes to ongoing MS observations, and provides services and applications such as land management, agriculture and forestry, disaster control, humanitarian relief operations, risk mapping, and security concerns. The main benefit of Sentinel-2 includes the improvement in land-use and land-cover (LULC) classification accuracy (Phiri et al., 2020). This mapping initiative also highlights several other benefits, such as improving production estimates that imply the need for machinery and human 2 Two common crops are cultivated in Juybar, including rice as an irrigated crop and citrus as an orchard crop. Rice is the first-ranked and citrus is the second-ranked crop among all crop types, which on average, annually cover about 32% and 20% of the county area (with a total area of 258.5 km 2 ), respectively. Moreover, about 35% and 20% of the whole croplands (with a total area of about 27,000 ha) in the study area are devoted to rice and citrus, respectively (JKMAZ, 2019).
resources for harvesting and providing valuable information for the entire supply chain (Saiz-Rubio and Rovira-Más, 2020). • Sentinel-1 SAR images: As supplementary data, Sentinel-1C-band SAR Ground Range Detected (GRD) is calibrated, and the ortho-corrected product is used. It is noteworthy that out of different possible single and dual polarizations, the vertical-transmit/vertical-receive (VV) and vertical-transmit/horizontal-receive (VH) bands are available for the study area, so we use the two mentioned bands. In the implementation phase, the instrument mode is set to "interferometric wide swath (IW)".   At the time of data collection, the most recent open-access highresolution (HR) data were limited to Sentinel-1 and Sentinel-2 (the resolution of Landsat was insufficient for this study). However, land-use changes (the eradication of citrus orchards) were observed in the study area during the second half of the last decade (2011− 2020), and fortunately, the oldest Sentinel data belonged to late June 2015, which covered the aforementioned time period. Because the entire dataset from the year is required for the time series analysis, only the 2016 dataset was usable (2015 was not usable). Thus, the study duration was limited to 2016-2019. The information about the data availability is shown here.

Methodology
One of the most important features of citrus trees is being "evergreen" and in order to know whether vegetation is evergreen or not (i.e., have a periodic cultivation calendar), we should extract the vegetation index of the study area for at least a whole year (12 vegetation indices). The Ever-greenness Index (EGI), Vegetation Dynamic Index (VDI), and some other spectral indices are proposed, which are produced using a series of images and by using Otsu thresholding and logical rules. To show that citrus trees have relatively constant NDVI while other crops like rice, wheat, etc. have periodic or non-periodic NDVI changes, the time series of NDVI is extracted from 2016 to 2019. In addition, sometimes, using the term "time-series" for the case in which a series of input images is processed is common (Farhadi et al., 2022;Saadat et al., 2019). In the current study, to obtain propitious features that could be applicable for citrus orchards extraction, analysis and processing of a series/set of input satellite images is required to be able to remove poor pixels (noisy pixels, cloud pixels, etc.). So, the "time-series analysis" sometimes refers to a series of primary images which provide us with a secondary dataset (by conducting preprocessing). For example, in the present study, a series of multispectral images from different dates are filtered in terms of cloudiness, and then the mean value or max value is calculated for corresponding pixels to a certain location but from different time, and this provides us with a clear pixel with high-quality and no cloud cover.
According to the workflow depicted in Fig. 3, all implementations can be divided into three steps: "data preparation and feature extraction", "classification and post-classification", and "validation", which are separately described as follows.

Data preparation and feature extraction
In the data preparation phase, the initial Sentinel-2 TOA reflectance, Sentinel-1, and Alos DSM images are called in the web-based cloud computing GEE platform, which is a professional portal for agricultural applications, vegetation/crop mapping, and monitoring (Mutanga and Kumar, 2019). As mentioned before, the ALOS DSM in the GEE repository belongs to late 2011, so the newest version of ALOS DSM is obtained from the JAXA website and then uploaded to the GEE platform. The call-backed data are pre-processed, i.e., they are filtered by date and boundary, and furthermore, for the Sentinel-2 case, the images with a less cloud percentage are selected out of a variety of cloudy scenes. No other radiometric, atmospheric, or geometric corrections are needed to be performed on the images. The images with the resolutions of 20 m and 30 m are resampled to 10 m and finally, the images are cropped by the shapefile of Juybar (Table 3).
In the feature extraction phase, three types of features are extracted out of original images. They include spectral, textural, and spatial features. The spectral features are obtained from the Sentinel-2 MS bands, and textural features are calculated using the 10 m-bands of Sentinel-2. By spatial feature, we mean the 30 m-DSM which is obtained from the ALOS satellite and discriminates the objects based on their elevations. Two water products that are obtained from the Sentinel-2 images are Normalized Difference Water Index (NDWI) and Modified NDWI (MNDWI) and are defined as Eqs. (1) and (2), respectively, and are used to better detect numerous water resources existing in Juybar. The other type of spectral products is VIs like NDVI (Rouse et al., 1974), SAVI (Huete, 1988), and EVI (Justice et al., 1998), which are calculated according to Eqs. (3)-(5), respectively. Another kind of spectral indices is monthly NDVI, which is produced for all 12 months of the year. For each individual month, the maximum NDVI value (for each pixel location) is extracted from the Sentinel-2-derived NDVI time series using the Maximum Value Composite (MVC) procedure (Holben, 1986). As the MVC method retains only the highest value for each pixel location among all pixels of a series of multi-temporal geo-referenced optical satellite data, it can drastically minimize the cloud effect and led to obtaining high-quality cloud-free optical images during the determined period of time. Due to the dense settlement in the northern part of Iran, especially in Juybar, some well-known built-up indices are calculated such as Normalized Difference Built-up Index (NDBI), Built-up Index (BUI), and Built-Up Area Extraction Index (BAEI) that help us to detect built-up areas. The mentioned built-up indices are defined as Eqs. (6)-(8), respectively. Using the Sentinel-2 bands, Bare Soil Index (BSI) is calculated based on Eq. (9) to be used both in defining other indices and also in classification.
Some other indices are proposed based on the above-mentioned primary indices. Prior to reaching the proposed indices, we describe the phenological behavior of citrus which helps discriminate citrus from other crops having significant phenological changes. In this regard, a scenario is set up in which the NDVI indices are produced for all 12 months. Every 4 months is collected in one sub-dataset named NDVI Jan-Apr , NDVI May-Aug , and NDVI Sep-Dec , which stand for the NDVI for January to April, May to August, and September to December, respectively. We do this because when we try to represent NDVI Jan-Apr , NDVI May-Aug , and NDVI Sep-Dec using a Red-Green-Blue (RGB) color composite system (the indices are devoted to the monitor's red, green, and blue channels, respectively), the citrus orchards are shown in pure white, with an R, G, and B values equal to 255 out of the 0-255 range. This is because of the fact that citrus trees are evergreen. For a better understanding, we go through Fig. 4a which schematically shows how to separate NDVI indices and how to assign them to RGB channels. The stacked NDVI images are presented in Fig. 4b and the time series of their changes are extracted for the period of 2016-2019 ( Fig. 4c) and nine markers are marked on three different croplands. The red marker is marked on three white color lands, and the blue and green markers are marked on six non-white lands. As is shown in Fig. 4c, which shows the time series of NDVI of these three croplands for the period of 2016-2019, there are relatively high and constant NDVI values for samples 4-6 all over the 4 years indicating citrus orchards. On the other hand, due to the fluctuating time series of samples 1-3 and samples 7-9, it is seen that they have meaningful phenological changes (undulation) over the 4 years and they are only cultivated in the specific seasons. For example, the crops of samples 4-6 are cultivated from June to August (JKMAZ, 2016) and no other crops are raised in the mentioned cropland during the year. Our personal information about the growing rice calendar as well as ground truth information from the red-marked land indicates that it is definitely rice cropland and no other crops are raised in the mentioned cropland during the year except rice (Saadat et al., 2019). In this methodology, all other crops like rice that have periodic cultivation calendars are shown in non-white colors. Similarly, orange croplands are devoted to any type of crop (except citrus) which has a periodic cultivation calendar. Considering the above-mentioned descriptions, we propose a binary index that presents the ever-greenness feature of the citrus trees. The proposed indicator which we call it Ever-greenness Index (EGI) is defined as Eq. (10). To obtain EGI, a threshold is applied to monthly NDVIs. The threshold is determined by analyzing some extracted charts similar to Fig. 4 and multi-level thresholding (MLT). In the MLT scheme, NDVI < 0.1 corresponds to barren areas such as sand, rock, etc., 0.3 < NDVI < 0.5 shows sparse vegetation like grassland and shrub, and NDVI > 0.6 presents the dense vegetation. Citrus trees, due to their intrinsic characteristic (having numerous dense leaves), lie in the latest category. In this regard, for defining EGI, firstly, MLT is applied to each monthly NDVI based on the threshold that separates dense vegetation from sparse vegetation (T up = 0.6). We define T up as the upper threshold. The proposed EGI is subsequently utilized for selecting random sample points for image classification which will be described later. Furthermore, by manipulating different arithmetic combinations of NDVIs, it was seen that the summation of differences of each subsequent NDVIs can be a helping feature for discriminating citrus from other types of croplands and other land covers. We named the developed index Vegetation Dynamic Index (VDI), which is defined as Eq. (11). For calculating VDI, the mean NDVI for the whole year (NDVI m ) and the lower threshold (T down  11), the summation of differences between each two subsequent NDVIs is calculated, and the multiplication by the second bracket masks the non-vegetation area from the produced index. It is noteworthy that NDVI m > T down is a binary index that shows areas with any kind of vegetation by "1"s and the non-vegetation areas with "0"s. VDI is a masked grayscale image in which the area with high vegetation changes is shown by relatively bright color tones while the areas with low or no changes are shown by low VDI values, which provide us with a relatively dark color tone. Since citrus trees' leaves do not change in the leaf-off season, VDI is a good spectral feature for discriminating citrus evergreen trees from other vegetations.
where ⋃ and ∑ represent the AND boolean operator and the summation operation, respectively. The "|…|" operator returns the absolute values of NDVI differences. NDVI m is the mean NDVI for the whole year. Respectively, T down and T up discriminate any type of vegetation and dense vegetation from the bare lands. Furthermore, t is a subscript for each individual month (t : 1 − 12 corresponds to January-December, respectively).
In another part of the feature extraction phase, winter-MNDWI and winter-NDVI are calculated. These two indices are from the winter season because in winter, we have the least crops in the croplands, and furthermore, the trees are in their leaf-off mode and all the rice paddy farms are full of water. The Otsu thresholding method, which is commonly used for segmenting images by exhaustively searching for the threshold that minimizes the intra-class variance (Eq. (12)), is utilized for extracting water and bare lands (lands with no vegetation).
where δ 0 2 (T) and δ 1 2 (T) are the variances of the two classes, and w 0 = ∑ T− 1 i=0 p(i) and w 1 = ∑ L− 1 i=0 p(i) are the probabilities of the two classes separated by a threshold T.
By applying the obtained Otsu threshold to the MNDWI and NDVI, water masks and vegetation masks are obtained. Furthermore, bare soil area is obtained by applying the Otsu thresholding to the BSI. Then a binary mask is made by an arithmetic computation between water mask, vegetation mask, and bare soil mask, which defines areas that contain water-covered or no-vegetated areas (at least one of them). This index which we call water-covered or no-vegetated (WCNV) area is defined as Eq. (13).
where OT 1 , OT 2 , and OT 3 denote the Otsu thresholds for producing water bodies and non-vegetated and bare soil coverages. Furthermore, ⋂ is the OR logical operator. The other category of features that are used in this paper is the textural indices group which are calculated based on the average of three 10 m-RGB bands of the Sentinel-2 dataset. We use eight textural descriptors, including entropy, mean, variance, homogeneity, contrast, dissimilarity, second momentary, and correlation, which are calculated from Gray-Level Co-occurrence Matrices (GLCM) (Haralick et al., 1973). The formulas related to the textural indices are collected in Table 4. These features help us to detect citrus orchards and are very useful for discriminating different types of land coverages from each other.
Since the backscattered radar waves are sensitive to the spatial characteristics of phenomena (e.g., objects' height, shape, etc.), the Sentinel-1 images are useful in separating citrus trees from other lowelevation vegetation. As there are numerous urban and residential areas in the study site, Sentinel-1 SAR data are also decided to be used to extract them very precisely (Semenzato et al., 2020). In addition, it helps us to detect water resources (Bioresita et al., 2018). It is noteworthy that the Sentinel-1 data are all from the winter season when the other trees are in the leaf-off mode and rice paddy cropland are full of water, so the Sentinel-1 Radar backscatter provides us with a good feature that is useful for discriminating citrus orchards from other objects.
By calculating all the mentioned features, three feature spaces are obtained which are then merged in the "feature fusion" step. All the above-mentioned primary imageries and indices (features) are stacked to provide us with a hybrid feature space out of multi-modal data. The mentioned feature space which contains all the initial images and the calculated features is then entered into the subsequent phase.

Classification and post-classification
To perform classification, first, we should prepare some training and validation samples for training the ML-based classifier and then its validation. Here we describe the proposed methodology for automatic and semi-automatic sample selection in this paper. Training samples are pointed on the satellite data based on the concept that is described in Fig. 4. In this color composition visualization, which covers all months of the year, citrus orchards are shown in pure white due to the fact that they are evergreen all over the year. Furthermore, given that the remaining lands and farms contain seasonal crops, they are depicted in a non-white color. This is an important feature for distinguishing citrus from other crop types. Regarding this, we proposed EGI and WCNV that help us to collect samples automatically and semi-automatically. To do this, first, EGI and WCNV are converted to shapefiles and then a variety of random samples are automatically created in the study area with the constraint that the selected points should be led into EGI and WCNV shapefiles. The samples are categorized into two classes named "citrus" and "non-citrus". Note that these two class samples are made using EGI and WNCY, separately. To increase the accuracy of sample selection, after the automatic selection is finished, the human operator checks the produced sample feature collection. Specifically, with the aid of GE's images (which are simultaneous with satellite images), the selected samples are checked rapidly and if a sample is collected inappropriately, it will be deleted from the sample collection. Furthermore, to improve the robustness and confidence of sample selection, a variety of samples are collected manually based on the GEE satellite view. The first methodology is automatic while the last one is semi-automatic (automatic selection + human post-revision). As citrus orchards have a relatively regular geometric shape and spatial patterns, this information can be exerted to the classification by the human supervisor as well. A few points (about 40 samples) are also collected in field surveying which is used for visual assessment of classification performance. Afterward, each category of samples (citrus and non-citrus) is randomly split up into two groups of training and validation in a 60%-40% ratio, respectively. Totally, we have "citrus-training", "citrus-validation", "non-citrustraining", and "non-citrus-validation" collections. Of the two defined classes, the "citrus" class includes oranges and mandarin orchards (the appearance of both is similar). The "non-citrus" class encompasses all BSI Rikimaru et al.
water-covered surfaces (such as lakes, rivers, reservoirs, etc.), the general area of the city (consisting of man-made objects such as buildings, roads, etc.), and coverages like bare land, grassland, crops except for citrus, etc. In Section 4, we discuss the accuracy of the proposed sample section method.
As a part of the classification process, we use the supervised Support Vector Machine (SVM) classifier, which is proved to be one of the most successful algorithms to classify the layer-stacked multi-modal dataset (Cortes and Vapnik, 1995;Faris et al., 2018;Joshi and Miller, 2021). There are several classification algorithms used in machine learning; however, SVM is better than most of them since it produces more accurate results. SVM classifiers are more computationally complex than other classifiers, even when the number of positive and negative samples is not equal. SVM could also normalize the data or project it into the Fig. 4. A technique for discriminating ever-green and seasonal crops: (a) Preparation of seasonal and monthly NDVIs (seasonal NDVIs are assigned to the computer monitor RGB color system to provide a visual presentation, and monthly NDVIs are used for time series extraction); (b) Ever-green citrus trees are shown in pure white color when NDVI Jan− Apr , NDVI May− Aug , and NDVI Sep− Dec are assigned to the red, green, and blue channels, respectively 1 ; (c) The high and constant value of NDVI for a citrus orchard vs. other farms 2 which have periodic cultivation schedule. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 1 For clear presentation, the histogram of the figure is manipulated by a 3σ stretch. Also, other visualization parameters were left as those were set by default of GEE, i. e., Range = [− 1,+1]; γ = 1; opaticy = 1. 2 Ground surveying indicated that rice is cultivated in green croplands. Also, orange crops are related to those crops which are cultivated in the winter season.
space of the decision boundary that separates the two classes (Yang et al., 2018). SVM also has the advantage of being able to work in ndimensional space, which makes it superior to other algorithms. Furthermore, SVM does not influence the findings that were expected before or after making minor changes to the feature extracted data (Almansour et al., 2019;Gola et al., 2019). SVM parameters are set based on the default values of GEE's guideline (Radial Basis Function (RBF) was chosen as the kernel type, and the Gamma (γ) and cost parameters were set to 0.5 and 10, respectively). SVM is trained by training samples, and then the hybrid feature space (including the Sentinel-2 MS bands, Sentinel-2-derived VIs, water indices, bare soil indices, built-up indices, Sentinel-1 VV and VH bands, eight textural features, and the DSM) is fed to the classifier.
By conducting the classification operation, a preliminary classified map (PCM) is obtained. After classification, the post-classification operations are used to make the classified image ready for extracting information. The post-processing is done on two levels including (1) overlaying WCNV to the PCM and then (2) applying morphological operations, i.e., opening and closing. In the first post-classification step, the WCNV binary map is overlayed to the PCM according to Eq. (24). According to this equation, a logical operation is performed between the PCM and WCNY binary maps. Let the pixels in the PCM that are related to citrus and non-citrus classes be presented by 1 and 0, respectively. In addition, the 1 s pixels in the WCNV are related to the pixels that are definitely non-citrus (because the 1 s pixels in WNCV are related to the water-covered or non-vegetated areas). By this assumption, the refined classified map (RCM) is a binary image in which the misclassified pixels (the pixels that are wrongly classified as "citrus" class) are corrected.
where RCM and PCM are the refined and the preliminary classified maps, respectively. Moreover, '|' is conditional equality, i.e., it means that RCM is the same as PCM with the condition that the statement in the braces is exerted to the PCM. In the second post-classification step according to Eq. (25), the RCM is firstly post-processed using morphological opening, i.e., performing an erosion followed by dilation using a structural element. This operation is called "opening". The opening operation removes very small misclassified pixels or speckles. Furthermore, the output of the opened image is entered into the "closing" operation, i.e., a dilation followed by an erosion using a structural element. This closes the gaps and holes in the classified map. Generally, the two mentioned morphological operations increase spatial coherency. It is noteworthy that we perform the mentioned chain operations iteratively (an iteration loop of >1) to reach our desired goal.
where FCM and RCM are the refined and the final classified maps, respectively. The opening, closing, erosion, and dilation are shown by "o, •, ⊖, and ⊕", respectively. Moreover, SE denotes the structural element.

Validation
The map cartography is conducted in ArcGIS 10.5 software where the final citrus orchard map is produced as a simple LULC map, and the area of citrus orchards is calculated precisely. Primarily, to precisely calculate the area, the final raster map is converted to a vector file and the numerous dispersed features of citrus are aggregated using the "dissolve" operation.
In the validation step, the accuracy of the produced maps in the previous step is evaluated both numerically and visually in threefold levels including (1) using the confusion (error) matrix, (2) comparing the amount of the calculated area with the agriculture department reports, and (3) visual comparison with the GE images. The threefold levels are described as follows:  Table 5 be a sample for the confusion matrix, then the OA and KC are defined as Eqs. (26) and (27), respectively.
• Visual assessment: A variety of image subsets are obtained from corresponding GE images and the final citrus map to have a visual evaluation of the classification process. The coordinates of the cropped subsets are also surveyed in the field operation. • Validation of the computed citrus area: The obtained findings are compared to some statistics provided by the Ministry of Agriculture, which are collected traditionally using ground-based manners.

Results
After preparing the initial data and filtering them by boundary, date, cloud coverage, resampling, and clipping, we apply an Otsu thresholding in MATLAB to segment the water bodies, non-vegetated areas, and bare lands. The Otsu thresholds for winter-MNDWI, winter-NDVI, and winter-BSI are calculated as about 0.3294, 0.2392, and 0.0275, respectively. Fig. 5 shows the output of Otsu thresholding and the obtained water mask, vegetation mask, and bare soil mask. It is seen that the desired water, non-vegetation, and bare soil coverages are extracted precisely. On the other side, monthly NDVIs are aggregated to produce EGI, which all of them are presented in Fig. 6. As citrus trees are sort of dense vegetation, we exerted a threshold to all monthly NDVIs by an upper threshold, i.e., T up = 0.6.
Eight textural descriptors including entropy, mean, variance, homogeneity, contrast, dissimilarity, second momentary, and correlation are obtained out of the mean band of the three Sentinel-2 visible channels. They are calculated by setting the parameters as follows: the size of neighborhood is set to 1, and a 3 × 3 square kernel is used, resulting in four GLCMs with the offsets (− 1,− 1), (0,− 1), (1,− 1), and Haralick et al.

(1973) Mean
a P(i, j): the relative frequency with which two pixels occur within a given neighborhood; i and j: pixels intensity values; G: number of gray levels.
A. Toosi et al.   (− 1,0). The parameters are set according to the GEE's default. Fig. 7 shows the obtained spectral and textural features. According to this figure, it is seen that different land coverages have different texture patterns and this is useful in classification. The results are then uploaded to GEE. As the first step of the classification process, a variety of samples were collected based on EGI and post-supervision on the GEE's background satellite view (GE images). As mentioned before, by postsupervision, we mean that the human expert checks all the samples rapidly and if a sample is mis-collected, it will be removed from the feature collection. Furthermore, it is noteworthy that the EGI and WCNV are first converted to a shapefile and then they are defined as two regions, and the samples should be randomly and automatically collected in these regions of interest (ROI). Furthermore, we collect a variety of samples manually to be sure our feature collection is robust enough. By doing this, we have two main groups of samples named "citrus" and "non-citrus". Finally, each of the two category sample collections is randomly and automatically divided into training and validation groups according to a 60%-40% ratio. Table 6 presents the details about the total number of samples, the training, and the validation samples. The distribution of samples is represented in Fig. 8. It is worthy to note that after selecting the samples, we checked all the samples carefully and as we see in Fig. 8 (as an example); surprisingly, no mis-collected samples were detected in the automatically collected samples. GE images are used as a reference for the collected sample points and their georeferencing accuracy can vary from 1 m to tens of meters in different locations around the world (Mohammed et al., 2013;Potere, 2008). In this study, its horizontal accuracy was evaluated using ten ground control points (GCPs) with a distribution that covers all over Juybar. The coordinates of GCPs are surveyed in fieldwork using a GNSS receiver in online real-time kinematic (ORTK) mode and with an average horizontal root mean square error (HRMSE) of about 0.029 m. The corresponding of each GCP is also marked on GE Pro software and its coordinate is extracted. Fig. 9b reports the details of coordinates and the accuracies. The HRMSE is then calculated according to Eq. (28). In this study, the total HRMSE is estimated as 2.705 m and this indicates that the GE images have sufficient accuracy to be used in our work, i.e., after all the points are revised by humans on GE images, the sample may have a maximum of about 2.7 m deviations (in the worst situation) from the real world.
where GCP and GE represent the points collected in the field and on the GE images, respectively. Moreover, n is the number of GCPs (n = 10). SVM classifier is trained using training samples and then the hybrid dataset is fed into it. Here, first, we want to analyze the impact of different features in the classification process, so we categorize the dataset into different batches and in each step, only a part of the batches is fed into the classifier. The batches are defined as follows: • batch 1: visible (red, green, and blue) bands of Sentinel-2 data • batch 2: twelve MS bands of Sentinel-2 data • batch 3: MS + spectral indices • batch 4: MS + spectral features + Sentinel-1 bands • batch 5: MS + spectral features + textural features + Sentinel-1 bands • batch 6: MS + spectral features + textural features + Sentinel-1 bands + DSM Six classified maps are produced out of the six mentioned batches. Due to the limited online computation capability and as we want to compare all the six classified maps simultaneously with the same sample points, the classification is conducted using fewer sample points (25,000 samples including 5000 citrus samples and 20,000 non-citrus samples). The training and validation points are produced based on a 60%-40% division rule. The results of the PCM are shown in Fig. 10. Furthermore, the confusion matrix related to the classified maps is presented in Fig. 11 (sub-figures a to f). Moreover, OA and KC of different batches are compared with each other in Fig. 12.
The results of the classification and the confusion matrices show that using only visible bands, we have relatively the poorest classification, but adding invisible bands of Sentinel-2 increases the OA and KC drastically. It is also seen that in this case study, adding other features like spectral indices, textural features, Radar bands, and DSM just made a slight improvement to the accuracy as the MS bands extracted the citrus orchards very accurately. In the poorest situation (batch 1), OA is 92%, which is more than the USGS standard limit (an OA equal to 85%), so all the maps are acceptable . Considering the UA of the citrus class, it is deduced that the majority of the pixels (with an average probability of 94.3%) that were classified as the citrus class actually represent the citrus on the ground; thus, the results are reliable. Furthermore, the citrus class average PA shows that the places on the ground that are citrus orchards were correctly classified on the map with an average probability of 97.6%.
In the second part, we classified batch 6 with the highest possible sample points that the computation capability let us (due to the Internet bandwidth). We do this using 60,000 samples including 15,000 "citrus" samples and 45,000 "non-citrus" samples. The results of the final citrus map are shown in Fig. 13. This map is the FCM which is refined by morphological operators and also overlayed by the WCNV. The confusion matrix related to the mentioned classified map is reported in Fig. 11g, and its OA and KC are estimated as 99.7% and 0.992, respectively. Considering the figure, it can be seen that if we go to the southernmost part of Juybar, the density of citrus orchards will increase. Moreover, planting citrus orchards is more prevalent in the western part compared to the eastern part of the county.
To visually evaluate the quality of the proposed method, a variety of paired image subsets from both GE's HR images and the final orchard maps were provided in Fig. 14. It should be noted that the locations of ground truth data are also checked and surveyed in the field operation. According to Fig. 14, it can be seen that citrus orchards have been precisely extracted considering the spatial resolution of Sentinel data (10 m).
According to the results, the area of citrus orchards for the end of the year 2019 is calculated as 4351.2 ha. It is worth pointing out that the latest statistical information from the Agriculture Department of Juybar, which is related to the end of the year 2020, confirms the findings of this paper. Using the ground-based information collection strategy, they estimated that by the end of the year 2020, about 4700 ha of Juybar's croplands are cultivated by citrus (both small fruitless seedlings and sparse citrus trees) (JKMAZ, 2021), and this is in good accordance with the area calculated in this study using the proposed RS-based technique.
The citrus map is also produced for the year 2016 (using the  methodology that was explained for the year 2019), and the classification is performed by an OA and KC equal to 98.56% and 0.967, respectively (the confusion matrix is provided in Fig. 11h). In this case, GE historical images were downloaded from Google Earth Pro software and were uploaded to GEE as the reference for classification validation.
In the year 2016, about 3716.8 ha of croplands were devoted to citrus orchards. Based on Fig. 15, which shows the change map of citrus orchards from 2016 to 2019, about 1386.3 ha of non-citrus land covers were changed into citrus orchards, and simultaneously, 754.2 ha of citrus orchards were converted into other LULC types. Totally, an increment of about 682.1 ha occurred in the area of citrus orchards.

Comparison of the manual and automatic citrus mapping method
Although the manual mapping of citrus orchards may give better results, the main aim of this research was to detect and extract orchards automatically using the spectral, phenological, textural, and spatial features of citrus. It is considered that citrus is an evergreen plant and has no considerable visible phenological changes within the year, unlike other crops such as rice, wheat, etc., which exist on farms in certain seasons of the year, and in some seasons, the fields are empty of crops, so their VI considerably changes within a year. It is also noteworthy that the manual extraction scheme for large areas such as counties, provinces, and even countries is very time-consuming, tedious, costly, and in some cases, impossible. Furthermore, in some cases, human-based errors may emerge. Thus, we believe that the automatic crop mapping seems to be more reasonable and it is already practiced in several similar studies (e.g., Akbari et al., 2020;Chen et al., 2018;Chen et al., 2019aChen et al., , 2019bLiu et al., 2020). Fig. 16 shows a variety of citrus orchards in GE images and on the other hand, some other croplands are also presented. The interesting thing is that the citrus orchards, due to their special texture and pattern, are quite distinguishable from other types of croplands. Furthermore, as they are evergreen, they can be detected all over the year in the GE images. This provides us with a good chance to select samples via GE images. After selecting samples, the samples were evaluated in the postsupervision phase in which all the points were collected correctly, and surprisingly, no mis-collected samples were detected. This semiautomatic manner, i.e., automatic sample selection and then postrefinement by the human expert, is not that time-consuming. In this study, <45 min were spent collecting data from about 60,000 samples, which is far superior to fieldwork in terms of finance and time savings. The proposed method of this paper can also be used in similar studies which needed to collect plenty of samples as fast as possible. This method has its own pros and cons. The ability to collect large numbers of samples in a very short time is one of the major advantages of this method. On the other hand, as its cons, we can mention that a few points may be collected incorrectly, and the human expert corrects them in the post-supervision step.

Point-wise vs. polygon-based sample selection strategies
It should be noted that due to the speed of the internet and the computational load issue, as well as the more accurate selection of training and validation samples and avoiding the selection of mixed pixels instead of pure pixels, a point-by-point sample selection strategy was used instead of the polygon-based manner, hence leading to a relatively faster classification process and more reliable results. The classified map obtained using the point-wise training and validation samples selection is rationally expected to have relatively more accuracy in comparison with that obtained by a polygon-based sample selection strategy. This is probably because the objects in high-and midresolution images are heterogeneous, so the polygons mostly contain a number of non-intended (mixed) pixels of the objective class. Choosing a random sample of training data from across a landscape ensures that the sample's class proportions are typical of the landscape's real class proportions. The analyst must make methodological decisions that present tradeoffs in data quality (i.e., class representativeness) and quantity while creating training data for RS image classification. The ideal situation is for the samples to be surveyed in the field. However, because these data are sometimes difficult to gather, researchers may be tempted to employ fewer sample points. Time and access constraints in croplands and other remote environments may force researchers to obtain training and validation sample points from imagery rather than field validation, i.e., through image interpretation, as has been done here, which benefits sample selection from both fields survey and GE HR images. Although the certainty of a training sample's class may be poor in some situations, this method allows for the quick and efficient gathering of a higher number of training samples, allowing for a bigger sample collection size.

Analysis of the impact of different feature batches using the McNemars' test
McNemar's non-parametric test (McNemar, 1947) is a popular test that is used to check the superiority of one classifier compared to another (Manandhar et al., 2009;Samadzadegan et al., 2017). In this paper, besides the OA, КC, etc. the McNemar's test has been implemented on the two confusion matrices of batch-1 and batch-6 (subfigures a and f of Fig. 11) based on a binary distinction between correct and incorrect class allocations collected in the contingency table (Table 7) and the standardized normal test statistic defined as Eq. (29). The null hypothesis is that the probabilities p(b) and p(c) are the same; in other words, none of the two models (in our case classified maps) is better than the other. The alternative hypothesis is that the Fig. 10. The PCMs produced out of six batches using the same training sample and the same classifier ((a) to (f) correspond to batches 1 to 6, respectively).
performances of the desired classifier for two classified maps are not equal. Thus, the null hypothesis and the alternative hypothesis are considered as "H 0 : p b = p c " and "H A : p b ∕ = p c ", respectively. The square of Z follows a Chi-squared statistical distribution with one degree of freedom (DoF = 1), so the test statistic is defined as Eq. (30) (Rozenstein and Karnieli, 2011). By setting a 95% confidence interval (α = 0.05 is the significance level), χ 2 is obtained to be 0.99, which is not greater than the critical value (χ 2 (α/2,DoF) = 3.84), i.e., χ 2 < χ 2 (0.025,1) , and also the probability value (p-value) is calculated to be <0.050. This accepts H 0 , so from a statistical point of view, McNemer's test results represent no sign of superiority of batch-6-derived classified map over the batch-1derived map (χ 2 =0.99, p-value < 0.050). It is noteworthy that although the Kappa Z-test is commonly used to find the superiority of one classifier or LULC over another, it is not appropriate if the same case studies are used in the comparison (Manandhar et al., 2009). . 11. The confusion matrices for PCMs produced out of six batches using the same training sample and the same classifier ((a) to (f) correspond to batches 1 to 6, respectively that are produced using 20,000 samples. Also, (g) corresponds to batch 6 which is produced by 60,000 samples); (h) confusion matrix for the year 2016.

Discussion on the computation speed of the proposed method and the role of GEE
The implementation was conducted on a computer with a 2.50 GHz Pentium Intel Core i7 CPU, NVIDIA GeForce MX 130 GPU, and 16 GB of RAM. However, the internet bandwidth is more important than the hardware. With a 16 Mb/s internet speed, GEE scripts consume about 20 min to be run for the case study of the current paper, so the performance on national and continental scales would be time-consuming having about 60,000 sample points and a hybrid feature space containing >30 bands. In the case of web-based client-server platforms like GEE where all computational affairs are conducted on the server-side and only the results are sent to the clients or users, the overall process gets accelerated considerably compared to the desktop-based geospatial data processing platforms.
GEE as a cloud-based platform makes it easy to access highperformance computing resources for processing big geospatial data, without having to suffer the IT pains currently surrounding them (Liang et al., 2020). Once an algorithm has been developed on GEE, users can produce systematic data products or deploy interactive applications backed by GEE's resources (Aryal, 2020;Wu et al., 2020). GEE consists of a multi-petabyte analysis-ready data catalog co-located with a highperformance, intrinsically parallel computation service. It is accessed and controlled through an internet-accessible application programming interface (API) and an associated web-based interactive development environment (IDE) that enables rapid prototyping and visualization of results (Gorelick et al., 2017). All these features made the GEE be used as a powerful tool for producing citrus orchards mapping for extensive case studies as fast as possible and with a minimum computational burden. Another important point is the internet bandwidth which intensely affects the GEE's performance. Iran as a developing country uses the fourth-generation (4G) internet and this has a devastating effect on the computational time when the study area increases. The fifth generation (5G) of the internet opens bright horizons for RS researchers working on cloud-based platforms.

Analysis of thematic change maps from 2016 to 2019
The RS-based method of this paper revealed a 17% total increase (about 632.1 ha) in the area of citrus orchards from 2016 to 2019. The majority of the changes have occurred in Juybar's central and southern areas. Besides the total changes, 754.2 ha of citrus orchards were converted into non-citrus LULC. The main reasons for this significant change can be many different factors, including poor economic situation, the low price of citrus in the last few years, and lack of governmental support, which made the farmers go for other alternatives such as turning citrus orchards into paddy fields.

Highlighting the research perspectives and the open research questions
In line with our research perspective, the concepts of the proposed method can be similarly used to extract other agricultural products, for example, as the most prevalent types of crops grown in the north of Iran. Therefore, in the rice paddies extraction task, a similar strategy can be used based on the phenology cycle of rice, defining different masks based on monthly VIs, and automatic sample selection. In the following, and in comparison with the findings of this study, we describe some of In line with our research findings and with the aim of extracting the citrus orchards, Xu et al. (2018a) have quantitatively compared the input features importance for such classifications done by RF. Their utilized features are as follows: topographic features (including elevation and slope), spectral features (including MS surface reflectance bands, NDVI, SAVI, and normalized difference moisture index (NWMI)), spatial features, and textural features. They found that the two topographic features have shown the most dominant effect which is followed by the features of the spectral domain, while textural features are the least important features in comparison with the aforementioned ones. They have also figured out that the elimination of a particular input feature can result in a 0.33-2.45% decrement in the OA of the classifications. As our findings showed that the spectral indices were adopted as our classification's efficient features, Xu et al. (2018a) proved that spectral indices have a greater influence on the OA in comparison with most spectral bands and textural indices. In light of the features used in related studies, such as those by Xu et al. (2018a), Rad et al. (2019), and Ali et al. (2020), and in accordance with our findings, it could be noted that feature selection for the classification of citrus orchards and many other crops highly relies upon the type of land covers existing in the intended region. This means that prior knowledge about the general nature of the case study affects the researchers' decision about the type of adopted features. To discuss the research question, our prior knowledge of a high percentage of urban areas, water, bare land, and the only evergreen crop existing in the case study region (citrus) has led us to employ our current features. Our findings showed that if the region is covered by pine forests or other tall trees or evergreen crops, we need supplementary data and features (especially HR DSM). In line with our research findings, the suitability of particular features for discriminating a specific crop from other crops has also been investigated and proved by Akbari et al. (2020).
In line with our research perspective, temporal filtering performed by Xu et al. (2018a) has also shown a significant increase in the OA of annual classifications and some degrees of decrease in the standard deviation, hence reducing illogical errors in land cover transitions. According to the findings of Saadat et al. (2019) as well as our findings, methods utilized in the SAR time-series can allow monitoring of rice cultivations in relation to rice paddies expanse, agricultural calendars, and the amount of cropping. They inferred that separating rice paddies from other crops is possible by using the C band of the SAR images due to  its sensitivity to rice growth and its dynamic range. Rad et al. (2019) have also reported that NDVI is a successful index for rice crop monitoring in the analysis of time-series data. These points can be helpful hints for better and more accurate class separation of crops in the north of Iran where citrus and rice have always been the main cultivated crops. Xu et al. (2018b) declared that seasonal MS images accompanied with the spectral, spatial, and topographic features provide them with good OA and KC in citrus orchard mapping for a case study of a middlesize county (they implemented a random forest classification in GEE). In this regard, our research is in accordance with the mentioned study, as in our research, the seasonal Sentinel MS images and hybrid feature set led to very high OA and KC. Furthermore, as the hybrid feature set provides us with relatively better accuracy compared with the visible or MS bands, this is in accordance with the study by Samadzadegan et al. (2017).

Perspectives for future research
Although this research is conducted in Juybar, it can be generalized to a coarser spatial scale, i.e., to the whole of Iran. Having in mind that the methodology is the same, only the input data may differ slightly. For instance, more coarse-resolution images from NOAA-AVHRR, SPOTvegetation, and MODIS may be propitious. In this study, we utilized ALOS DSM with a spatial resolution of 30 m, as the best available elevation data. To more efficiently discriminate citrus orchards from other evergreen land coverages such as lawns, some sorts of bushes, pine trees, oak trees, etc. normalized DSM (Eq. (31)) derived from HR DSM and digital terrain model (DTM), which can be as profitable as a height feature since orchard trees lie in a specified range of elevation (~ 2-3 m). For future research, to better extract citrus orchards, HR nDSMs from Light Detection and Ranging (LiDAR) point clouds or the point clouds obtained from the structure from motion (SfM) of UAV images, as elevation features, can be imported into the classification process as well, if available.
In addition, citrus orchards in terms of geometric patterns, texture, leaves density, and elevation can have a unique behavior in the radar range of the electromagnetic spectrum, so by using HR Radar images, one can have better performance in mapping citrus orchards. In areas with other kinds of evergreen trees, UAV or VHR aerial/satellite images can help to extract single citrus trees or orchards by providing highquality GLCM texture descriptors (Chen et al., 2019a(Chen et al., , 2019b.
The combination of Sentinel-2-derived VIs time series in GEE is recommended when there are continuous gaps in the time series, and more elaborated phonological retrieval methods, such as gap-filling and robust smoothing techniques, should be used instead. Despite this, sensor harmonization and temporal smoothing may improve land surface phenology retrievals, but at the cost of increasing preprocessing complexity, which impedes the implementation and fast computation in GEE. However, general errors that occurred in the classification of OA can be related to uncertainties of input data and inherent limitations of the classifier approach; some other challenges can also be mentioned. An imperative challenge may pertain to the mixed pixel problem. Because of the complicated and diverse land cover types in coastal areas, the presence of mixed pixels is a probable incident that can lead to inaccurate estimation of classified maps. In this research, we found that selecting samples in a point-wise manner decreases the risk of mixed pixels danger. To address the issue of mixed pixels in future studies, subpixel analysis can presumably be beneficial too. The deficiency of ground reference samples can be the other issue. To compensate for the deficiency, they can be adopted from HR imageries. In general, isolated parts of a specific land cover type in a long time series are the other problematic points. Addressing this issue with well-suited filters can be a study in the future as well.

Conclusions and recommendations
This research proposed a classification-based method for automatic mapping of citrus orchards in Juybar. The proposed method is based on analyzing the time-series images and using multi-modal satellite imageries. The developed framework can be used for monitoring citrus orchard dynamics. When validated using reference data, it produced a yearly land-cover classification sequence with reasonable accuracy, according to the conclusions of this study. The mapping performance was increased even more by using a multi-modal feature set. Because of the relevance of measurement, spectral characteristics were found to be more important determinants in citrus orchard identification. The possibility of citrus orchard mapping at a higher temporal frequency with RS time series was established in this study, which can be used as a beneficial reference for sustainable land-use management.
Numerical and visual assessment of the results shows the high potential of our method for extracting citrus orchards and producing citrus orchard maps. The produced maps can be included in various programs such as land use status and biodiversity conservation. In addition, the method presented in this study will be an efficient auxiliary method for automatically producing up-to-date land cover maps throughout the country which enables users (e.g., land managers and policy builders) to examine the dynamics of land cover over longer periods of time. The results also show that the GEE platform can provide an opportunity for controlling high computational needs and extracting spectral-temporal properties from large time-series images to a large spatial extent.
This method is especially recommended for areas with a variety of land cover classes and frequent cloud cover. The proposed method had a good performance for a study area like Juybar, which has different types of heterogeneous vegetation and it is usable for other areas. On a large scale, the finding of this research can help relevant agri-rural organizations and departments, especially the Ministry of Agriculture, by providing them with updated maps based on which they could make decisions given this accurate geospatial product. The policymakers, considering the mapping (and change monitoring) strategy of this paper, can have a better view of what is taking place in the world of agriculture, both in the short-time and long-term, and can take suitable actions.
Besides the benefits of these maps for the administrations, they could be also applicable on a fine scale. Establishing a database of citrus traits at the leaf, canopy, and orchard level is very important for future agriculture to help farmers with (1) site-specific orchard management, (2) planning farm management practices, (3) selection of best cultivars, and (4) scheduling precise and proper application of inputs in citrus production paradigm. As a result, increasing net returns and optimizing resource use are important goals. Such precise cropland maps would help the real estate agents have a GIS for the agriculture sector and better manage the affairs. As another applied example, it could be mentioned that nowadays, citrus orchards are usually sprayed by UAVs (mechanized spraying) and a precise map could help UAVs find a group  Table 7 is defined as follows: a: number of pixels with correct classification in both maps. b: number of cases that are correctly classified in batch-1 but incorrectly classified in batch-6. c: number of cases that are correctly classified in batch-6 but incorrectly classified in batch-1. d: number of pixels with the wrong classification in both maps. correct and incorrect: the number of correct and incorrect classified pixels for the classified maps of batch-1 and batch-6.
of citrus orchards and automatically spray them one by one. This helps the farmers and the stakeholders because if each farmer wants to spray his/her orchard individually and in traditional manners, it will be timeconsuming and costly due to renting the ground-based sprayer, hiring labor, buying much more poisons, etc.

Funding
This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Availability of data and material
Raw data were generated at Tehran University. We confirm that the data, models, or methodology used in the research are proprietary, and derived data supporting the findings of this study are available from the first author on request.

Declararion of Competing Interest
We have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Declaration of Competing Interest
None.

Data availability
Data will be made available on request.