Performance and the Optimal Integration of Sentinel-1/2 Time-Series Features for Crop Classiﬁcation in Northern Mongolia

: Accurate and early crop-type maps are essential for agricultural policy development and food production assessment at regional and national levels. This study aims to produce a crop-type map with acceptable accuracy and spatial resolution in northern Mongolia by optimizing the combination of Sentinel-1 (S1) and Sentinel-2 (S2) images with the Google Earth Engine (GEE) environment. A total of three satellite data combination scenarios are set, including S1 alone, S2 alone, and the combination of S1 and S2. In order to avoid the impact of data gaps caused by clouds on crop classiﬁcation, this study reconstructed the time series of S1 and S2 with a 10-day interval using the median composite method, linear moving interpolation, and Savitzky–Golay (SG) ﬁlter. Our results indicated that crop-type classiﬁcation accuracy increased with the increase in data length to all three data combination scenarios. S2 alone has higher accuracy than S1 alone and the combination of S1 and S2. The crop-type map with the highest accuracy was generated using S2 data from 150 days of the year (DOY) (11 May) to 260 DOY (18 September). The OA and kappa were 0.93 and 0.78, respectively, and the F1-score for spring wheat and rapeseed were 0.96 and 0.80, respectively. The classiﬁcation accuracy of the crop increased rapidly from 210 DOY (end of July) to 260 DOY (August to mid-September), and then it remained stable after 260 DOY. Based on our analysis, we ﬁlled the gap of the crop-type map with 10 m spatial resolution in northern Mongolia, revealing the best satellite combination and the best period for crop-type classiﬁcation, which can beneﬁt the achievement of sustainable development goals 2 (SDGs2).


Introduction
Agriculture is a vital sector in the economy of Mongolia, which accounted for 10.8% of the Gross Domestic Product (GDP) in 2018 [1]. The northern part of Mongolia is a agriculturally dominated region, and is vulnerable to weather conditions since about 80% of the cultivated land is rain-fed. As the principal cereal crop of Mongolia, wheat production has drastically declined due to policy changes, investment reduction, and successive severe droughts since the mid-1990s [2]. The accurate and timely mapping of early crop-type identification is also critical for Mongolia's crop production. Therefore, it is necessary to determine the most suitable time windows for crop type classification in Mongolia.
This study takes northern Mongolia as the study area. This study aims to design a crop classification framework and fill the crop-type map gap with acceptable spatial resolution. In addition, this research expects to answer the following scientific questions: (1) Which satellite combination is best for crop-type classification in northern Mongolia? (2) What are the suitable time windows for crop classification in northern Mongolia? (3) What is the most effective and sensitive metric for crop classification?

Study Area
The provinces of Darkhan and Selenge are located in the north-central part of Mongolia ( Figure 1). The study area has the most favorable natural conditions for crop growth and produced approximately 42% of Mongolia's national wheat in 2018 [1]. Almost all crops in the study area are rainfed crops, which are very sensitive to changes in weather conditions. The study area is dominated by a temperate continental climate with a dry and cold winter and a humid and warm summer. According to meteorological data obtained from the Information and Research Institute of Meteorology, Hydrology, and Environment (IRIMHE) in Mongolia, annual precipitation in 2018 was 313 mm, with approximately 71% (224 mm) of the annual precipitation occurring in June-August. The average annual temperature is 2 • C, being −22.5 • C in January and 20 • C in July ( Figure 2). Due to the constraints of the accumulated temperature, crops in the area are dominated by cold-tolerant and early maturing spring wheat varieties and rapeseed. Spring wheat and rapeseed are the major crops of the study area. Figure 3 shows the phenology of spring wheat and rapeseed. In general, spring wheat is sown in early May, the critical growing period is from June to August, and the harvest period is from late September to early October. Rapeseed has similar phenology to spring wheat, which is sown in May and harvested in September. The length of the growing period of crops in the region is 90 to 140 days.

