High-Resolution Mapping of Winter Cereals in Europe by Time Series Landsat and Sentinel Images for 2016–2020

: Winter cereals, including wheat, rye, barley, and triticale, are important food crops, and it is crucial to identify the distribution of winter cereals for monitoring crop growth and predicting yield. The production and plating area of winter cereals in Europe both contribute 12.57% to the total global cereal production and plating area in 2020. However, the distribution maps of winter cereals with high spatial resolution are scarce in Europe. Here, we ﬁrst used synthetic aperture radar (SAR) data from Sentinel-1 A/B, in the Interferometric Wide (IW) swath mode, to distinguish rapeseed and winter cereals; we then used a time-weighted dynamic time warping (TWDTW) method to discriminate winter cereals from other crops by comparing the similarity of seasonal changes in the Normalized Difference Vegetation Index (NDVI) from Landsat and Sentinel-2 images. We generated winter cereal maps for 2016–2020 that cover 32 European countries with 30 m spatial resolution. Validation using ﬁeld samples obtained from the Google Earth Engine (GEE) platform show that the producer’s and user’s accuracies are 91% ± 7.8% and 89% ± 10.3%, respectively, averaged over 32 countries in Europe. The winter cereal map agrees well with agricultural census data for planted winter cereal areas at municipal and country levels, with the averaged coefﬁcient of determination R 2 as 0.77 ± 0.15 for 2016–2019. In addition, our method can identify the distribution of winter cereals two months before harvest, with an overall accuracy of 88.4%, indicating that TWDTW is an effective method for timely crop growth monitoring and identiﬁcation at the continent level. The winter cereal maps in Europe are available via an open-data repository. methodology, X.H. and W.Y.; software, X.H.; validation, J.W.; formal analysis, J.D.; resources, W.Y.; writing—original draft preparation, X.H. and W.Y.; writing—review and editing, Y.F., B.P., Y.Z., S.S. and W.Y.; visualization, X.H. All authors have read the published version of the manuscript.


Introduction
Winter cereals, such as wheat, rye, barley, and triticale, are important cereal crops, and approximately 278 million hectares are planted in 2020, accounting for 30% of the world's grain-crop area and 41% of global grain production (https://www.fao.org/faostat/en/ #data, accessed on 10 May 2020). Global grain prices fluctuate during the wheat-growing season, as supply expectations shift and change significantly in response to the amount of planted acreage, the weather, and growing conditions [1]. Winter cereals are also perhaps the most political commodity in the world because they are the main ingredients of bread, which is the most basic food. The Arab Spring of 2010 started as a direct result of bread resolution crop maps (e.g., winter cereals) at national and local sites, which were supported by the "Sentinel-2 for agricultural" project [30,33,34]. Under the framework of the agricultural policy of the European Union, the Land Parcel Identification System (LPIS) opened in France was developed to identify all declared agricultural parcels through aerial photography and a geographic information system (GIS) at a scale of 1:5000 [35]. However, there are still some hindrances in the existing winter cereal maps. Firstly, it is difficult to classify on a larger scale due to the fixity of sites [30]. Secondly, the generation of crop maps based on the LPIS platform is time-consuming [36], so it is difficult to provide real-time maps.
Crop maps produced before the harvest, especially in areas with water resources problems, are requested by policy and decision makers for use in agriculture monitoring, simulating crop water use, calculating statistics, and satisfying the timeliness of yield prediction demands [10,37]. The frequency and intensity of drought are projected to increase in Europe in the future as the climate changes, leading to a risk of crop yield losses and even crop failure [36,38,39]. The early-season (as early as possible) mapping of winter cereals is therefore urgently needed for yield forecasting and assessments of food security [5,40]. However, due to the limited input information, it is more difficult to identify the distribution of crops before harvest than at the end of the growing season [5].
In this study, we used a phenology-based method (TWDTW) to detect the geographical locations of winter cereals in Europe and mapped winter cereal distribution for 2016-2020 with the resolution of 30 m. Furthermore, we assessed the performance of the TWDTW method in identifying winter cereals at the early stage of the growing season in Europe. The winter cereal map can be updated annually and used as evidence for agricultural management and policy making.

