Two-Stepwise Hierarchical Adaptive Threshold Method for Automatic Rapeseed Mapping over Jiangsu Using Harmonized Landsat/Sentinel-2

Rapeseed distribution mapping is a crucial issue for food and oil security, entertainment, and tourism development. Previous studies have used various remote sensing approaches to map rapeseed. However, the time-consuming and labor-intensive sample data used in these supervised classification methods greatly limit the development of large-scale mapping in rapeseed studies. Regarding threshold methods, some empirical thresholding methods still need sample data to select the optimal threshold value, and their accuracies decrease when a fixed threshold is applied in complex and diverse environments. This study first developed the Normalized Difference Rapeseed Index (NDRI), defined as the difference in green and short-wave infrared bands divided by their sum, to find a suitable feature to distinguish rapeseed from other types of crops. Next, a twostepwise hierarchical adaptive thresholding (THAT) algorithm requiring no training data was used to automatically extract rapeseed in Xinghua. Finally, two adaptive thresholding methods of the standalone Otsu and Otsu with Canny Edge Detection (OCED) were used to extract rapeseed across Jiangsu province. The results show that (1) NDRI can separate rapeseed from other vegetation well; (2) the OCED-THAT method can accurately map rapeseed in Jiangsu with an overall accuracy (OA) of 0.9559 and a Kappa coefficient of 0.8569, and it performed better than the Otsu-THAT method; (3) the OCED-THAT method had a lower but acceptable accuracy than the Random Forest method (OA = 0.9806 and Kappa = 0.9391). This study indicates that the THAT model is a promising automatic method for mapping rapeseed.