Collection of Crop Reference Data
Data were collected between the 6 to 26 September 2018, close to the harvest season of the crop (Figure 4). Information collected included geographic location (longitude and latitude), crop-type name, and photos. Samples were collected by the GPS-Video-GIS application (GVG) [53], an all-in-one system of location, photo, and attribute tagging providing quick access to photo information with geolocation and crop-type information. The GVG software is a smartphone application that can be downloaded for free from: https://gvgserver.cropwatch.com.cn/download, accessed on 5 April 2022. A total of 721 samples were collected, including 508 spring wheat samples, 76 rapeseed samples, and 137 fallow field samples during a 3000 km-long field trip. After completing the field survey, we checked and corrected the position of each sample using high-resolution imagery. Finally, the samples were split into training and validation samples in a ratio of 60:40 for crop classification and accuracy evaluation.

Satellite Data Processing
In this study, the S1 synthetic aperture radar (SAR) ground range detected (GRD) and S2 top-of-atmosphere (TOA) reflectance images (Level-1C) were used for crop type identification. S1 and S2 were launched by the European Space Agency (ESA). GEE has archived all S1 and S2 data for free use. The S2 optical images cannot provide helpful information in cloudy and rainy weather, while the SAR satellite of S1 can obtain usable images in cloud-covered conditions. Therefore, we used a combination of S1 and S2 images to find the best combination of features to classify crop types and extract crop types in the study area.

Sentinel-2 Data
S2 is a 14-band multispectral acceptable spatial resolution optical image with 4 bands of 10 m (blue, green, red, near-infrared), 6 bands of 20 m (red edge 1, red edge 2, red edge 3, red edge 4, swir1, swir2) and 4 bands of 60 m (aerosol, water vapor, cirrus, cloud cover) spatial resolution. In this study, 20 m reflectances were resampled to 10 m by the nearest-neighbor approach [54]. S2 consists of 2 satellites (2A and 2B) with ten days of revisit time for each satellite and an interception/overlapping period of 5 days for both satellites. S2 TOA Level-1C images were acquired from 11 May to 31 October 2018, followed by cloud-free time-series reconstruction processing. S2 TOA Level-1C data are widely used in crop-type classification and have a minimal impact on crop classification results, such as paddy [27], maize [29], winter wheat [30], and mixed crop types' identification [31].
(1) Cloud free The cloud-masking band of S2 provides information on cloud cover, but studies have shown that it is not as effective when used to remove pixels contaminated by cloud cover [40]. In this study, a modified Luo-Trishchenko-Khlopenkov (LTK) scene identification algorithm [55] designed for Moderate Resolution Imaging Spectroradiometer (MODIS) data was used to clear cloud and shadow of Landsat [55] and S2 [40] data. You et al. [31,40] explained the processing principle of the modified LTK for cloud detection of S2. It uses aerosol, blue, green, and red bands of S2, normalized difference moisture index (NDMI), and normalized difference snow accumulation index (NDSI) to calculate cloud volume and detect cloud layers. This study also employed modified LTK to remove the pixels of S2 contaminated by clouds.
(2) Time-series reconstruction of metric Swir1 and Swir2 reflectance bands have a good ability to separate crops [25,31]. Enhanced vegetation index (EVI) [56], green chlorophyll vegetation index (GCVI) [27,57] normalized differential senescent vegetation index (NDSVI) [58], normalized vegetation index (NDVI) [59], red-edge NDVI (RENDVI) [35] and red-edge position (REP) [35] indices were used as features to extract crop types from the study sites. NDVI is an index widely used to identify the distribution of crop types; EVI has a strong ability to overcome the saturation phenomenon during the growing season; RENDVI and REP are red-edge-based indices, both sensitive to crop conditions; and the GCVI index proved to have a strong ability to separate crops. The formulae for NDVI, EVI, GCVI, RENDVI, REP, and NDSVI are summarized below. (1) where Nir, Red, Blue, Green, Red1, Red2, Red3, and Swir1 are the reflection bands of S2. The above features are discontinuous due to the cloud cover, increasing crop differentiation difficulty. In this study, the median composite method, linear moving interpolation [38,54] and Savitzky-Golay (SG) filter [39], were utilized to reconstruct the continuous S2 features within a 10-day interval.
In step 1, the time series with a 10-day interval of S2 features are generated by the median composite method, and a total of 18 images for each feature were generated by the median composite method from May to October 2018. The median composite method is easily applied and proved to be superior to the mean composite method by many previous studies [40,60,61].
In step 2, the gap existing in the S2 composite time series is filled by the linear temporal moving interpolation method [31,38,40]. After analyzing the length of the data gaps, we found that the maximum missing data length is 30 days; therefore, the half moving window is set to 3 (equals 30 days, Equation (7)): where Index i is the index existing gap, and Index k and Index j the good observation before and after Index i . d k , d i , and d j and day of a year of kth, ith and jth.
Step 3: After filling the gaps in the S2 composite time series, the noise of the S2 time series generated by the linear shift interpolation method is smoothed with an SG filter with a window at 3 × 3. In order to avoid the memory computational error, the smoothed time series is first imported into assets of GEE and then exported to the javascript code.

