Developing a phenology- and pixel-based algorithm for mapping rapeseed at 10m spatial resolution using multi-source data

As a major oilseed crop, large-scale and high-resolution maps of rapeseed (Brassica napus L.) are critical for predicting annual production and ensuring global energy security. However, such free maps are still unavailable in large areas. We designed a new pixeland phenology-based algorithm and produced a new data product for rapeseed planting area (20172019) over 33 countries at 10-m spatial resolution based on the multiple data. The product showed a good consistence (R=0.88) 10 with the official statistics (Food and Agricultural Organization of the United Nations, FAO) at national level. Rapeseed maps achieved at least 0.81 F1-scores of spatial consistency when comparing with the Cropland Data Layer (CDL) of America, Annual Crop Inventory (ACI) in Canada and Crop Map of England (CROME) in England. Moreover, their F1-scores ranged 0.84-0.92 based on the independent validation samples, implying a good consistency with ground truth. Furthermore, we found that rapeseed crop rotation is ≥2 years in almost all countries. Our derived maps with high accuracy suggest the robustness of 15 pixeland phenology-based algorithm in identifying rapeseed over large regions with various climate and landscapes. The derived rapeseed planting areas freely downloaded can be applied to predict rapeseed production and optimize planting structure. The product is available publicly at http://dx.doi.org/10.17632/ydf3m7pd4j.3 (Han et al., 2021).

Considering the unique phenological characteristics of crops, many studies have indicated the potential ways based on phenology for crop identifying in large areas (Ashourloo et al., 2019;Dong et al., 2016;Liu et al., 2018bLiu et al., , 2020aZhang et al., 2020). The algorithm based on phenology develops classification rules through analyzing the unique characteristic of the crop, which have been successfully applied to map rice (Dong et al., 2016), soybean (Zhong et al., 2014), corn (Zhong et al., 2016), 65 and sugarcane (Wang et al., 2020a), but not yet for rapeseed. Rapeseed has unique reflectance and scattering characteristics (Ashourloo et al., 2019;McNairn et al., 2018;Long, 2015, 2016), and undergoes three canopy morphologies (Ashourloo et al., 2019;Rondanini et al., 2014), including leaves, yellow petals, and pods/branches. Each canopy shape strongly influences how solar radiation intercepts (Sulik and Long, 2016). Compared with other crops, rapeseed is more easily to be identified, of which the yellow flowers significantly increase the reflectance of red and green bands (Pan et al., 2013;70 Sulik and Long, 2015). Additionally, when rapeseed grows, the backscatter signal increases because the rough canopy structure formed by the intertwined pods. (McNairn et al., 2018;Mercier et al., 2020;Tian et al., 2019;Veloso et al., 2017). Thus we are sure the specific featrues of reflectance values and scattering coefficients of rapeseed on S1/2 data will provide enough information for automatic rapeseed mapping over larger areas and with a finer resolution.
Thus, we integrated multi-source data to: 1) develop new method for identifying rapeseed ; 2) apply the new method to generate 75 rapeseed maps with a spatial resolution of 10m in three years (2017-2019) across the main planting areas of 33 countries; 3) analyze the geographic characteristics of rapeseed planting and crop rotation. The proposed algorithm and its derived products may benefit decision-makers, scientists and local farmers to ensure food and ennerge security.

Study area 80
This work identified rapeseed planting areas for 33 countries in three continents (North America, South America, and Europe) as they are the main rapeseed producers in the world (Fig. 1). The planting areas and production of rapeseed are large in Canada and the European Union (EU) (Carré and Pouzet, 2014;van Duren et al., 2015;Rondanini et al., 2012). According to the report S1c). Rapeseed planting seasons are different because of differences in natural conditions (such as climate and temperature) in different countries (Singha et al., 2019;Wang et al., 2018), which brings great challenges to map rapeseed. 95 Figure 1. The locations of 33 countries and the sample blocks for phenological monitoring with a radius of 10km.