Introduction
Detailed and fine rapeseed distribution mapping is of great importance for three reasons. Firstly, rapeseed map products are not only significant for food and oil security (rapeseed is an oily crop providing edible oil for the increasing population [1,2]), but they are also meaningful for animal husbandry (rapeseed contains rich protein for animal feed [3]). Secondly, the distribution information of rapeseed is pivotal for further research, such as crop yield prediction and the assessment of climate change effects on agriculture, promoting the development of digital agriculture [4]. Thirdly, the beautiful yellow flowers of rapeseed fields provide resources for entertainment and tourism [5]; thus, rapeseed maps offer lots of details for decision makers to plan for tourism development.
Remote sensing techniques have been widely used in crop mapping because they are timely, labor-saving, and inexpensive in contrast with traditional field surveys [6,7].
Previous studies showed that the mapping of rapeseed attracts less attention than that of its counterparts (e.g., winter wheat and paddy rice) [2]. The challenges of rapeseed mapping include the following: similarity in spectral response between rapeseed and other vegetation [2], spectral mixing pixels due to patchy rapeseed fields [8], cloud contamination reducing effect observation times [5], and expensive training samples limiting large-scale studies [5,9,10].
Existing studies concerning rapeseed mapping can be divided into two categories. The first category comprises studies employing ensemble classification techniques, such as Random Forest (RF) and Support Vector Machine (SVM), using images taken during the flowering period [11][12][13][14]. For example, Meng et al. (2020) used all the cloud-free Sentinel-2 images to select the optimal temporal phases for winter wheat and rapeseed extraction in Zhongxiang, China [11]. This method has a high accuracy but depends on large amounts of training samples.
The second category comprises studies employing phenological-based threshold segmentation [9,15,16]. This method has been receiving increasing attention for automatic crop mapping because it has a promising accuracy with an empirical threshold [9,[17][18][19]. Nevertheless, a fixed threshold requires some samples to choose an optimal threshold for the purpose of improving accuracy, which, in some studies, was not enough to meet the requirement of zero training samples [9,[16][17][18], and in other studies, it was challenging to apply it in other areas because the accuracy decreased when the environment changed [20,21]. Therefore, in some cases, the threshold was manually intervened with to achieve better performance [22]. For instance, Ashourloo (2019) developed the Canola Index (CI) concerning near-infrared (NIR), red, and green bands for automatic rapeseed mapping, while its threshold was selected using an empirical method, which needed sample data to choose an optimal threshold when extending to testing sites [9]. Some adaptive threshold methods, such as the Otsu thresholding method originally used in image segmentation [23], have been introduced in other research fields. The Otsu thresholding method has become popular in remote sensing areas. For example, Zhang (2021) and Singh (2020) employed the Otsu with Canny Edge Detection (OCED) thresholding method to extract aquatic vegetation from water [20,22]. These methods are usually used to extract open water from land, aquatic vegetation from water, and glaciers from non-glaciers because the images can be simply divided into two parts: the object and the background [20,22]. Furthermore, objects and backgrounds have huge differences in spectral response, so a threshold can separate objects and backgrounds easily [22]. To the best of our knowledge, the adaptive thresholding method has not been applied in crop mapping to date, because crops share similarities in spectra, making it a challenge to determine a suitable threshold.
Regarding classification features, spectral and phenological characteristics form the hypothetical bases to discriminate crops [6,11,24]. The Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), and Soil-Adjusted Vegetation Index (SAVI) are major features in most staple crop-mapping studies [2,25]. Because the canopy of rapeseed differs greatly from that of other crops during the flowering stage, new taxonomic features are often mined for rapeseed in studies of rapeseed extraction. Wang (2018) developed the Ratio Oilseed Rape Colorimetric Index (RRCI) based on HSV transformation according to the principle that the color of rapeseed landscapes differs from that of other landscapes during the flowering period, and then a decision tree was employed for automatic rapeseed extraction using a uniform threshold in Hubei province [15]. D'Andrimont (2020) used the Normalized Difference Yellow Index (NDYI) concerning green and blue bands to extract rapeseed fields based on RF and then tracked the flowering date [10]. Veloso (2017) explored the temporal profile of NDVI and the radar backscatter of various crops, including rapeseed, to evaluate phenological characteristics for potential applications [26]. The results of their study suggest that radar and visible-near infrared are the major features in crop mapping. However, there is a significant lack of short-wave infrared (SWIR) bands in phenological features, with focus usually placed on visible-near infrared (vis-NIR) bands. The SWIR band is sensitive to water, and it can improve our understanding of vegetation moisture content; thus, it may provide more information to identify specific crops [27,28].
The objective of this study was to automatically map rapeseed using a two-stepwise hierarchical adaptive threshold method. This was undertaken by developing an adaptivethreshold-based algorithm, which includes two key points: the optimal feature and the adaptive threshold. Accordingly, this study aimed to (1) determine the optimal rapeseed discrimination feature; (2) find a suitable adaptive threshold method; and (3) develop an automatic rapeseed mapping method based on the feature and the adaptive threshold method mentioned above and apply this method in Jiangsu Province on the Google Earth Engine (GEE) platform.

Study Area
As shown in Figure 1, Xinghua County and Jiangsu Province were chosen as the study sites. Xinghua County was the training area because the field data used to develop the method were collected in Xinghua County, and Jiangsu Province was the testing area in order to evaluate the robustness of the proposed approach in a large-scale region with a more complicated environment.   [29], eastern China, covering an area of 102,600 km 2 [30]. Xinghua County has an area of 2393.35 km 2 and lies in the middle of Jiangsu Province, sharing a climatic environment and a crop planting structure with Jiangsu Province.
The climate type is a transition zone between subtropical and temperate monsoon climate zones, where the annual average temperature ranges from 13.6 to 16.1 • C, and the average annual precipitation is 1000 mm, which provides a suitable environment for agricultural development [31]. The summer crops in the study area are paddy rice, corn, and soybean; the winter crops are winter wheat and rapeseed. According to the situ survey, rapeseed is planted in a fragmented manner, mostly distributed on the ridge and polder. The phenological calendar of rapeseed presented in Table 1 was collected from the Meteorological Bureau of Jiangsu Province.  II III I  II III I  II III I  II III I  II III I  II III I  II III I  II III I  II  Phase Seedling: blue; tillering: green; overwintering: gray; bolting: orange; flowering: yellow; pod filling: purple; harvesting: pink.

Sample Data Acquisition
The sample datasets were collected by conducting a field survey in Xinghua County in April and by visually interpreting high-resolution images, including those from Sentinel-2 and Landsat-8, on the GEE platform [32]. The samples were divided into rapeseed, wheat, trees, shrublands, river, aqua-farm ponds, built-up, greenhouses, and solar panels. These landscapes were split into three groups: rapeseed formed a group; wheat, trees, and shrublands were grouped into other vegetation; and river, aqua-farm ponds, built-up, greenhouses, and solar panels were grouped into non-vegetation. The number of samples is shown in Table 2, and their distribution is shown in Figure 2. Notably, all the samples in Jiangsu Province were only used for validation, while the samples in Xinghua County were split into training and validation samples. The former was used to develop the rapeseed automatic method by analyzing spectral temporal behaviors; the latter was used for validation.

Remote Sensing Data
Remote sensing images were retrieved from the GEE platform. The GEE cloudcomputed platform provides large remotely sensed images and vast computing power. Considering the access and the consistency of Sentinel-2, Landsat-7, and Landsat-8, several bands, including red, green, blue, NIR, SWIR1, and SWIR2, of these three sensors were selected [33]. In particular, the B8A band of Sentinel-2 was selected for the NIR band because it is a better match for Landsat imagery [4]. In this study, images taken during the rapeseed growth cycle were used for spectral analysis, and images taken during the rapeseed flowering period were used for rapeseed mapping.

Data Preprocessing
Images in the GEE platform are well preprocessed with radiometric calibration and registration [16,33]. Thus, the image preparation process included cloud removal, harmonization, and temporal aggregation.
Cloud removal was the first preprocessing step. Jiangsu Province is always cloudy during the rapeseed life cycle. Due to cloud and shadow contamination, it is a challenge to obtain cloud-free images during a specific period [34]. To reduce the effect of clouds, different cloud removal methods were used for different datasets [35]. The Quality Assessment (QA) band was used to identify clouds, and then the mask function was applied to remove pixels whose cloud cover was less than 50% in the Landsat series [16,36,37], while the Sentinel-2 probability dataset was used to detect clouds in Sentinel-2 imagery [38].
The second preprocessing step was harmonization. After cloud removal, there were some missing data in the imagery. Before using temporal aggregation to fill these missing pixels, harmonization needs to be carried out to reduce the differences in resolution, spectral bandwidth, spectral response function, and the field of view of the different sensors because Jiangsu Province is such a cloudy region that it is a challenge to obtain a completely clear image during the rapeseed growing period with standalone sensor imagery. Here, the widely used linear regression was used to reduce the spectral difference between sensors [39][40][41][42]. The linear regression formula is shown in Table 3 [42]. ETM+ and OLI imagery were harmonized to the reference of MSI data [40]. Thereafter, these imageries were compiled into one image dataset called the Harmonized Landsat and Sentinel-2 (HLS) dataset [43]. Temporal aggregation was the final preprocessing step. It was carried out using the median algorithm on the GEE platform. As a result, a median value layer stack of every pixel in a time-series image dataset was calculated and used to fill in the missing pixels caused by cloud removal [4,13]. Moreover, the median value layer stack prevented oversaturation and lower values; thus, the time-series data could be smoothed [16,44]. The time window was not the same for every temporal aggregation application. A shorter one may lead to noise; a longer one may cause the loss of critical phenological information. There are options for time intervals for temporal aggregation when performing diverse functions. Here, temporal aggregation was used three times, and each time had a different time window. First, the HLS dataset was daily aggregated to analyze the temporal behaviors of all the selected bands so that the high temporal resolution could keep tracking the entire trajectory of the landscapes. Moreover, daily aggregation allows for the median value to be calculated in overlapping images for the same day, paving the way for analysis. Second, the time-series Normalized Difference Rapeseed Index (NDRI) and NDVI datasets were monthly aggregated because the monthly temporal curve was enough to capture key information. Third, the flowering phase (from 1 March to 30 May) phenological window was used to produce a seamless cloud-free composite image covering Jiangsu Province. The number of images available is shown in Table 4. All images listed in Table 4 were used in spectral temporal analysis; images taken during March and April were used for rapeseed mapping.

Two-Stepwise Hierarchical Adaptive-Threshold-Based Rapeseed Mapping Model
The two-stepwise hierarchical adaptive-threshold-based method (THAT) was introduced to automatically extract rapeseed in Xinghua County, and it was tested in both Xinghua and Jiangsu. The workflow is shown in Figure 3. First, cloud detection, imagery harmonization, and temporal aggregation were applied to prepare an HLS dataset. Thereafter, the HLS dataset was smoothed using the Savitzky-Golay (S-G) filter, and then it was used to analyze the temporal characteristics of the different bands of the various landcovers. A new vegetation index named the Normalized Difference Rapeseed Index (NDRI) was developed according to the temporal spectral analysis. A two-stage hierarchical classification was proposed to first extract vegetation and then rapeseed [20,45]. Using the NDVI and the Otsu thresholding method (NDVI-Otsu), the non-vegetation landcover was masked because this study focused only on vegetation. Next, NDRI was combined with the adaptive threshold to automatically extract rapeseed. When using NDRI and the Otsu thresholding method, it is called NDRI-Otsu; when using NDRI and the OCED thresholding method, it is called NDRI-OCED. Finally, validation was used to evaluate the potential of the developed method.

Construction of Normalized Difference Rapeseed Index Using Temporal Spectral Analysis
To select optimal features, temporal spectral analysis was conducted to understand the temporal spectral trajectory of different landcovers during the rapeseed growth cycle [16,26,46]. The spectral reflectance of vegetation is seasonal, expressing phenological details at the different growth stages that can be used to extract certain crops from other crops that share a similar spectrum [7,13]. Before temporal analysis, the S-G filter was adopted to smooth the time-series spectral bands [34]. The optimal bands were selected based on one premise: if there are maximum and minimum bands in a rapeseed field in a certain period, these bands are the optimal bands that can be used to determine the features that distinguish rapeseed from other vegetation [17].
As a further step, the Jeffries-Matusita (JM) distance was calculated to quantitatively compare the indices commonly used in rapeseed mapping [5]. Here, the sample data of trees, shrublands, and wheat were grouped into the other category, and then the JM distance was calculated between the rapeseed (object) and other (background) sample data.

NDVI-Otsu Thresholding Method for Non-Vegetation Masking
SWIR is sensitive to water information, including subtle differences in soil moisture and leaf water [27]. To avoid the effects of rivers, lakes, ponds, and weather conditions (such as rain), NDVI was calculated to erase the non-vegetation area, including built-up, river, aqua-farm ponds, panels, and greenhouses, where NDVI was usually lower than vegetation [5,22,47]. The masking threshold was obtained using the Otsu algorithm (see Section 3.3) [23]. The formula of NDVI is described in Equation (1): where ρ nir and ρ red represent the reflectance of the NIR band and red band, respectively.

Adaptive Thresholding Method to Extract Rapeseed
As mentioned in Section 1, the commonly used empirical threshold method suffers from low applicability. Hence, this study adopted an adaptive thresholding method to overcome the weakness of uniform thresholds. Here, this study introduced two kinds of adaptive thresholding methods and developed a threshold-based rapeseed mapping method. To evaluate the threshold selected using the adaptive thresholding method, the empirical thresholding value was utilized as a comparison. The optimal empirical threshold value was chosen according to the overall accuracy (OA) calculated at every threshold value with steps of 0.01 from −0.30 to −0.10.
The Otsu algorithm is widely used in simple image segmentation where objects and backgrounds have huge differences [48]. It is a threshold selection method solely based on gray-level histograms with no prior knowledge [23]. The interclass variance is iteratively obtained from the histogram of the image for each threshold (see Equation (2)). The maximum interclass variance corresponds to the optimal threshold [22].
where BSS means the between sum-of-squares, describing the interclass variance quantitatively; p represents the number of categories, which is defined as 2 in this study; V is the value of the band used for segmentation; and category k is assigned by every V less than a specific threshold.
Although the Otsu algorithm provides a changeable threshold value for different environments, it has a limitation in that its accuracy becomes lower when facing a complicated image that contains several categories to classify [23] or the background is much larger than the object pixels [22]. To solve this problem, a combination of the Otsu thresholding method and Canny Edge Detection (OCED) was introduced [20,22]. There was a drastic difference at the boundary between rapeseed and non-rapeseed in NDRI, which could be detected by Canny Edge Detection (CED). It is noteworthy that, after the OCED process, the redundant pixels were cut, leaving the informative boundary pixels preserved. Thus, the histogram distribution represented the boundary of the rapeseed and other vegetation, containing a more balanced proportion of object and background. Thus, the accuracy of the Otsu algorithm improved [20,22].
In conclusion, the THAT automatic rapeseed mapping model can be described as follows: where P rapeseed means the pixels considered as rapeseed; P i means pixel i in the image; and NDV I i and NDRI i mean the NDVI and NDRI values of pixel P i , respectively. T 1 is the threshold derived from the Otsu thresholding algorithm with the whole NDVI image, while T 2 is the threshold derived from the adaptive algorithm (Otsu method or OCED model) with the NDRI image masked NDVI > T 1 .

Validation to Evaluate the Ability of the THAT Method
The accuracy assessment included two kinds of validation methods. First, samples were used to calculate the confused matrix, OA, user accuracy (UA), producer accuracy (PA), and kappa coefficient [49]. Second, the rapeseed mapping result based on OCED-THAT was compared with that based on the Otsu method, as well as widely used Random Forest classifiers; Random Forest is a stable, efficient, easily parameterized ensemble learning method and fits well in remote sensing classification [50]. Notably, only NDVI and NDRI were used in RF in order to maintain consistency with the THAT method.

Development of New Feature for Rapeseed Extraction
According to the spectral temporal analysis shown in Figure 4, rapeseed achieved the maximum value among all the vegetation in the visible red and green bands during the flowering period, as Davoud (2019) found in their study [9]. In contrast, the SWIR bands achieved the minimum value during the same period. Therefore, rapeseed achieved the maximum and minimum values during the flowering period, which follows the premise mentioned in Section 3.1. The green band was the most striking out of the visible bands, and SWIR1 was the most striking out of the SWIR bands. Accordingly, the green band and the SWIR1 band were selected as the optimal bands to extract rapeseed. The normalized difference between these two bands was computed so that the result could be confined to the range [−1, 1] [28,51]. In this way, the Normalized Difference Rapeseed Index (NDRI) [21] was developed. The formula of NDRI is described in Equation (4): where ρ green and ρ swir1 represent the reflectance of the green band and the SWIR1 band, respectively.   Figure 5a, a separation gap between rapeseed and other vegetation was found during the flowering period. The mean NDRI value of rapeseed was about −0.1, while the mean values of other vegetation were clustered around −0.35, theoretically forming a histogram with a bimodal pattern that can be segmented using the thresholding algorithm. Concerning the NDYI displayed in Figure 5c, the mean NDYI value of rapeseed was about 0.25 while that of wheat was around 0.2, and the other had a mean value of 0.12; therefore, a theoretical three-peaked histogram can be generated, which is challenge to segment with one single threshold. As for the NDVI and RRCI values illustrated in Figure 5b,d, thresholding segmentation was not an ideal method to extract rapeseed from these features.
As Table 5 illustrates, the JM distance indicates that NDRI outperformed other vegetation indices, including NDVI, NDYI, and RRCI, in rapeseed identification.

Threshold Selection for Rapeseed Extraction
For the non-vegetation mask, the NDVI threshold values were 0.4456 in Xinghua County and 0.3671 in Jiangsu Province.
For the extraction of rapeseed, the NDRI threshold values are displayed in Figure 6. In Xinghua County, the threshold values were −0.2267 and −0.1485 corresponding to the Otsu and OCED thresholding algorithms, respectively, as shown as Figure 6a,b. In Jiangsu Province, the threshold values were −0.2893 and −0.1648 corresponding to the Otsu and OCED thresholding methods, respectively, as shown as Figure 6c,d. As seen in Figure 7, the threshold value of −0.13 was the appropriate empirical threshold corresponding to the peak accuracy in Xinghua County. Conversely, this value was not appropriate in Jiangsu Province. The threshold calculated using the OCED algorithm outperformed the empirical threshold and the Otsu threshold.   Figure 8 depicts the overall extent of rapeseed distribution in both Xinghua and Jiangsu in 2020 using three methods. Figure 9 displays zoom-in view of sites A and B in Figure 8. In the optical true-color composite images of sites A and B displayed in Figure 9(1), rapeseed was yellow-green, and the other crops are dark green. By zooming in on three case regions of sites A and B (Figure 9a,c), it was found that the results of the OCED (demonstrated in Figure 9(2)) and Random Forest (shown in Figure 9(4)) had a high spatial consistency with the true-color optical Sentinel-2 image (Figure 9(1)), while the result of the Otsu suffered from noticeable errors as displayed in Figure 9(3). As Figure 6 shows, the NDRI threshold based on the Otsu was lower than that based on the OCED method. Conversely, as depicted in Figure 5a, rapeseed achieved a higher NDRI value than other vegetation. Therefore, a lower threshold derived from the Otsu method would lead to more commission errors. Table 6 illustrates that the Otsu-based rapeseed map in Xinghua County in 2020 ( Figure 8) had a high accuracy with the OA and kappa coefficient in the order of 0.9724 and 0.9327, respectively. However, when this method was used in Jiangsu Province, where the environment is larger and more complicated, the OA and kappa coefficient decreased to 0.9310 and 0.8089, respectively. The kappa coefficient was slightly improved in Xinghua County and remarkedly increased in Jiangsu Province based on the OCED thresholding algorithm. Compared with the RF-based method, the kappa coefficient of the THAT-based method was close to that in Xinghua County while slightly lower than that in Jiangsu Province. A comparison between the RF-based and the OCED-based rapeseed maps showed that the spatial distribution of these two products ( Figure 9) had a high consistency in Jiangsu Province. The OCED-based method was more likely to cause underestimation but at an acceptable level. This suggests the potential of the OCED-based method to extract rapeseed from other vegetation with optically sensed images.

Discussion
Several factors influence the accuracy of classification, namely, separability in spectral features between objects and backgrounds, the optimum threshold, and the availability of cloudless datasets [45].

The Potential of NDRI in Rapeseed Mapping
This study verified the potential of the NDRI vegetation index for rapeseed discrimination. The fact that rapeseed undergoes two obvious changes during the flowering stage provides theory support for the application of NDRI to identify rapeseed. Firstly, blooming makes the reflectance of the rapeseed field turn from foliar information to a mixture of foliar and floral signals [1]. Therefore, the intensity of the green band increases corresponding to the yellow-green blossom, since the yellow color is a mixture of green and red rays [1]. Secondly, during the reproductive stage, rapeseed needs more water, and evapotranspiration also increases [52,53]. As a result, the spectral reflection of the rapeseed field contains changes in water content. The SWIR1 band of rapeseed shows a significant reduction during the flowering stage, correlating to the increase in evapotranspiration, whereas other vegetation remains stable. Not only does NDRI show a huge separation between rapeseed and other vegetation, but it also keeps the value of rapeseed at its highest among all the vegetation. Thus, rapeseed and other vegetation form a bimodal image, which provides a means to apply the adaptive thresholding algorithm in crop mapping. It is noteworthy that Xu (2006) developed the Modified Normalized Difference Water Index (MNDWI) for water extraction, as it shared exactly the same band composite with NDRI [21].
However, the red-edge wavelength was not considered in this study because Landsat imagery does not contain red-edge waveband information, making harmonization with Sentinel-2 imagery challenging. Future studies may focus on the difference in red-edge information between rapeseed and other vegetation [8]. In addition to spectral features, the application of point clouds and plant height features [54] could also be considered in the future.

Adaptive Threshold Method in Crop Mapping
In this paper, the adaptive thresholding method was used to select an optimal NDVI threshold for non-vegetation masking and an optimal NDRI threshold for rapeseed extraction. The adaptive threshold method is widely used in algae monitoring since the spectral difference between water bodies and algae is very large. In contrast, there is a significant lack of adaptive threshold applications in crop mapping because the histograms of vegetation with the most commonly used indices are multimodal, which decreases the accuracy of the adaptive thresholding method [20,22,23]. Consequently, the threshold is usually empirically selected in previous studies [17,55], which affects the accuracy of the crop map products [7,44]. The results of our study indicate that the OCED method had good performance, with an OA of 0.9559 and a kappa of 0.8569.
However, we can envision more challenges when our method is replicated in large regional studies. Since GEE often needs zoning processing due to the user's memory limitations, some pseudo-rapeseed pixels will still be extracted when using the OCED method in the case of large-scale zoning processing if the rapeseed is not well distributed and there happens to be no rapeseed in some sub-regions. Thus, improvements should be made to detect whether there is any rapeseed. Only if there are any rapeseed pixels existing can the automatic rapeseed mapping algorithm operate.

Solution for Cloudy Regions
In this research, multi-sensor temporal aggregation was applied to improve temporal resolution in order to reduce cloud contamination. In the southern region of China, it is often cloudy and rainy, and optical remote sensing images are susceptible to cloud contamination, causing a reduction in the number of available images [13,22,56]. Additionally, temporal resolution is more critical for crop extraction because an insufficient temporal resolution can easily lose the information of crucial crop growth periods [55]. In particular, NDRI contains the SWIR band, which is sensitive to clouds. Accordingly, the influence of clouds should be considered [57]. High spatial and temporal resolution satellite datasets may be a good solution to reduce the influence of clouds [7]. Thus, in this study, Landsat 7, Landsat 8, and Sentinel 2 images were compiled into one dataset to improve the number of good-quality observations during the flowering period to obtain higher accuracy rapeseed mapping in a cloudy region.
As a cash crop, rapeseed is cultivated in patches in southern China, leading to a spectral mixture in remote sensing images [17]. Furthermore, rapeseed mapping requires a high spatial resolution image, which is in contrast to the mapping of bulk crops (e.g., wheat and paddy rice). However, the HLS image had a coarser spatial resolution (30 m) than the standalone Sentinel-2 image (10 m). Recently, the application of deep learning was considered in a study to extract rapeseed in cloudy regions, and it attained a satisfactory result [5]. Future research may focus on combining the deep learning technique and the threshold method [5]. In addition to optical remote sensing, combinations with radar data are expected to improve the accuracy of crop extraction in cloudy and rainy areas [8,13,26,43,56,58].

Conclusions
This study developed a new index for rapeseed identification called the Normalized Difference Rapeseed Index (NDRI) by analyzing the temporal profile of wavebands. The Jeffries-Matusita (JM) distance was utilized to evaluate the separability of rapeseed and other vegetation, and it was found that NDRI outperformed NDVI, NDYI, and RRCI in rapeseed mapping. Moreover, adaptive thresholding methods, including Otsu and Otsu with Canny Edge Detection (OCED), were conducted to choose an optimum threshold value. When comparing these methods, the OCED proved to be the best among the Otsu and empirical threshold methods. Based on NDRI and the adaptive thresholding method, a two-stepwise hierarchical adaptive thresholding method for rapeseed discrimination was developed in Xinghua County and extended to Jiangsu Province, China. The proposed method had a satisfying accuracy without training data, keeping a high consistency with pseudo-true-color composite images. In future research, more classification features should be fully utilized, such as red-edge features.

Conflicts of Interest:
The authors declare no conflict of interest.