Sentinel-1 Data
We collected all S1 data from 11 May to 31 October 2018 to determine the crop type for this study. S1 in C band imagery (frequency = 5.4 GHz) globally had a 12 or 6 day revisit cycle depending on 1A and 1B imagery availability. The IW modes of S1 provided dual-polarization with vertical transmit and vertical receive (VV), vertical transmit, and horizontal receive (VH). VV and VH of S1 images with 10 m are widely used in the crop-type classification [34]. S1 in GEE was proceeded by Sentinel-1 Toolbox using the thermal noise removal, radiometric calibration, and terrain correction. The edge of S1 is removed using a 'connected Components' method in GEE, which mask groups of neighboring pixels processed with values lower than −25 dB in the VV polarization [34,62]. The noise of VV and VH is filtered by the Refined Lee filter [63]. The backscattering coefficients are converted to decibels (dB) via log scaling (dB = 10 × log10(X)) in GEE. Some studies indicated that the ratio of VH and VV (here defined as VH/VV) is a good index to distinguish crops, but it should be calculated in the natural domain [34]; therefore, this study converts the backscattering coefficients of VH and VV from dB value into the natural value by X = 10 × pow (dB/10), and then calculates VH/VV value. Finally, we used the same approaches introduced in Section 2.3.1 to generate time series of VV, VH, and VH/VV with 10-day intervals. The smoothed VV, VH, and VH/VV time series are also imported into assets of GEE and then exported to the javascript code.

Methodology
The crop-type identification framework is illustrated in Figure 5. This framework consists of five modules: satellite processing and indices' preparation, metrics' time series preparation, scenario selection, reference samples' preparation, classification, and accuracy assessment.

Satellites' Processing and Indices' Preparation
The primary role of this module is to process S1 and S2 images and provide them for metrics' generation. Modified LTK methods eliminate the cloudy pixels in S2 image collection, and then the Swri1 and Swir2 reflectance bands are selected for metrics' generation. Vegetation indices, including EVI, GCVI, NDSVI, NDVI, RENDVI, and REP, are generated by equations introduced in Section 2.3.1. Because other reflectance bands are used to calculate vegetation indices, they are excluded in crop-type classification to avoid information redundancy. The edge of S1 with IW mode was removed by 'connected Components', and then the noise was suppressed by Refined Lee Filter, calculating the VH/VV in the natural domain. Finally, VH, VV bands in the log domain, and VH/VV, were prepared for metrics generation.

Metrics' Time-Series Preparation
The function of this module is to provide metrics for crop classification in different cases. Eight S2 features and three S1 features are processed by median synthesis, linear shift interpolation, and SG smoothing methods to generate metrics with an interval of 10 days. The method of time-series reconstruction is described in Section 2.3.1. The method had successfully been applied to crop-type identification in northeast China [31,40]. This module generated 198 features (11 × 18) between 11 May and 31 October 2018.
The M-index [64,65] can determine the ability of a specific characteristic time series to distinguish between two categories. The formula of the M-index can be described as

Scenario Selection
Three scenarios were set in the module, namely crop classification with S1 alone, crop classification with S2 alone, and crop classification with S1/S2 together. The performance of the three scenarios will be compared to select the best combination and the earliest and optimal period for crop-type classification in the study areas.

Reference Samples' Preparation
The role of this module is to provide training samples for the classifier and to provide validated samples for result accuracy assessment. In this module, the samples of spring wheat, rapeseed, and fallow land are divided into training and validating groups at 60:40. We had 305 samples (spring wheat: 240, rapeseed: 51, fallow land: 14) for training and 208 samples (spring wheat: 168, rapeseed: 32, fallow land: 8) for validation.