Remote sensing data
We collected the Sentinel-2 (S2) and Sentinel-1 (S1) imagery (Table 1). The S1/2 satellites are launched by the European Space Agency (ESA) (Drusch et al., 2012;Torres et al., 2012). The highest spatial resolution of S2 images is 10m. We used 100 red (b4), green (b3), and blue (b2) spectral bands with 10 m spatial resolution Top-Of-Atmosphere (TOA) reflectance observations. The S2 TOA product includes the Quality Assessment (QA) band, which was used to remove most of the badquality images (e.g. clouds information) in this study. However, it is difficult to remove all clouds due to the quality of the QA band (Wang et al., 2020a;Zhu et al., 2015). The S1 includes four modes: Stripmap (SM), Interferometric Wide Swath (IW), Extra Wide Swath (EW), and Wave (WV) (Torres et al., 2012). We used the IW mode, which provides dual-band cross-105 https://doi.org/10.5194/essd-2021-34 polarization (VV) and vertical transmit/horizontal receive (VH) with a 12 day or 6-day repeat cycle and 10m space resolution (Torres et al., 2012). The S 1/2 images were obtained on GEE. Also, we used QA bands to remove most of the bad-quality images on GEE (Sample code can be found at https://code.earthengine.google.com/?scriptPath=Examples%3ADatasets%2FCOPERNICUS_S2). See Table 1 for more details 110

Digital elevation model
We used a spatial resolution of one arc-second (approximately 30m) elevation data from the Space Shuttle Radar Terrain Mission (Table 1) (Farr et al., 2007). Then we calculated the spatial distribution of slop on GEE (Sample code can be found at https://code.earthengine.google.com/?scriptPath=Examples%3ADatasets%2FUSGS_SRTMGL1_003) ( Fig. S1d-f). Later, we extracted areas with a slope less than 10° to mask hilly terrain where rapeseed is unlikely to be planted (Jarasiunas, 2016). 115

Cropland and agricultural statistics data
In this study, cropland data from the GFSAD30 were used to identify major farming areas in different countries (Phalke et al., 2020;Xiong et al., 2017). The existing crop data products containing rapeseed information include: 1) the 30-m Annual Crop Inventory (ACI) in Canada (Fisette et al., 2013). 2) the 30-m Cropland Data Layer (CDL) was generated in America (Boryan et al., 2011). CDL and ACI layers are downloaded on GEE. 3) the Crop Map of England (CROME) was generated in America 120 in England. These three crop layer products are generated based on satellite images and a large number of training sample collections. In this study, rapeseed maps in ACI, CDL and CROME were used for accuracy verification at the pixel level. The FAO releases annual statistics on the area for major crops in different countries or regions. We selected the statistics from FAO for accuracy verification. Please see Table 1 for more details