Study Area
Our study area covers 32 European countries ( Figure 1). The area has highly diverse climatic and topographic characteristics. The temperate maritime climate and the Mediterranean climate are the main climates in Europe, and the agricultural landscape is characterized by medium sized farms (2-25 ha) [41].   We used surface reflectance (SR) products, after atmospheric corrections, from the Enhanced Thematic Mapper Plus (ETM+) sensor onboard Landsat 7, the Operational Land Imager (OLI) sensor onboard Landsat 8 (United States Geological Survey), and the Multispectral Instrument (MSI) sensor onboard Sentinel-2 (European Space Agency), to obtain the NDVI from 2016 to 2020. Firstly, to reduce the impact of clouds, the quality bands sourced in Landsat and the QA60 band provided by Sentinel-2 were employed to remove the clouds in their respective images. Then, we used the nearest neighbor method to resample the NDVI of Sentinel-2 to 30 m, so as to keep the same spatial resolution with Landsat. Furthermore, we obtained all valid NDVI values (Landsat 7, Landsat 8, and Sentinel-2) within one month with cloud-free, and chose the maximum value for monthly synthesis for each pixel in each image, which was proved to be effective for crop mapping and to show the growing state. In addition, the monthly cloud-free image frequencies from October to the following July, over 2016-2019, were calculated for each pixel and shown in Figure 2. The monthly maximum composite NDVI was implemented on GEE and then downloaded to the local server. For Sentinel-2, the Level-2A (SR) product cannot be obtained from either GEE or ESA prior to 2018 due to missing data. We used surface reflectance (SR) products, after atmospheric corrections, from the Enhanced Thematic Mapper Plus (ETM+) sensor onboard Landsat 7, the Operational Land Imager (OLI) sensor onboard Landsat 8 (United States Geological Survey), and the Multispectral Instrument (MSI) sensor onboard Sentinel-2 (European Space Agency), to obtain the NDVI from 2016 to 2020. Firstly, to reduce the impact of clouds, the quality bands sourced in Landsat and the QA60 band provided by Sentinel-2 were employed to remove the clouds in their respective images. Then, we used the nearest neighbor method to resample the NDVI of Sentinel-2 to 30 m, so as to keep the same spatial resolution with Landsat. Furthermore, we obtained all valid NDVI values (Landsat 7, Landsat 8, and Sentinel-2) within one month with cloud-free, and chose the maximum value for monthly synthesis for each pixel in each image, which was proved to be effective for crop mapping and to show the growing state. In addition, the monthly cloud-free image frequencies from October to the following July, over 2016-2019, were calculated for each pixel and shown in Figure 2. The monthly maximum composite NDVI was implemented on GEE and then downloaded to the local server. For Sentinel-2, the Level-2A (SR) product cannot be obtained from either GEE or ESA prior to 2018 due to missing data. Moreover, we used synthetic aperture radar (SAR) data from Sentinel-1 A/B, in the Interferometric Wide (IW) swath mode, to obtain the radar response to the physiological Moreover, we used synthetic aperture radar (SAR) data from Sentinel-1 A/B, in the Interferometric Wide (IW) swath mode, to obtain the radar response to the physiological Remote Sens. 2022, 14, 2120 5 of 16 structural development of crops [42]. The SAR data offered vertical transmittance/vertical receival (VV) and vertical transmittance/horizontal receival (VH) polarization bands with 10 m spatial resolution. Similarly, we resampled the SAR data to 30 m based on the nearest neighbor method and obtained the monthly maximum synthetized values of VH. The above processes were run on the GEE platform.

Field Samples and Agricultural Census Data
Winter cereals are biennial cereal crops sown in the autumn and harvested in the following July. Wheat, winter barley, rye, and triticale are typical winter cereals in Europe and they had similar seasonal change curves of NDVI and phenological characteristics [43]. In this study, we first obtained 50-200 winter cereal field samples for each country in 2018. These samples were selected randomly by visually interpreting the color and textures of images on Google Earth [5] and comparing the NDVI thresholds with empirical values [44]. Further, these samples were evenly distributed in space. Then, the standard NDVI seasonal change curves of winter cereals in each country were obtained by calculating the average NDVI value of field samples in the corresponding countries ( Figure 3). structural development of crops [42]. The SAR data offered vertical transmittance/vertical receival (VV) and vertical transmittance/horizontal receival (VH) polarization bands with 10 m spatial resolution. Similarly, we resampled the SAR data to 30 m based on the nearest neighbor method and obtained the monthly maximum synthetized values of VH. The above processes were run on the GEE platform.