Classifier: Random Forest
The random forest (RF) classifier was used to classify the crop types. RF builds multiple trees based on random bootstrapped training data samples [66]. Each tree contains a subset of training data through replacement (a bagging approach), and the nodes are split using the best split variables. Typically, about two-thirds of the samples (in-bag samples) are used for the training tree, while the remaining one-third (out-of-the-bag samples) are used for internal cross-validation techniques to estimate the performance of the RF algorithm [66]. The error estimating process of the RF model is called the out-of-bag (OOB) error. Each decision tree is generated independently without any pruning, and each node is split using a user-defined number of features selected randomly [66]. As many previous studies have proven the excellent performance of RF [42,67], we employ the GEE platform to train and apply an RF classifier to map crop types due to its high performance. We extracted all S1 and S2 features for analysis of crop classification between May and October. To reduce the effect of non-cropping information on the classifier's crop classification effectiveness, we removed non-cropping satellite imagery with a cropping mask to improve crop classification accuracy [68,69]. The crop cover mask is provided by the IRIMHE in Mongolia [70].

Accuracy Assessment
This study used the RF classifier and the collected field crop-type data for crop-type classification. For each class, 60% of the total points were randomly selected for training, and the remaining 40% of field survey points were used to assess the accuracy. The confusion metrics were used to assess crop-type classification accuracy, including user accuracy (UA), producer accuracy (PA), overall accuracy (OA), Kappa coefficient, and F1-score. These accuracy calculations were presented in the equations below (more details in R. G. Congalton) [71].
where X ii is the number of observations in row i and column i, X i+ and X +i are the marginal totals of row i and column j, and N is the total number of observations. Figure 6 shows the 10-day intervals of EVI, GCVI, NDSVI, NDVI, RENDVI, REP, SWIR1, SWIR2, VH, VV, and VH/VV curves for spring wheat, rapeseed, and fallow land. EVI, GCVI, NDVI, and VH/VV curves indicate that spring wheat and rapeseed reached their growth peaks at day 210 of the year (DOY) (late July) and 230 DOY (late August). The difference in NDVI between rapeseed and spring wheat indicates that rapeseed was greener than spring wheat in the study area. The differences in EVI, GCVI, NDVI, and VH/V indicate that they could separate spring wheat and rapeseed, while NDSVI and REP are difficult to separate due to similar curves. To estimate the separation between spring wheat and rapeseed, we calculated the M values of the 12 indices for spring wheat, rapeseed, and fallow land (Figure 7). The M value between spring wheat and rapeseed was much lower than 1, indicating that separating spring wheat and rapeseed with a single index was challenging. Compared with spring wheat and rapeseed, the M values between fallow and non-fallow lands were relatively high, and the separation between them was better than that between spring wheat and rapeseed. Figure 7. The M value between spring wheat, rapeseed, and fallow land for 12 features (W, R, and F represent spring wheat, rapeseed, and fallow land, respectively). Figure 8 shows the OA, kappa, and F1-score variation for spring wheat and rapeseed. The OA, kappa, and F1-score increased with the extension of monitoring time length. S2based crop classification performed best in all three scenarios with the highest OA, kappa, and F1-score. The results of OA, kappa, and F1-score values are shown in Appendix A (Tables A1 and A2). The combination of S1 and S2 was superior to S1 alone for crop classification. A study in Heilongjiang province also found that the accuracy of using S1 only was much lower than the results using S2 and S1 [40]. The change in OA based on S2 indicates that OA increased rapidly from day 190 to day 260, and the largest increase rate (0.029) occurred from 210 DOY to 220 DOY. The same trend was observed in the change in the kappa coefficient. Spring wheat was the dominant crop in the study area. The F1-score of spring wheat was significantly higher than that of rapeseed, and it increased slowly from 0.89 to 0.96; the F1-score remained unchanged after 240 DOY. The F1-score of rapeseed maintained a rapid increase from 140 DOY (F1 = 0.22) to 260 DOY (F1 = 0.80). From 210 DOY (end of July) to 260 DOY (August to mid-September), the F1-score of rapeseed maintained a rapid growth. The curves of Figure 6 indicate that the difference in features increased from 210 DOY to 260 DOY, which can explain the changes in OA, Kappa, and F1-score.