Crop calendars 125
We used two crop phenological data sets to assist in extracting rapeseed phenological parameters: (1) crop calendars in different countries (https://ipad.fas.usda.gov/ogamaps/cropcalendar.aspx). The crop calendars come from the United States Department of Agriculture (USDA) which only record the planting and harvest time of rapeseed in some countries (Table S1).
(2) field records of the crop phenology in Germany. In-situ observations come crop phenological records shared by the Deutsche Wetterdienst (DWD) in Germany (Kaspar et al., 2015). The DWD provides field observations of crop phenological 130 periods following the Biologische Bundesanstalt, Bundessortenamt, and Chemical (BBCH) scale throughout Germany (Table   1). DWD records the start date and the end date of rapeseed flowering (d'Andrimont et al., 2020;Kaspar et al., 2015). Note that both crop calendars and DWD do not contain information on the peak flowering dates of rapeseed. We used all stations that fully recorded the start and end of the flowering period from 2017 to . Finally, 281, 269, and 253 stations available in 2017 (the spatial distribution of the DWD rapeseed stations can be found in Fig. S2

Optical and SAR characteristics during the growing period of rapeseed
We selected available rapeseed parcels and in-situ observations of DWD from different climate regions and years to analyze the optical (reflectance and vegetation index) and SAR (VV, VH) characteristics of rapeseed along time. For example, Fig. 2 shows the time series of one rapeseed parcel around the DWD station with an id equal to 13126 in 2018. The rapeseed parcel 140 shows unique visual characteristics during the flowering period (Fig. S3). The flower is becoming yellow when rapeseed is approaching peak flowering (Pan et al., 2013;Tao et al., 2020;Wang et al., 2018). Rapeseed is yellow-green on the true color images of S2 and Google Earth during the flowering period (Fig. S4). The reflectance of green band and red band separately increased from 0.09 and 0.06 (2018/4/17, before flowering) to 0.16 and 0.14 (2018/5/7, peak flowering), and decreased after flowering (Fig. 2a). The reflectance of the blue band is lower than red and green bands during flowering. The reflectance of 145 the red band increased again and higher than the green band during the rapeseed harvest period. This is similar to previous research results ( Ashourloo et al., 2019;Sulik and Long, 2015). The Normalized Difference Yellow Index (NDYI, Eq.1) can capture the increasing yellowness in the time-series (d' Andrimont et al., 2020;Sulik and Long, 2016). Also, the NDYI increased from -0.03 on April 17 to 0.21 on May 7 (Fig. 2b). NDYI reaches a peak during the flowering time of rapeseed. This unique spectral feature of rapeseed in the flowering period is caused by the yellow petals. 150 where green is the TOA reflectance of the green band (b3) of the S2 imagery, blue is the blue band (b2) reflectance. S1 backscattering changes with the growth of rapeseed. We used the VV and VH time series smoothed by the Savitzky-Golay

Sample blocks collected for phenological monitoring 170
As a pre-requisite to map rapeseed at a large-scale, the phenology of rapeseed in different countries needs to be identified and 10km sample block over rapeseed producing areas of different countries and randomly sampled 10 rapeseed parcels uniformly for each block (Fig. S5). The rapeseed plots were identified by phenological characteristics obtained from visual-interpretation 175 and reference data: (1) high-resolution images available in S2 and Google Earth. It should be noted that the Google Earth images during rapeseed flowering were used to assist the visual interpretation of rapeseed parcels. (2) Spectral reflectance (red band and green band) and spectral index (NDYI) and scattering coefficient profiles (VV and VH) from S1/2 time series. The rapeseed parcels without high-quality available time-series imagery because of data quality issues such as clouds were omitted.
Finally, 75 sample blocks, 84 sample blocks, 84 sample blocks in 2017, 2018, 2019 were uniformly and randomly collected. 180 The sample blocks are shown in Fig. 1. We extracted the growth phenology information of rapeseed by calculating the average of the pixels of all rapeseed parcels in each block in different regions.

Flowering and pod phase detection in different countries
The phenology of rapeseed is different among regions. To find out the flowering dates of rapeseed in different countries, we evaluated each phenological sample block from 2017 to 2019. First, we calculated the time-series average value of S1/2 of all 185 sample rapeseed parcels pixels during the rapeseed growth period for each block in conjunction with the crop calendar. VV and VH time series for each sample rapeseed parcels were smoothed using the SG filter. Second, these S1/2 peak flowering dates and pod dates were derived for all sample blocks based on the method in Section 2.3.1. We found the peak flowering dates of rapeseed have an explicit latitude gradient, especially in Europe (Fig. S6). Also, we found the signal with the maximum of VH is within 45 days after the peak date of flowering (Fig. S7). Then we calculated the difference in the peak flowering 190 date of each sample block in different years. The results showed that the error were within 10 days of advance or delay (Fig.   3d). Therefore, it is reasonable to use the same time for rapeseed identification in different years in the same area in this study.
Previous studies and field observation records show that the flowering period of rapeseed is about 30 days (d' Andrimont et al., 2020;Chen et al., 2019;Kaspar et al., 2015;She et al., 2015). Therefore, we divide each month into two-time periods (the 15th is the dividing line). Two consecutive half-months are defined as suitable classified flowering periods. Finally, we 195 designated the flowering period for each sample block based on the peak flowering dates (Fig. 3a-c).

Developing phenology-and pixel-based algorithm for mapping rapeseed
The temporal profile analysis at rapeseed fields in this study and many previous studies indicated that the spectrum at flowering stage and the scattering signal at pod stage are the key features to identify rapeseed (Ashourloo et al., 2019;Bargiel, 2017;Han et al., 2020;Mercier et al., 2020;Sulik and Long, 2015;Veloso et al., 2017). Previous studies have found the high 205 reflectance values of green band and red band at the flowering stage are the main spectral factors to distinguish other crops (Ashourloo et al., 2019). We developed phenology-and pixel-based rapeseed mapping algorithm for rapeseed mapping using four features: including 1) spectral bands (red and green), 2) spectral indices (NDYI), 3) polarization bands (VH), and 4) terrain (slope). And four steps were orderly conducted for mapping annual planting areas (Fig. 5). Step 1: determining the threshold of the feature indicators. Thresholds of indicators are the key parameters to determine the 210 area accurancy identified. We analyzed the histograms of the random samples selected from different countries as the previous studies (Zou et al., 2018) suggested. We found consistently the normal distributions of green, blue, and NDYI (Fig. S9) for all samples during the flowering stage from different regions. 98% of the rapeseed pixels show red > 0.07, green > 0.11, and NDYI > 0.05. However, we found some pixels with a higher NDYI, while in fact they are polluted by the cloud with a "rainbow" appearance, 215 would be misclassified into rapeseed. Such misclassifications caused by some bad-quality observations from the S2 image cann't be removed due to the limited quality of the QA band and simple cloud score algorithm (Wang et al., 2020c;Zhu et al., 2015). The "rainbow" of the cloud is come from the push-broom design of S2 (Fig. 4a) and spectral misregistration (For more details, please refer to ESA, 2015a, andESA, 2015b). Based on the principle of the relative displacement of each spectral channel sensor in the S2 push-broom design (Frantz et al., 2018;Liu et al., 2020b;Zhao et al., 2018), we developed a new 220 spectral index (NRGBI) to reduce the influence of "rainbow" (Eq.2). The scatter plot of NDYI and NRGBI of rapeseed field samples and "rainbow" around cloud samples (visual interpretation) show that the NRGBI (threshold is -0.05) can effectively distinguish rapeseed from the "rainbow" (Fig. 4h).  by quality assurance band (QA60)) for Sentinel-2 TOA image, with the red arrows indicating the cloud "rainbow" appearance at high altitude in the S2 image (image source: Copernicus Sentinel-2 data); (g) Sentinel-2 TOA image of rapeseed at the flowering stage, with the yellow arrow for the rapeseed fields; (h) scatter plots of NDYI and NRGBI of rapeseed field samples and "rainbow" around clouds samples from the S2 images, with the color density for the number of pixels. 235 Step 2: identifying all rapeseed pixels from different images during the flowering period and aggregating them into annual rapessed planting areas (Fig. 5). Because the peak flowering dates and the number of available images of rapeseed fields in a region are different (Fig. S10), rapeseed classifications based on single image could fail in capturing rapeseed flowering dynamics (Ashourloo et al., 2019). To avoid the misclassification from bad-quality observations during rapeseed flowering stage, the aggregate approach from all the available S2 images during the flowering period was used. Hence, a larger number 240 of images will result in better performance (Fig. S10). Step 3: combining optical with SAR images to ensure the accuracy of the rapeseed maps. The high VH values during the rapeseed pod stage is another distinct feature that distinguishes it from other crops (Mercier et al., 2020;Tian et al., 2019;Van Tricht et al., 2018;Veloso et al., 2017). Considering the variability of flowering in different fields and the duration of pod stage (Section 2.3.2), we calculated the maximum VH during the period between the second half of the flowering stage and 245 30 days after it (~45 days). Within 45 days, at least three S1 satellite images are available in the study areas. Also, the areas with a slope >10° (where rapeseed is unlikely to be planted) were removed (Jarasiunas, 2016). All pixels that meet the requirements are defined as rapeseed.
Step 4: removing the "salt and pepper" noise according to the connected domain threshold of the pixels and filling the gaps inside the parcels (Hirayama et al., 2019). The thresholds of different indicators in different regions can be found in Table S3. 250