Field Samples and Agricultural Census Data
Winter cereals are biennial cereal crops sown in the autumn and harvested in the following July. Wheat, winter barley, rye, and triticale are typical winter cereals in Europe and they had similar seasonal change curves of NDVI and phenological characteristics [43]. In this study, we first obtained 50-200 winter cereal field samples for each country in 2018. These samples were selected randomly by visually interpreting the color and textures of images on Google Earth [5] and comparing the NDVI thresholds with empirical values [44]. Further, these samples were evenly distributed in space. Then, the standard NDVI seasonal change curves of winter cereals in each country were obtained by calculating the average NDVI value of field samples in the corresponding countries ( Figure 3). In this study, we used the same standard curve for these latter countries because they are geographically adjacent.
In addition, we selected 3356 field samples in 32 countries (1502 samples for winter cereals and 1854 samples for other, non-winter cereal crops) in 2018 for validating the identification accuracy at the pixel level. The method of generating these validation samples is similar to the above procedure. In addition, the winter cereal plating area at country and municipal levels in the study area were obtained from the National Bureau of Statistics in each country to assess the accuracy at the regional level (https://www.fao.org/faostat/en/#data, accessed on 25 May 2020). Overall, we obtained the country-level agricultural census data for all 32 countries from 2016 to 2020 and municipal-level data for 11 In this study, we used the same standard curve for these latter countries because they are geographically adjacent.
In addition, we selected 3356 field samples in 32 countries (1502 samples for winter cereals and 1854 samples for other, non-winter cereal crops) in 2018 for validating the identification accuracy at the pixel level. The method of generating these validation samples is similar to the above procedure. In addition, the winter cereal plating area at country and municipal levels in the study area were obtained from the National Bureau of Statistics in each country to assess the accuracy at the regional level (https://www.fao. org/faostat/en/#data, accessed on 25 May 2020). Overall, we obtained the country-level agricultural census data for all 32 countries from 2016 to 2020 and municipal-level data for Remote Sens. 2022, 14, 2120 6 of 16 11 countries from 2016 to 2019. Therefore, our validation with respect to the municipal-level data was implemented for 11 countries.

Method
The methodology used in this study is shown in the flow chart ( Figure 4). The workflow consists of the following steps: (1) constructing monthly maximum composite NDVI images; (2) using SAR data to distinguish winter rapeseed at each pixel; (3) preparing the standard seasonal change in NDVI-based winter cereal samples; (4) winter cereal identification with the TWDTW method; (5) evaluation of classification accuracies.
Remote Sens. 2022, 14, x FOR PEER REVIEW 6 of 18 countries from 2016 to 2019. Therefore, our validation with respect to the municipal-level data was implemented for 11 countries.

Method
The methodology used in this study is shown in the flow chart ( Figure 4). The workflow consists of the following steps: (1) constructing monthly maximum composite NDVI images; (2) using SAR data to distinguish winter rapeseed at each pixel; (3) preparing the standard seasonal change in NDVI-based winter cereal samples; (4) winter cereal identification with the TWDTW method; (5) evaluation of classification accuracies.

Time-Weighted Dynamic Time Warping (TWDTW)
We used the Time-Weighted Dynamic Time Warping (TWDTW) method to identify winter cereals in Europe. The TWDTW method is an improved algorithm of Dynamic Time Warping (DTW), which was originally developed for speech recognition [19]. The details about TWDTW can be found in previous studies [5,24,26].
First, we determined the standard NDVI seasonal change curves [5,25] of winter cereals from the average NDVI values of field samples for each country in 2018. Then, we calculated the distance between the time series of each image pixel and the standard NDVI seasonal change curve of winter cereals by using the TWDTW method. All distance values were ranked from small to large and pixels with low distances were considered to have a higher probability of being winter cereals. In addition, we used the agricultural census area of winter cereals for each country to obtain distance thresholds, which means the total area of winter cereal pixels were consistent with the agricultural census area of each country. Moreover, we presumed that the standard NDVI seasonal change curves for each country did not change from year to year. We employed the curves in 2018 to identify the planting area of winter cereals from 2016 to 2017 and 2019 to 2020, and assessed the potential of the TWDTW method in temporal transferability.

Removing the Disturbance of Winter Rapeseed
Winter rapeseed is sown in early September and harvested the following July or August, which is similar to winter cereals. It is difficult to distinguish between winter cereals and winter rapeseed in optical imagery due to the similarity of the seasonal changes in NDVI. However, rapeseed is taller than other winter cereals and has randomly oriented branches, resulting in high volume scattering and lower attenuation of the signal from the ground [42]. Thus, the high and increasing VH backscatter of rapeseed at late growth stages can be used to discriminate between winter cereals and winter rapeseed [42]. In this study, we referred to Veloso et al. and Dong et al. [5,42],and considered that the VH backscatter values for winter rapeseed in May were above −15.5, and that of other winter cereals were less than −15.5. Accordingly, we assigned higher distance values to the pixels

Time-Weighted Dynamic Time Warping (TWDTW)
We used the Time-Weighted Dynamic Time Warping (TWDTW) method to identify winter cereals in Europe. The TWDTW method is an improved algorithm of Dynamic Time Warping (DTW), which was originally developed for speech recognition [19]. The details about TWDTW can be found in previous studies [5,24,26].
First, we determined the standard NDVI seasonal change curves [5,25] of winter cereals from the average NDVI values of field samples for each country in 2018. Then, we calculated the distance between the time series of each image pixel and the standard NDVI seasonal change curve of winter cereals by using the TWDTW method. All distance values were ranked from small to large and pixels with low distances were considered to have a higher probability of being winter cereals. In addition, we used the agricultural census area of winter cereals for each country to obtain distance thresholds, which means the total area of winter cereal pixels were consistent with the agricultural census area of each country. Moreover, we presumed that the standard NDVI seasonal change curves for each country did not change from year to year. We employed the curves in 2018 to identify the planting area of winter cereals from 2016 to 2017 and 2019 to 2020, and assessed the potential of the TWDTW method in temporal transferability.

Removing the Disturbance of Winter Rapeseed
Winter rapeseed is sown in early September and harvested the following July or August, which is similar to winter cereals. It is difficult to distinguish between winter cereals and winter rapeseed in optical imagery due to the similarity of the seasonal changes in NDVI. However, rapeseed is taller than other winter cereals and has randomly oriented branches, resulting in high volume scattering and lower attenuation of the signal from the ground [42]. Thus, the high and increasing VH backscatter of rapeseed at late growth stages can be used to discriminate between winter cereals and winter rapeseed [42]. In this study, we referred to Veloso et al. and Dong et al. [5,42], and considered that the VH backscatter values for winter rapeseed in May were above −15.5, and that of other winter cereals were less than −15.5. Accordingly, we assigned higher distance values to the pixels with VH greater than −15.5, which means that those pixels are less likely to be misclassified as winter cereal.

Accuracy Assessment
We evaluated the accuracy of identification of winter cereals at regional and pixel levels by using the census dataset and the field samples. At the regional (municipal) level, the coefficient of determination (R 2 ), relative error (RE), and root mean squared error (RMSE) were used to compare the identified planted area with the agricultural census data for winter cereals. At the pixel level (3393 field samples), the producer's accuracy (PA), user's accuracy (UA), overall accuracy (OA), and Kappa coefficient (KC) [45] were used to assess the accuracy of the identified maps. The PA is the percentage of the ground-sampled winter cereal fields that were correctly identified as winter cereals. The UA is the percentage of the identified winter cereal fields that correspond to winter cereal fields in the field samples. The OA shows the overall effectiveness of the identification. In addition, for the Kappa coefficient, higher values indicate better classification accuracy.

Results
We examined how early winter cereals can be identified with an overall accuracy reaching 88%. The TWDTW method was used to identify the winter cereals based on the NDVI time series with different lengths starting at November, then increased by monthly increments thereafter. All 3356 validation field samples were used to assess the accuracies of the early-season mapping of winter cereals. With the increase in the length of the time series, the classification accuracy (i.e., OA, PA, and UA) also increased ( Table 1). The overall accuracy reached 88.4% in May, close to the maximum overall accuracy (91.8%) at the end of season, with a stable accuracy. These results imply that winter cereals can be identified two months before harvest (i.e., in May) by using the TWDTW method. We produced early-season winter cereal maps for 2016-2020 with 30 m resolution based on the monthly maximum NDVI from the Sentinel-2 and Landsat data. Figure 5 illustrates the early-season winter cereal distribution for May 2018 for all 32 European countries. Based on the municipal-level agricultural census data, we evaluated the performance of the TWDTW method at the municipal-level for identifying winter cereals. Since the standard NDVI curves for winter cereals were derived from field samples collected in 2018, we examined the accuracy of the method for identifying winter cereals in 2018. The R 2 values between the identified planted area and the agricultural census data ranged from 0.53 to 0.99 for the 11 countries, and the slope of the regression lines between the identified areas and the census areas ranged from 0.6 to 1.35, from 2016 to 2020. Meanwhile, the RE ranged from 9% to 85%, and RMSE ranged from 2 × 10 4 ha to 23 × 10 4 ha, demonstrating that the values of RE and RMSE changed greatly in these 11 countries ( Figure 6). We also examined the potential of this approach for other years (i.e., 2016, 2017, and 2019), and the identified winter cereal areas were compared with the agricultural census area for the other three years (2016, 2017, 2019), except 2020 due to the lack of agricultural census. The results showed that R 2 and the slope changed little compared with 2018 in most countries, except for France and Denmark, which implied that the TWDTW method had the ability to apply the standard NDVI seasonal change curve from one year to other years ( Figure 6).
In addition, we examined the accuracy of the identification at the pixel level based on the 3356 field samples of winter cereals and other land cover types that were obtained from GEE for 32 countries in Europe ( Table 2). The overall accuracies vary from 71.7% to 98.73% in 32 countries, and the Kappa coefficients range from 0.42 to 0.97. The average overall accuracy is above 90% when the 32 countries in Europe are considered together, and the Kappa coefficients are higher than 0.7 in most countries, indicating that the TWDTW method has the potential to identify planted areas of winter cereals in our study area. However, lower identification accuracy still exists in some countries. For example, Bulgaria has the lowest overall accuracy (71.7%) and Kappa coefficient (0.42). In addition, Figure 7 shows zoomed-in images, with rich spatial details for the winter cereals in France and Romania in 2018.
Remote Sens. 2022, 14, x FOR PEER REVIEW 8 of 1 from 0.53 to 0.99 for the 11 countries, and the slope of the regression lines between the identified areas and the census areas ranged from 0.6 to 1.35, from 2016 to 2020. Mean while, the RE ranged from 9% to 85%, and RMSE ranged from 2 × 10 4 ha to 23 × 10 4 ha demonstrating that the values of RE and RMSE changed greatly in these 11 countries (Fig  ure 6). We also examined the potential of this approach for other years (i.e., 2016, 2017 and 2019), and the identified winter cereal areas were compared with the agricultural cen sus area for the other three years (2016, 2017, 2019), except 2020 due to the lack of agricul tural census. The results showed that R 2 and the slope changed little compared with 2018 in most countries, except for France and Denmark, which implied that the TWDTW method had the ability to apply the standard NDVI seasonal change curve from one yea to other years ( Figure 6).  In addition, we examined the accuracy of the identification at the pixel level based on the 3356 field samples of winter cereals and other land cover types that were obtained from GEE for 32 countries in Europe ( Table 2). The overall accuracies vary from 71.7% to 98.73% in 32 countries, and the Kappa coefficients range from 0.42 to 0.97. The average overall accuracy is above 90% when the 32 countries in Europe are considered together, and the Kappa coefficients are higher than 0.7 in most countries, indicating that the TWDTW method has the potential to identify planted areas of winter cereals in our study area. However, lower identification accuracy still exists in some countries. For example, Bulgaria has the lowest overall accuracy (71.7%) and Kappa coefficient (0.42). In addition, Figure 7 shows zoomed-in images, with rich spatial details for the winter cereals in France and Romania in 2018. To examine the effect of SAR data in the identification, we took Germany as an example to identify the winter cereals with and without using SAR data during the implementation of the TWDTW method. The relationship between the estimated and agricultural statistical planted areas of winter cereals with and without SAR data was illustrated in Figure 8. The performance of the TWDTW method with SAR data is better than that without using SAR data. The R 2 value increased from 0.49 to 0.81 after using SAR data.
We also used the random forest method to identify winter cereals with the same field samples for further comparison with the TWDTW method. We compared the accuracy of the TWDTW method and random forest method in the classification of winter cereals in a total of 32 European countries. For a fair comparison, the same training samples were used to run the random forest method in GEE and the TWDTW method in the local server, respectively. Moreover, the same field samples were employed to evaluate the accuracy of classification. Overall, the accuracy of the TWDTW method is higher than random forest in identifying winter cereals ( Table 3). The overall accuracies are 90.82% and 80.84% for TWDTW and random forest, respectively, while the Kappa coefficients are 0.58 and 0.70, respectively. Our results demonstrated that the TWDTW method had higher accuracy than the random forest method in classification.   To examine the effect of SAR data in the identification, we took Germany as an example to identify the winter cereals with and without using SAR data during the implementation of the TWDTW method. The relationship between the estimated and agricultural statistical planted areas of winter cereals with and without SAR data was illustrated in Figure 8. The performance of the TWDTW method with SAR data is better than that without using SAR data. The R 2 value increased from 0.49 to 0.81 after using SAR data. We also used the random forest method to identify winter cereals with the same field samples for further comparison with the TWDTW method. We compared the accuracy of the TWDTW method and random forest method in the classification of winter cereals in  To further examine the details of classification with the TWDTW method and random forest, we selected a small place in Romania to check the classification (Figure 9). Compared with the TWDTW method, random forest had a large misclassification of non-winter cereals, which was identified for winter cereals. For random forest, the overall accuracy is 85.44%, and kappa coefficient is 0.75, while for TWDTW, the overall accuracy is 93.86%, and kappa coefficient is 0.88.  To further examine the details of classification with the TWDTW method and random forest, we selected a small place in Romania to check the classification (Figure 9). Compared with the TWDTW method, random forest had a large misclassification of nonwinter cereals, which was identified for winter cereals. For random forest, the overall accuracy is 85.44%, and kappa coefficient is 0.75, while for TWDTW, the overall accuracy is 93.86%, and kappa coefficient is 0.88.

Discussion
In this study, we used a phenology-based TWDTW method for identifying winter cereals by comparing the similarity of seasonal changes in NDVI at a given pixel with a standard NDVI curve for winter cereals. We produced the distribution maps for winter cereals with 30 m spatial resolution for 32 countries in Europe and the average overall accuracy was 91.8%. We also used the random forest method to identify winter cereals with the same field samples for further comparison with the TWDTW method.
Our results demonstrated that the random forest method had lower accuracy than the TWDTW method in the classification of winter cereals, although random forest is faster in GEE than the TWDTW method in local server. A small number of samples was the main cause of low accuracy for random forest and less vegetation types of non-winter-

Discussion
In this study, we used a phenology-based TWDTW method for identifying winter cereals by comparing the similarity of seasonal changes in NDVI at a given pixel with a standard NDVI curve for winter cereals. We produced the distribution maps for winter cereals with 30 m spatial resolution for 32 countries in Europe and the average overall accuracy was 91.8%. We also used the random forest method to identify winter cereals with the same field samples for further comparison with the TWDTW method.
Our results demonstrated that the random forest method had lower accuracy than the TWDTW method in the classification of winter cereals, although random forest is faster in GEE than the TWDTW method in local server. A small number of samples was the main cause of low accuracy for random forest and less vegetation types of non-winter-cereals directly led to the misclassification that non-winter cereals were identified as winter cereals. Compared with the machine learning algorithms, the TWDTW method requires only a few training samples and is beneficial for large scale crop identification, even if the samples are insufficient [5,46]. A previous study used 30,000 samples as the training dataset to classify the cropland in Europe by using the machine learning method [41]. However, the overall accuracy only ranged from 79.2% to 88.8% in most of the countries. Based on 58,178 input samples, a recent study used the random forest classifier to generate the crop type map in Europe. For wheat, the producer's accuracy and user's accuracy were 78.2% and 49.6%, respectively [47]. Moreover, the standard NDVI seasonal change curves (calculated from the filed samples in 2018) were successfully applied to identify winter cereals for other years (2016, 2017, 2019) by using the TWDTW method. The R 2 between identified areas and agricultural census areas were greater than 0.5 in most of the countries, indicating that the TWDTW method was effective even though there are no training samples in other years [26]. A previous study trained the classifier based on samples from 2008 to 2011 and applied the trained classifier from 1894 to 2007 to identify crops by using the Classification and Regression Tree (CART) algorithm [48]. However, the results showed that the R 2 between identified areas and census areas only ranged from 0.142 to 0.192 in different crops.
The TWDTW method can identify winter cereals two months before harvest. This is an important benefit of the method that makes it useful for the early and continuous prediction of winter cereal production [49]. Knowing the distribution of crops ahead of the harvest, we could determine in a timely manner the regions that may suffer food crisis and provide information for local governments to take action and provide assistance [30,50]. Here, we compared the accuracy of early season maps and end of season maps to further evaluate the performance of our method. The end of season (i.e., July) winter cereal maps meant that all growing season images were used to generate the maps and the overall accuracy was 91.8%. The overall accuracy of early season winter cereal maps (88.4%) was stable in May and only 3.4% lower than that of end of season maps. A recent study [5] has also examined the ability of the TWDTW method for early-season mapping and revealed that winter wheat can be identified at the end of April with an overall accuracy of 89.88%, which is comparable to our results. In addition, some studies have conducted early-season winter crop mapping based on other methods. Skakun et al. [8] used a phenology-based method for winter crop identification and demonstrated that the winter crop maps could be achieved 1.5-2 months before harvest with an overall accuracy of 90%, but the spatial resolution (250 m) was relatively coarse. Tian et al. [51] employed the decision tree classification model to map winter crops in some counties and found that winter crops could be distinguished at an average of three months before harvest with an accuracy greater than 90%. However, this requires a large amount of training samples, which will limit the areas lacking samples and is not conducive for large scale mapping [52]. In addition, previous studies have noted that computing time is one of the challenges when using the TWDTW method [24,26]. In this study, the TWDTW method was run on a configuration using 8 cores, with 2.70 GHz and 512 GB memory. In Finland (the area was approximately 33.8 million ha), the processing time was 38 min when classifying winter cereal in 2018 (10 images). This process is significantly faster than that recorded in the previous study, which consumed 9 h when performing crop mapping with 13 images (4 cores and the area of study region was 24,256 ha) [26]. Therefore, the calculation procedure of the TWDTW method used in this study can provide support for real-time crop mapping. Furthermore, incorporating SAR data into the TWDTW method can efficiently remove the rapeseed from other winter cereals and the resulting classification depended on the rapeseed planting area.
In this study, the accuracy of the TWDTW method is comparable for all countries except Bulgaria, where it has a relatively poor overall accuracy of 71%. This is probably because Bulgaria has a Mediterranean climate in the south and a continental climate in the north, creating variability in the phenological features and growing status of winter cereals. In addition, there are some limitations for the TWDTW method. First, the threshold of Euclidean distance used to judge whether a pixel is a winter cereal relies strongly on statistical agricultural census data. Therefore, the accuracy of the method is closely related to the reliability of agricultural census data. However, the agricultural census data usually have hysteresis, which is still a challenge for real-time crop identification [5,53]. Secondly, the effectiveness of retrieving seasonal changes in crop growth is largely determined by the quantity of cloud-free image data, which directly affects the accuracy of identification. Although we used the monthly maximum merged NDVI values from Landsat and Sentinel-2 to obtain more effective observation data, the TWDTW method performed poorly in identifying winter cereals when the cloud percentage was up to 60%. A previous study has demonstrated that the accuracy of crop identification will be reduced when lacking effective satellite data [54]. In addition, the difference in the wavelength between the Sentinel-2 and Landsat sensors may affect the quality of the synthetic NDVI data. It is still a challenge to reduce the impact from this difference [5,55]. In the future, development of the TWDTW method without the limitation of agricultural census area and a higher applicability of effective satellite data will help to improve the identification accuracy of crops.

Conclusions
In this study, we used the TWDTW method to produce winter cereal maps with 30 m spatial resolution in Europe for the period of 2016 to 2020 by combining Landsat and Sentinel data. The TWDTW method performed well in Europe for winter cereal identification, with an average overall accuracy of 91.8% and producer's and user's accuracies of 91% ± 7.8% and 89% ± 10.3%, respectively. The Kappa coefficients are more than 0.7 in most countries. Compared with the agricultural census area, the TWDTW method explained 77% ± 15% of spatial variability in the planted area at the municipal level for 11 countries. More importantly, the TWDTW method can effectively provide winter cereal maps two months before harvest with stable accuracy. In general, our study mapped the planted area of winter cereals with a high spatial resolution in Europe, which will provide a resource for decision makers to timely monitor the growth of winter cereals.