Accuracy Assessment
OA, kappa, and F1-score indicators show that the best classification results could be achieved between DOY 130 and DOY 260 using S2 data ( Figure 9). During this period, OA and kappa were 0.93 and 0.78, respectively; the F1-scores for spring wheat and rapeseed were 0.96 and 0.80, respectively (Table 1). PA was 0.98, 0.75, and 0.75 for spring wheat, rapeseed, and fallow. UA was 0.94, 0.86, and 1.00 for spring wheat, rapeseed, and fallow. The largest mixed classification occurred in rapeseed and fallow, one-quarter of which were misclassified as spring wheat.

Classification Results
We produced and assessed the crop-type map using only the optical images (S2), only SAR (S1), and the combination of SAR (S1) and optical (S2) images in northern Mongolia. The highest accuracy with the crop-type map using S2 satellite data for Selenge and Darkhan Provinces in Mongolia is shown in Figure 10. According to our designed method, the details of crop-type distribution were better extracted ( Figure 11). Spring wheat is the main crop in this region. The classification results are divided into three classes, i.e., wheat, rapeseed, and fallow, and the highest OA of the three classes was 93%, respectively. From the analysis, 86.6% of the total cropland area was spring wheat, rapeseed accounted for only 5.7% of the cropland area, and fallow land accounted for 7.7% of the cropland in 2018.

Discussion
This study demonstrated 10 m crop-type classification using S1 and S2 data using the RF algorithm in the GEE cloud computing platform in northern Mongolia in 2018. The produced high-resolution crop-type map is essential for sustainable agricultural development, regional food security, crop planting areas, and spatial distribution of crops. In Mongolia, a previous study mapped only croplands at the country level, with a 30 m spatial resolution [70]. However, the classification of the specific crop types is still challenging. This paper proposed and presented an RF classifier approach based on the optimal data combination of S1 and S2. Finally, we generated the crop types map at 10 m spatial resolution in northern Mongolia using S2 alone, which generated the first 10 m spatial resolution crop types in this area.

Classification Accuracy Compared with Others
The classification results using S2 satellite data have the highest performance. The PA and UA were 0.94 and 0.98, 0.86 and 0.75, and 1 and 0.75 for spring wheat, rapeseed, and fallow land, respectively. The OA and kappa of crop-type classification reached 0.93 and 0.78, respectively. There are many studies on crop classification. For instance, R. Sonobe et al. [72] used 82 vegetation metrics of S2 to identify crop spatial distribution, and the OA obtained by random forest was from 90.2% to 92.2%. R. Saini and S. K. Ghosh et al. [73] have used the S2 optical images to classify crops. The OAs of the classification achieved by RF and SVM using S2 imagery were 84.22% and 81.85%, respectively. Jiang et al. [74] classified crop mapping in China using S2-derived indicators (NDVI, REP, and NDRI). The final results show that the OA of the classified crop map was within 0.94. Immitzer et al. [75] performed crop-type classification using S2 data. The result showed overall accuracies ranging between 65% for tree species and 76% for crop types in Germany. Furthermore, numerous crop type classification and mapping studies are presented based on the Landsat and S2 optical images [32,74,76]. Song and Wang [77] identified the winter wheat planting area using S1 data in North China Plain. The OA reached 84%, and the Kappa coefficient was 0.77, lower than our result. Compared to classification accuracy reported by previous studies, it is evident that our results are acceptable.