Accuracy assessment
Total three ways were appled to assess the accuracy of rapeseed map acquired by our new phenology-and pixel-based algorithm. First, we compared the rapeseed areas retrieved by the new method with FAO statistics. Our rapeseed data is a 260 binary (0 or 1) map with a spatial resolution of 10m, equivalently as each pixel = 100 m 2 . We then can calculate the total area of rapeseed map in each country and compare them with the national rapeseed statistics. The RMSE (Eq.3), and the MAE (Eq.4), and the coefficient of determination (R 2 , Eq.5) are used to verify the accuracy of rapeseed mapping.

= √∑ ( − ) 2 =1
(3) where n is the total number of countries. is the mapped rapeseed planting areas, ̅ is the corresponding mean value, is the records rapeseed planting areas from FAO, ̅ is the corresponding mean value.
Also, we compared the rapeseed maps with two open access datasets that include rapeseed layers (ACI, CDL and CROME) in Canada, America and England on the pixel scale. We used them as the reference data for 2018, 2019 (Boryan et al., 2011;270 Fisette et al., 2013). To unify the spatial resolution of the rapeseed maps, CDL, ACI and CROME were resampled to 10m resolution for comparison. UA (Eq.6), PA (Eq.7), and F1 (Eq.8) were calculated based on confusion matrices (Table S2) to measure the classification accuracy.
Thirdly, we also selected randomly verification samples based on the previous studies (Pekel et al., 2016;Wang et al., 2020b) to validate the rapeseed map. A grid (0.2 latitude by 0.2 longitude) was generated within the rapeseed map in 2018 acquired 275 by our pixel-and phenology-based method. Two points (rapeseed and non-rapeseed) was generated randomly in each grid by interpreting visually images available in the S2 and Google Earth, together with spectral reflectance, spectral index and scattering coefficient profiles from S-1/2 time series. The confusion matrices were similarly used to assess the accurancy according to Eqs 6~8.

=
(6) 280 Where xij is the value of the i-th row and j-th column; xi is the sum of the i-th row; xj is the sum of the j-th column. Although the statistical data and existing products are not completely same as the real areas and locations of rapeseed planted on the ground, these datasets do benefit to validate the accurancy of rapeseed maps at different scales (national and pixels).

Accuracy assessment
We compared the derived rapeseed areas with those from the FAO statistics. The total planting areas of rapeseed are well consistent with the agricultural statistics at the national level, with RMSE of 1459.64 km 2 , MAE of 785.25 km 2 , amd R 2 of 0.88 (Fig. 6). We found the derived areas in 2017 and 2019 are larger than those in 2018, especially for the countries with the 290 relative smaller rapeseed areas (e.g. many European countries indicated by the small plot located in the bottom right and below of Fig. 6. The more availability of S2 images together with better quality in 2018 could contribute the larger rapeseed areas derived by the new method (Liu et al., 2020a). The comparsion of our rapeseed maps with those of America CDL in 2018, 2019, Canada ACI and England CROME in 2018 was consistently indicated by a higher accurancy according to the confusion matrices ( Table 2)  The confusion matrices (Table S4) based on random sampling points show that the accuracy of the rapeseed maps varies in different regions. We found Europe shows the highest accuracy (F1 0.91), followed by Chile (F1 0.9), and the lower F1 (0.84) indicated by North America. Such disparity in accuracy might be ascribed to the different availability of high-quality images in the studied areas. Nevertheless, all rapeseed maps derived by our method show a reasonable accuracy across such a larger 305 region.

More details of rapeseed maps derived by the new method
To show more details of rapeseed maps derived by our method, we captured some images in some areas of each country. The rapeseed maps show good spatial consistency with the actual rapeseed planted on the ground (Fig. 7). From the area densely planted by rapeseed in Canada ( Fig. 9-a) to relatively sparse planting ones such as in Chile ( Fig. 9-b) and European countries 310 (e.g. Fig. 9-c,d) (Lowder et al., 2016), from regular rectangles (e.g. Fig. 9-a, h) to irregular parcels ( Fig. 9-c, d), from temperate oceanic climate (Fig. 9-c-e) to temperate sub-continental ( Fig. 9-a, f), or even subtropic climate ( Fig. 9-b), all filed details were indicated clearly in our maps. Fragmentation of land in some European countries, especially in Eastern and Central Europe after land reform in 1989 (Hartvigsen, 2013(Hartvigsen, , 2014, such as Estonia (Fig.7f) (Jürgenson and Rasva, 2020;Looga et al., 2018). Although under different climate, terrain, landscapes and over very larger region, the algorithm proposed in our study 315 consistently shows a good classification accuracy across 33 countries. Thus, the rapeseed maps based on S1/2 data can effectively identify the fields in detail with high spatial resolutions and clear field boundaries. More rapeseed classification details can be found in Fig. S12 and Fig. S18. Furthermore, our results showed a consistent distributions between our rapeseed maps and the existing products at the pixel level (Fig. S13). The yellow grids (71%~80%) mean they are simultaneously identified as rapeseed areas both by our method 320 and ACI/CDL/CROME, while red grids for disagreement. The difference in accuracy might be caused by the number of highquality images available in different regions (Dong et al., 2016). Despite the various ground conditions, methods, images, and spatial resolutions among the products, the comparison results further verify the accuracy of our rapeseed map (Gong et al., 2020;Singha et al., 2019).

Spatial patterns of rapeseed planting areas 330
Canada shows the largest rapeseed planting area (Fig. 8, Fig. S14), with the total area of 118,489.73 km 2 in 2018, higher than those in Europe (106,814.67 km 2 ). France and Germany are two leading rapeseed growing countries in Europe, accounting for around 66.3% European rapeseed areas together with other three countries (United Kingdom, Poland, and Ukraine). The country-wide rapeseed areas in all 33 countries were further normalized to show clearly the spatial patterns (Fig. S14). The spatial patterns of three years (2017~2019) are consistent at the national level. Moreover, we also plotted the geographic 335 characteristics of rapeseed areas along latitude and longitude by three regions (Fig. 8). Rapeseed in Europe is widely planted in the countries with a latitude of 45~56°N and longitude of -2°4°, 9°19°, and 22°27°, with exception of the steep mountainous areas and the cold northern areas (Fig. 8a) (van Duren et al., 2015). In North America, the areas with latitude of 44~44.5°N, 51~55°N, 56~57°N and longitude of -118°-117°, -114°-98°are densely distributed by rapeseed (Fig. 8b).

Investigating the rapeseed rotation systems
We obtained three-year rapeseed maps at a 10m spatial resolution, and with a higher accurancy which was validated by annual 345 national statitisc books, open accessed public products and random sampling points at 0.2°x 0.2°grids. These rapeseed maps, with good quality for three consecutive years, provide a new opportunity to investigate rapeseed rotation systems (Liu et al., 2018a). Crop rotation information is considered as an important factor for yield, soil degradation, and greenhouse gas emissions (Harker et al., 2015;Liu et al., 2018a;Ren et al., 2015;Rudiyanto et al., 2019;Zhou et al., 2015). Thus, we selected 25 representative areas (Fig. S15) to analyze the rapeseed rotation patterns according to the three criterions: 1) high-quality images 350 available during the annual rapeseed flowering period from 2017 to 2019; 2) high rapeseed classification accuracy; 3) rapeseed is the main crop type and a large area is planted by rapeseed every year. The rapeseed rotation was calculated by the frequency in each rapeseed pixel ( Fig. S16-17).
We classfied the rapeseed rotation break into three types: ≥2 years, 1 year, 0 year ( Fig. 9 and Fig. S16-17. We found most countries show the rotation break ≥2 years (the highest ratios of green parts) (Fig. 9), especially for European countries (Fig.  355 9-b). The rotation break ≥2 years in Canada mounts to 70%, followed by 30% with 1-year break ( Fig. 9-a). The histogram futher confirmed that 20 locations have identified as the rotation break ≥2 years, with their rates of planting areas higher than 90% ( Fig. 9-d). Many previous studies have found that a two-or three-year rotation break will significantly reduce the number of spores, especially rhizomes and blacklegs, suggesting rotation system is an important step in controlling diseases (Gill, 2018;Harker et al., 2015;Ren et al., 2015;Zhou et al., 2015). Moreover, rapeseed rotation will also benefit yield, insects, moisture, 360 fertility, and reducing weeds (Bernard et al., 2012;Harker et al., 2015;Pardo et al., 2015;Peng et al., 2015;Ren et al., 2015).
Thus, more efforts should be input to produce longer time-series rapeseed maps, and obtain detailed rotation information in the future.

Uncertainty analysis
Generating annual high-resolution maps of a specific crop over a larger region is a big challenge (Gong et al., 2020;Liu et al., 2018bLiu et al., , 2020a. Pixel-and phenological-based algorithms, multisource remote sensing data, and the GEE are useful to map rapeseed with high-resolution and over larger areas. Besides, the algorithm does not need training sample data and reduces 370 disturbance from agronomy differences by combining images of multiple dates. However, the uncertainty are from the follows: 1) cropland layer. We used the GFSAD30 datasets to identify croplands. However, GFSAD30 has its limitations such as classification error (Phalke et al., 2020). 2) the number of satellite images available. Although our annual rapeseed maps are consistent with FAO statistics, and show higher accuracy comparing with existing products, the maps are limited by the goodquality observations during the critical growth stages. For example, Fig. 10a shows that the error in the area of France in 2017, 375 and error could be attributed to the lack of clear S2 images during rapeseed flowering period (Fig. 10b). Rapeseed flowering period is generally characterized by higher NDYI, red band, and green band reflectance, thus missing the peak will fail in https://doi.org/10.5194/essd-2021-34 identifying rapeseed pixels (Fig. 10c). 3) thresholds for different indicators. The threshold is the key for mapping crop (Ashourloo et al., 2019;Dong et al., 2016;Liu et al., 2020a;Wang et al., 2020a;Zhang et al., 2015). Although the reference thresholds for three regions are given in this study, it should be cautious when applying them into other regions. 4) the 380 complexity of ground environment. For example, landscape types might impact the accuracy of rapeseed map (Wang et al., 2020a).