Impact of SAR Information on Crop-Type Classification
Several studies utilized S1 in crop classification to overcome S2 limitations on cloudy days [42][43][44][45][46][47]. However, the OA of classification using only S1 has a lower accuracy than optical imagery-based classifications. In our study, the classification accuracy was only SARbased and significantly lower than that which was optical-based. Many previous studies have shown that the combination of optical and SAR data provides better discrimination of crops than optical data alone [32,78,79]. For example, K.Van Tricht et al. [32] used joint optical and SAR imagery to classify crop maps in Belgium. They found that the result based on optical and SAR data together outperformed optical or SAR data alone. However, our study achieved the best crop classification using only S2 features. To further validate our results, this study compared the crop classification accuracy of S2 with different combinations of SAR information (VV, VH, and VH/VV) using the best classification period from 11 May to 18 September 2018 (Table 2). Among the above combinations, the best crop classification results came from the S2 and VH/VV combinations with OA of 0.92, Kappa of 0.74, and F1-score of 0.95 for spring wheat F1-score of 0.76 for rapeseed. The worst crop results came from the S2, VV, VH, and VH/V combinations with OA of 0.91, Kappa of 0.69, and F1-score of 0.94 for spring wheat F1-score of 0.70 for rapeseed. This result shows that the combined use of optical and SAR data does not always yield the best classification results. Our findings are supported by the study of You et al. [40], who found a minimal contribution of S1 to the separation of maize and soybean in Heilongjiang Province, China. In short, S2-based is the best data combination among three data integration scenarios for crop-type classification in northern Mongolia. Our results with those of other studies indicate that crop-type accuracy varies with study area, and the best data combinations should be compared to each other to achieve the best crop-type classification accuracy.

Compare to Percentiles Composite of Time-Series Reconstruction
We reconstructed continuous time series of S2 by median synthesis, linear moving interpolation, and SG filters. These time series describe changes in the growth process of spring wheat and rapeseed and capture changes in phenology. Percentile composite is widely used in metric synthesis [27,29] as it is less sensitive to the length and continuity of satellite data [80]. In this study, we used percentile synthesis (percentiles equal 5,25,50,75,95) to generate metrics for S2 from 130 DOY (11 May) to 260 DOY (18 September) 2018. The OA and kappa were 0.90 and 0.64 based on the percentile approach, respectively, and OA was significantly lower than the 10-day-intervals time-series reconstruction results. This suggested that the 10-day-interval time-series reconstruction performs well in classifying crop types in Mongolia. This can be attributed to the dry climate of Mongolia, which means that the amount of available data is high, and there are fewer missing data values, allowing for better reconstruction of continuous time series.

Shortcomings
Although good results were achieved in crop-type classification, the F1-score for rapeseed was only 0.80, and about one-quarter of the rapeseed samples were misclassified as spring wheat, meaning that the accuracy of rapeseed classification needs to be improved in the future. Misclassification of wheat and rapeseed is a challenge that needs to be addressed. In the future, a better indicator should be devised to improve the separability of spring wheat and rapeseed. The time-series reconstruction method used in this study achieved good performance in crop-type classification due to the dry climate in the study area. However, its performance in humid regions should be tested robustly.

Conclusions
The current study fills a gap in high-spatial-resolution crop types in northern Mongolia and could provide strong support for yield estimation of spring wheat. This study used cloud computing technology to compare crop-type classification using three scenarios (S1 alone, S2 alone, S1 + S2) in northern Mongolia in 2018. The study results show that the S2-based classification results outperformed S1 + S2-and S1-based ones. The classification results improved with increasing monitoring time length, with 130 DOY (11 May) to 260 DOY (18 September) being the best period for crop-type classification, implying that good results can be achieved as early as one month before harvest. Thus, from our results, we recommend using S2 data between 130 DOY and 260 DOY for crop-type classification mapping in northern Mongolia. The crop classification effect improved rapidly from 210 DOY (end of July) to 260 DOY (August to mid-September). Our results also indicate that differences between optical and SAR satellite combinations need to be compared to select the best satellite combination in the crop classification process. Furthermore, our result also indicate that the 10-day-interval data integration outperformed the percentile composites. In short, the framework demonstration in this study can prioritize mitigation plans to support national food security in Mongolia.  Data Availability Statement: All relevant data are shown in the paper.

Acknowledgments:
We are very grateful to the Information and Research Institute of Meteorology, Hydrology, and Environment (IRIMHE) of Mongolia. Special thanks to G. Sarantuya for providing human resources and required instruments for fieldwork. We are very grateful to the CropWatch team and United Nations Economic and Social Commission for Asia and the Pacific (UNESCAP) for supporting this study. Finally, we would like to express our deep appreciation to the anonymous reviewers for reviewing the manuscript and providing comments to improve the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1. Overall accuracy (OA) and Kappa index values for the different classification scenarios.