Implications and improvements
Despite the above limitations, the new phenology-based method proposed by us has the potential to extend to other regions by revising the phenology metrics. Recently, the Harmonized Landsat and Sentinel-2 database improve the spatial resolution and revisit cycle of images (Claverie et al., 2018;Shang and Zhu, 2019). Similar or even higher rapeseed classification accuracy 390 can be expected. Furthermore, remote sensing data fusion algorithms have been continuously developed (e.g., STARFM and ESTARFM) (Zhu et al., 2010). Finnally, various deep learning models have been explored for classifying crops and lowering errors (Hu et al., 2019;Zhong et al., 2019). Integrating phenological metrics and deep learning models might further improve rapeseed mapping accuracy. Thus, such rapeseed products will objectively track the dynamics of rapeseed planting areas as well as agricultural management in the future. 395

Data availability
The rapeseed map produced is accessible at Mendeley Data (http://dx.doi.org/10.17632/ydf3m7pd4j.3) (Han et al., 2021). The rapeseed maps with 10 m resolution are provided in this study. The dataset includes a set of GeoTIFF images in the ESPG: 4326 spatial reference system. The values 1 and 0 represent rapeseed and non-rapeseed, respectively. We encourage users to https://doi.org/10.5194/essd-2021-34 independently verify the rapeseed map. In addition, S 1/2 images, CDL, ACI, and SRTM are available on GEE 400 (https://developers.google.com/earth-engine/datasets/). For more detailed information about the data collected in this work, please see Table 1.

Conclusions
Large-scale and high-resolution rapeseed maps are the basis for crop growth monitoring and production prediction. We designed a new method for mapping rapeseed based on the spectral and polarization feature and multi-source data. The new 405 algorithm have produced three annual rapeseed maps (2017~2019) at 10m spatial resolution in 33 countries. Three different verification methods indicate that our products have relatively higher accuracy. Our approach reduces disturbances from different planting time and bad-quality observations in some degree. The 10m rapeseed maps do provide more spatial details of rapeseed. Also, we have found that rapeseed crop rotation is 2 years or more in almost all countries. The rapeseed mapping method proposed in this study could be applied into other regions and other crops. The derived rapeseed data product is useful 410 for many purposes including crop growth monitoring and production, rotation systems planning.

Author contributions
ZZ and JH designed the research. JH and LY collected datasets. JH implemented the research and wrote the paper. ZZ, JC, LZ, JZ, and ZL revised the paper. 415

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.