Developing High-Resolution Crop Maps for Major Crops in the European Union Based on Transductive Transfer Learning and Limited Ground Data

: Precise and timely information on crop spatial distribution over large areas is paramount to agricultural monitoring, food security, and policy development. Currently, automatically classifying crop types at a large scale is challenging due to the scarcity of ground data. Although previous studies have indicated that transductive transfer learning (TTL) is a promising method to address this problem, it performs poorly within regions where crop compositions and phenology differ largely. Here we transferred random forest classiﬁers trained in limited regions with diversiﬁed growing conditions and land covers to the rest of the study area where ground data are scarce, with more than 130,000 Sentinel-2 images processed using the Google Earth Engine (GEE) platform. We established the 10 m crop maps for four major crops (i.e., maize, rapeseed, winter, and spring Triticeae crops) across 10 European Union (EU) countries from 2018 to 2019. The ﬁnal crop maps had a high accuracy with overall accuracy generally greater than 0.89, with user’s accuracy and producer’s accuracy ranging from 0.72 to 0.98. Moreover, the resulting maps were consistent with the NUTS-2 level ofﬁcial statistics, with R2 consistently greater than 0.9. We further analyzed the crop rotation patterns and found that the rotation intervals across these EU countries were generally at least one year. Maize was dominantly rotated with winter Triticeae crops or converted to other land covers in the following year. Rapeseed was generally grown in rotation with winter Triticeae crops, whereas the rotation patterns of winter and spring Triticeae crops were more diversiﬁed. Red Edge Position (REP) and Normalized Difference Yellow Index (NDYI) played signiﬁcant roles in crop classiﬁcation across the EU. This study highlights the potential of the developed TTL method for crop classiﬁcation over large spatial extents where labeled data are limited and the differences in crop compositions and phenology are relatively large.


Introduction
Unprecedented global population growth led to a continual increase in food demand and agricultural production intensification [1,2]. However, land degradation, shortages of cropland and water resources, and more frequent extreme weather events caused by global warming have hampered the increase in food production and threatened food security [3][4][5]. Accurate and near-real-time monitoring of crop spatial distribution over large areas has important implications for crop growth monitoring, agricultural land use and agronomic management optimization, and food security early warning [6].
With an enormous increase in the spatial and temporal information derived from satellites, remote sensing has been widely used in agricultural applications for decades [7][8][9][10]. Previous studies have primarily used the Normalized Difference Vegetation Index (NDVI) land cover map for Germany [17,37], maps of the major crops in Northeast China [38], and rice paddy maps in Bangladesh [39]. Therefore, datasets of high-resolution crop maps for the EU are particularly needed to facilitate the development of related scientific research on food security and policy making.
The objective of this study is to (1) propose an automatic approach to identifying staple crops with limited ground samples by transferring an RF classifier, based on Sentinel-2 imagery through the Google Earth Engine (GEE) cloud computing platform; (2) demonstrate the reliability of this approach in mapping maize, rapeseed, winter Triticeae crops (winter wheat, barley, and rye), and spring Triticeae crops (spring wheat, barley, and oat) from 2018 to 2019 across the major cultivation areas in the EU.

Study Area
The study area included 10 countries, i.e., England, the Netherlands, Germany, Denmark, France, Italy, Poland, Hungary, Slovakia, and the Czech Republic, covering approximately 195.8 Mha of land area across the EU (Figure 1). Triticeae crops, maize, and rapeseed are staple crops in the EU, accounting for 65.9% of the harvested area for primary crops [40]. We chose these countries for two reasons: (1) To cover the major planting areas of the EU. In particular, the harvested area of the above crops in these countries makes up a total of 62.7% across the EU, and in each country accounts for beyond 2%, except for the Netherlands, which is selected due to its detailed crop field data; (2) To span a widespread growing environment where the scalability and reliability of our methods are investigated, consisting of dry, temperate, continental, and polar and alpine climates [41].

Reference Data
Reference data were acquired from existing nationwide crop field datasets or land cover maps. The first data for England, provided by the Rural Payment Agency (https://ckan.publishing.service.gov.uk/dataset, accessed on 7 April 2022) and named crop map of England (CROME), were generated by an RF classifier combining Sentinel-1 and S2 images with an overall accuracy (OA) of 95.4% (81%) and a Kappa coefficient of 0.95 (0.8) for 2019 (2018). CROME included over 20 main crop types, grassland, and non-agricultural land coves (e.g., woodland and water), and the user's accuracy (UA) and producer's accuracy (PA) of major crop types were generally greater than 80%. The second data, 10 m land cover maps for France (https://www.theia-land.fr/en/product/land-cover-map, accessed on 7 April 2022), classified 23 land cover classes since 2018 based on S2 images by an RF classifier with an OA of 89% and F1 scores for crop classes exceeding 90% [42]. The third dataset was the crop field dataset retrieved from the Base Registration Crop Parcels (BRP) in the Netherlands, providing the cultivated crop information at the parcel level (i.e., polygon) in Agricultural Area Netherlands (www.PDOK.nl, accessed on 7 April 2022). In addition, the Land Use/Cover Area frame Survey (LUCAS) in situ data was used to directly validate the classification results for all EU countries in 2018, conducted every three years since 2006 [43,44]. The previous LUCAS core observation (i.e., a circle with a 1.5 m radius and area of 7.07 m 2 ) that corresponds to a fraction of an S1/2 pixel cannot be directly useful for such decametric sensors. The up-to-date LUCAS 2018 data delineated a polygon of 0.52 ha with homogeneous land cover by a new LUCAS Copernicus module, making it suitable for verifying high-resolution land cover maps.
The land cover maps for England and France produced by the RF classifier necessitated quality control to ensure the high-quality training and validation of the resultant crop maps. They were both used for sampling as more widespread ground samples covering diverse growing conditions over space perhaps made large-area crop mapping using TTL successful. The classification results in CROME were represented as hexagonal cells of 0.41 ha (approximate to a 64 m grid); however, the predicted confidence information for each cell was lacking. Thus, we first generated a 12 km-by-12 km grid across England and then selected parcels larger than 1 ha for each land cover, which corresponds to 3 cells with the same label. Finally, we randomly sampled 250 points in each grid cell and filtered out the points falling outside the parcels. Such sampling method was repeated for land cover maps of France; the differences were that (1) a high confidence mask with a threshold of 95% was used to restrict the sampling extent; (2) the "straw cereals" class was excluded from sampling as it contained both winter and spring Triticeae crops and was not directly usable; and (3) a 50 km grid cell was generated across France. In total, 56,877 and 55,144 training samples were acquired for 2018 and 2019, respectively, containing five land cover classes including other land covers (e.g., sugar beet, potato, sunflower, grassland), maize, rapeseed, winter Triticeae crops (i.e., winter wheat, barley, and rye) and spring Triticeae crops (i.e., spring wheat, barley, and oat) ( Table 1). The cultivation area of each crop type at a national (i.e., Nomenclature of Territorial Units for Statistics-1 regions, NUTS-1 regions) and subnational level (i.e., NUTS-2 regions) from 2018 to 2019 was acquired from Eurostat (https://ec.europa.eu/eurostat/databrowser, accessed on 7 April 2022). Such data were used as additional validation of the resultant maps, particularly in areas where ground samples were not available.

Methods
The flowchart for TTL-based crop mapping is presented schematically in Figure 2, consisting of image preprocessing, feature selection, transferring supervised classification, and accuracy assessment. The cultivation area of each crop type at a national (i.e., Nomenclature of Territorial Units for Statistics-1 regions, NUTS-1 regions) and subnational level (i.e., NUTS-2 regions) from 2018 to 2019 was acquired from Eurostat (https://ec.europa.eu/eurostat/databrowser, accessed on 7 April 2022). Such data were used as additional validation of the resultant maps, particularly in areas where ground samples were not available.

Methods
The flowchart for TTL-based crop mapping is presented schematically in Figure 2, consisting of image preprocessing, feature selection, transferring supervised classification, and accuracy assessment.

Image Processing
We first used the Landsat simple cloud score algorithm to mask out low-quality pixels due to cloud contamination [45]. Such a method assumed that: (1) clouds are reasonably bright in the blue and cirrus bands; (2) clouds are reasonably bright in all visible bands; (3) clouds are moist; (4) clouds are not snow. It calculated cloud score based on four spectral bands (Aerosols, Blue, Green, and Red band) and two spectral indices (Normalized Difference Moisture Index (NDMI) and Normalized Difference Snow Index (NDSI)). Clouds were detected more thoroughly using this method than using the QA60 band which could only mask out dense and cirrus clouds [18]. Moreover, cloud shadows were masked using the Temporal Dark Outlier Mask (TDOM) method [46]. Then, we reconstructed the S2 time series by a moving median composite method as the number of clear observations showed obvious spatiotemporal heterogeneity [18]. In particular, we first generated a time series of 15-day temporal interval composites by calculating the median

Image Processing
We first used the Landsat simple cloud score algorithm to mask out low-quality pixels due to cloud contamination [45]. Such a method assumed that: (1) clouds are reasonably bright in the blue and cirrus bands; (2) clouds are reasonably bright in all visible bands; (3) clouds are moist; (4) clouds are not snow. It calculated cloud score based on four spectral bands (Aerosols, Blue, Green, and Red band) and two spectral indices (Normalized Difference Moisture Index (NDMI) and Normalized Difference Snow Index (NDSI)). Clouds were detected more thoroughly using this method than using the QA60 band which could only mask out dense and cirrus clouds [18]. Moreover, cloud shadows were masked using the Temporal Dark Outlier Mask (TDOM) method [46]. Then, we reconstructed the S2 time series by a moving median composite method as the number of clear observations showed obvious spatiotemporal heterogeneity [18]. In particular, we first generated a time series of 15-day temporal interval composites by calculating the median value of observations at each interval. The reason for choosing the 15-day composite was that there was no substantial shift in plant growth within two weeks and the overall classification performance was similar between 10-day and 15-day composites [18,47]. Second, a gap-filling method was conducted on the composite time series in which a given observation was replaced by the median value of three adjacent observations (i.e., previous, current, and subsequent observations). Such a method is readily realized in the GEE platform and can generate homogenous results with high computational efficiency [48]. It is also not affected by outliers in time series compared to the mean value composite method.

Transferring Supervised Classification
We used RF classifier to map crop extent, which has been widely applied in land cover classification across different geographical areas [24,37,48,51]. Moreover, previous studies have demonstrated that RF was immune to high data dimensionality and multicollinearity and robust to noisy training data [7,52]. RF is an ensemble classifier consisting of multiple decision trees in which two-thirds of the training samples are trained and the remaining one-third is used for cross-validation [53]. The final classification output is selected as the class with the maximum votes against all decision trees. The number of decision trees (numberOfTrees) was set at 100 according to the trade-off between the classification accuracy and computational efficiency, while other parameters were set to the default in the GEE platform.
We first trained the RF classifier on training samples derived from the land cover maps of England and France (mentioned in Section 2.2.2) using a total of 126 spectral features (14 composite periods with 9 features in each period) ( Table 2). Our supervised models classified land covers into five classes, i.e., "others", "maize", "rapeseed", "winter Triticeae crops", and "spring Triticeae crops". We tried a range of time windows (e.g., 1 January to 31 December or 1 February to 31 December) and determined to select 1 April

Transferring Supervised Classification
We used RF classifier to map crop extent, which has been widely applied in land cover classification across different geographical areas [24,37,48,51]. Moreover, previous studies have demonstrated that RF was immune to high data dimensionality and multicollinearity and robust to noisy training data [7,52]. RF is an ensemble classifier consisting of multiple decision trees in which two-thirds of the training samples are trained and the remaining one-third is used for cross-validation [53]. The final classification output is selected as the class with the maximum votes against all decision trees. The number of decision trees (numberOfTrees) was set at 100 according to the trade-off between the classification accuracy and computational efficiency, while other parameters were set to the default in the GEE platform.
We first trained the RF classifier on training samples derived from the land cover maps of England and France (mentioned in Section 2.2.2) using a total of 126 spectral features (14 composite periods with 9 features in each period) ( Table 2). Our supervised models classified land covers into five classes, i.e., "others", "maize", "rapeseed", "winter Triticeae crops", and "spring Triticeae crops". We tried a range of time windows (e.g., 1 January to 31 December or 1 February to 31 December) and determined to select 1 April to 31 October considering the trade-off between the number of good satellite observations and classification accuracy. The model was then transferred to other countries where labeled data were lacking, to test its spatial scalability and robustness in large-area crop mapping with limited training data.

Accuracy Assessment
First, point-based validation was implemented in two ways. One was performed based on crop field datasets or land cover maps in England, France, and the Netherlands for two years, and the other was based on LUCAS in situ data covering the whole studied areas only for 2018 due to data unavailability in 2019 (Table 2). In particular, we first adopted a stratified random sampling design, 5000 samples were randomly distributed across each country and the number of each class's samples was determined according to the areal ratio in each country [54]. We used a sampling unit consisting of 3 × 3 pixels to mitigate the adverse impact of geolocation inaccuracy and then labeled these samples according to reference data. The samples located in homogeneous parcels and larger than 1 ha or corresponding confidence greater than 95% were retained to ensure high-quality validation. Overall, 4953 (4968), 4991 (4984), and 2620 (2634) validation samples were reserved in 2018 (2019) for England, France, and the Netherlands, respectively. Additionally, the LUCAS Copernicus data were linked with LUCAS core points to augment ground truth information based on the conditions that the LUCAS theoretical geolocation fell within the LUCAS polygon and the reported Copernicus level-2 land cover was consistent. As a result, 28,882 LUCAS Copernicus points were retained for validation. Finally, the accuracy was summarized by the confusion matrix, including four metrics (i.e., PA, UA, OA, and F1 score) which were calculated as: where N c is the total number of correctly classified validation samples, n represents total number of validation samples, N i denotes the number of correctly identified validation samples of class i, R i donates the number of validation samples of class i, C i implies the number of validation samples classified as class i.
In addition, area comparison using statistical data was a commonly used method for accuracy assessments on large scales [12,35,37]. The area estimates for each class from the resultant crop maps were compared to the NUTS-2 level statistics, except for spring Triticeae crops that were only recorded at the NUTS-1 level. The coefficient of determination (R 2 ) and the root mean square error (RMSE) were calculated to evaluate the consistency in the area between classification maps and statistical data, which were defined as: where O i and E i are the statistical and estimated area of a NUTS-2 (NUTS-1) geopolitical entity i, respectively. n represents the total number of countries.

Accuracy Assessment of the Developed TTL-Based Method
To examine whether RF models trained in limited regions with plentiful training data could be applied to classify crops over large areas where ground data were lacking, we conducted rigorous accuracy assessments combining both point-based validation and subnational-scale area comparison (Figures 4 and 5, Tables 3, 4 and S2-S8). In total, the OA was generally greater than 0.89, whereas PA and UA varied between different countries or crops. The crop maps for three countries with abundant validation samples showed higher OA ranging from 0.93 to 0.96. The best-performing crops were winter Triticeae crops in England and maize in the Netherlands with PA and UA above 0.9. Winter Triticeae crops generally showed higher accuracy with PA and UA > 0.88 than spring Triticeae crops (PA and UA > 0.7). Rapeseed achieved moderately high accuracy (PA and UA > 0.84). Since the land cover map for France and the LUCAS data did not distinguish between winter and spring Triticeae crops, we defined them as a single class named Triticeae crops. For France, Triticeae crops were mapped with high accuracy with an F1 score of 0.9, followed by rapeseed and maize (F1 ≥ 0.86). Overall, validating the crop map for 2018 using the LUCAS data covering the whole study area indicated that the classification performance was the best for maize with an F1 score of 0.81, followed by rapeseed and Triticeae crops (F1 ≥ 0.77).
Triticeae crops generally showed higher accuracy with PA and UA > 0.88 than spring Triticeae crops (PA and UA > 0.7). Rapeseed achieved moderately high accuracy (PA and UA > 0.84). Since the land cover map for France and the LUCAS data did not distinguish between winter and spring Triticeae crops, we defined them as a single class named Triticeae crops. For France, Triticeae crops were mapped with high accuracy with an F1 score of 0.9, followed by rapeseed and maize (F1 ≥ 0.86). Overall, validating the crop map for 2018 using the LUCAS data covering the whole study area indicated that the classification performance was the best for maize with an F1 score of 0.81, followed by rapeseed and Triticeae crops (F1 ≥ 0.77).     The comparisons showed that the estimated and statistical areas for the winter Triticeae crops, maize, and rapeseed were evenly distributed around the 1:1 line, whereas those for the spring Triticeae crops appeared to be more scattered, and the corresponding mapped area exhibited a slight underestimation, especially in 2018 compared with the other three crops ( Figure 5). The area estimates were consistent with the statistics with R 2 The comparisons showed that the estimated and statistical areas for the winter Triticeae crops, maize, and rapeseed were evenly distributed around the 1:1 line, whereas those for the spring Triticeae crops appeared to be more scattered, and the corresponding mapped area exhibited a slight underestimation, especially in 2018 compared with the other three crops ( Figure 5). The area estimates were consistent with the statistics with R 2 consistently greater than 0.9 and RMSE less than 126.39 Kha, indicating that our TTL-based method succeeded in large-area crop mapping where ground data were limited.

Crop Type Classification Results
The resultant crop maps for 10 EU countries in 2018 and 2019 at 10 m spatial resolutions were presented in Figures 6 and S1. Overall the outlines of individual agricultural parcels were distinctly identified in the current maps, indicating that both inter-field heterogeneity and intra-field homogeneity were well characterized. The agricultural lands in the EU were dominated by medium-sized fields (2-25 ha) and the crop maps captured such diversified agricultural parcels fairly well. Although the field size was relatively small (even < 0.2 ha) in the Netherlands, a 10 m pixel was enough and the classification results clearly delineated these small parcels.
We further calculated the important values of the features in the RF models to investigate the critical factors for crop mapping across EU countries ( Figure S2). NDVI during August was consistently ranked as the top feature, followed by REP, red edge 1, NDYI, and swir1 in unfixed order, highlighting the important roles of the REP, red edge bands, and NDYI in classifying the major crop types across EU countries. parcels were distinctly identified in the current maps, indicating that both inter-field heterogeneity and intra-field homogeneity were well characterized. The agricultural lands in the EU were dominated by medium-sized fields (2-25 ha) and the crop maps captured such diversified agricultural parcels fairly well. Although the field size was relatively small (even < 0.2 ha) in the Netherlands, a 10 m pixel was enough and the classification results clearly delineated these small parcels. We further calculated the important values of the features in the RF models to investigate the critical factors for crop mapping across EU countries ( Figure S2). NDVI during August was consistently ranked as the top feature, followed by REP, red edge 1, NDYI, and swir1 in unfixed order, highlighting the important roles of the REP, red edge bands, and NDYI in classifying the major crop types across EU countries.
A visual comparison of the current crop maps and the reference dataset for France and the Netherlands revealed that crop spatial distributions were correctly captured across numerous varieties of agricultural landscapes (Figures 7 and S3). Our products outperformed existing land cover maps for England and France and filled the gaps in highresolution crop maps for EU countries. One drawback of CROME was that the classification results were represented as hexagonal cells of 0.41 ha (approximate to a 64 m grid), which could contain several land cover types. By contrast, our results successfully made clear distinctions among the different crops. Furthermore, we distinguished winter and spring Triticeae crops even though the land cover maps for France incorporated them as one class. A visual comparison of the current crop maps and the reference dataset for France and the Netherlands revealed that crop spatial distributions were correctly captured across numerous varieties of agricultural landscapes (Figures 7 and S3). Our products outperformed existing land cover maps for England and France and filled the gaps in high-resolution crop maps for EU countries. One drawback of CROME was that the classification results were represented as hexagonal cells of 0.41 ha (approximate to a 64 m grid), which could contain several land cover types. By contrast, our results successfully made clear distinctions among the different crops. Furthermore, we distinguished winter and spring Triticeae crops even though the land cover maps for France incorporated them as one class.
Winter Triticeae crops were the most widely cultivated crops across the whole study area (accounting for 44.3% of the total four crop areas), followed by maize (26.4%), rapeseed (15.6%), and spring Triticeae crops (13.7%). They showed similarity in the spatial pattern along latitudes and were mainly distributed with latitudes ranging from 47-56 • N ( Figure S4a-d). More specifically, winter Triticeae crops were also cultivated south of 47 • N, whereas spring Triticeae crops were rarely planted, and maize was mostly distributed with latitudes ranging from 47 to 49 and 51 to 54 • N (roughly accounting for 54% of the total maize area). However, these crops differed largely in the spatial distribution along longitudes. Winter Triticeae crops were widely distributed with longitudes ranging from 2 • W to 5 • E and 8 to 21 • E, whereas the distribution of spring Triticeae crops was more concentrated with longitudes ranging from 8 to 12 • E and 15 to 23 • E ( Figure S4e,h). The distribution of maize was relatively widespread with longitudes ranging from 3 • W to 2 • E, 5 to 14 • E, and 16 to 23 • E. Rapeseed was distributed at longitudes of 2 • W-5 • E and 9-20 • E ( Figure S4f,g).

Crop Rotation Analysis
The rotation patterns of each crop were further analyzed based on high-accuracy crop maps from 2018 to 2019 (Figures 8-10). The rotation pattern in a certain crop pixel was considered as no rotation (≥1-year break) if the crop was identified in the pixel in both  (Figure 8). Rapeseed was the most remarkable crop with all the countries consistently characterized by a rotation break of ≥1 year in more than 99% of pixels, followed by spring Triticeae crops (90.9%), winter Triticeae crops (84.1%), and maize (78.7%) (Figure 8). Denmark had the lowest proportion of pixels (<80%) with a rotation for more than one year. Moreover, a relatively less proportion of rotation was found in France, Italy, the Netherlands, and Poland for maize, and Germany and England for winter Triticeae crops. The spatial details of the rotation patterns further demonstrated the differences in rotation patterns among different crops (Figure 9). Interestingly, the crop parcels from the previous year were distributed in adjacent fields in the following year ( Figure 9). Winter Triticeae crops were the most widely cultivated crops across the whole study area (accounting for 44.3% of the total four crop areas), followed by maize (26.4%), rapeseed (15.6%), and spring Triticeae crops (13.7%). They showed similarity in the spatial pattern along latitudes and were mainly distributed with latitudes ranging from 47-56°N ( Figure S4a-d). More specifically, winter Triticeae crops were also cultivated south of 47°N, whereas spring Triticeae crops were rarely planted, and maize was mostly distributed with latitudes ranging from 47 to 49 and 51 to 54°N (roughly accounting for 54% of the total maize area). However, these crops differed largely in the spatial distribution along longitudes. Winter Triticeae crops were widely distributed with longitudes ranging from 2°W to 5°E and 8 to 21°E, whereas the distribution of spring Triticeae crops was more concentrated with longitudes ranging from 8 to 12°E and 15 to 23°E (Figure S4e,h). The distribution of maize was relatively widespread with longitudes ranging from 3°W to 2°E, 5 to 14°E, and 16 to 23°E. Rapeseed was distributed at longitudes of 2°W-5°E and 9-20°E mark had the lowest proportion of pixels (<80%) with a rotation for more than one year. Moreover, a relatively less proportion of rotation was found in France, Italy, the Netherlands, and Poland for maize, and Germany and England for winter Triticeae crops. The spatial details of the rotation patterns further demonstrated the differences in rotation patterns among different crops (Figure 9). Interestingly, the crop parcels from the previous year were distributed in adjacent fields in the following year (Figure 9). We also investigated the land cover type used for crop rotation in the following year by calculating the proportions of pixels for a certain crop in 2018 converted to other land cover types in 2019 ( Figure 10). Maize was generally grown in rotation with winter Triticeae crops (50.7% of the total pixels had undergone this rotation) or converted to other land covers (40.7%), whereas the rotation of rapeseed and winter Triticeae crops was the dominant pattern for rapeseed (86.5%). The rotation patterns of the winter and spring Triticeae crops were more complicated. The former was mainly rotated with other land covers (44.1%), rapeseed (23.8%), maize (20.5%), and a little with the latter (11.6%). On the contrary, spring Triticeae crops were primarily rotated with winter Triticeae crops (44.6%), followed by other land covers (31.7%) and rapeseed (14.9%). Additionally, the land cover type used for crop rotation (except for rapeseed) varied between countries. In Czechia and Denmark, maize was also rotated with spring Triticeae crops (>27%).Unlike the predominant rotation of maize and winter Triticeae crops in Germany, France, and England (>66%), maize was grown in rotation with other land covers in Hungary, Italy, and the Netherlands (>69%). Moreover, the winter Triticeae crops were largely rotated with the spring Triticeae crops in Denmark (>45%). In addition, the rotation patterns of the three crops other than rapeseed in the Netherlands were consistent, that is, dominantly rotated with other land covers (>77%).     We also investigated the land cover type used for crop rotation in the following year by calculating the proportions of pixels for a certain crop in 2018 converted to other land cover types in 2019 ( Figure 10). Maize was generally grown in rotation with winter Triticeae crops (50.7% of the total pixels had undergone this rotation) or converted to other land covers (40.7%), whereas the rotation of rapeseed and winter Triticeae crops was the dominant pattern for rapeseed (86.5%). The rotation patterns of the winter and spring Triticeae crops were more complicated. The former was mainly rotated with other land covers (44.1%), rapeseed (23.8%), maize (20.5%), and a little with the latter (11.6%). On the contrary, spring Triticeae crops were primarily rotated with winter Triticeae crops (44.6%), followed by other land covers (31.7%) and rapeseed (14.9%). Additionally, the land cover type used for crop rotation (except for rapeseed) varied between countries. In Czechia and Denmark, maize was also rotated with spring Triticeae crops (>27%). Unlike the predominant rotation of maize and winter Triticeae crops in Germany, France, and England (>66%), maize was grown in rotation with other land covers in Hungary, Italy, and the Netherlands (>69%). Moreover, the winter Triticeae crops were largely rotated with the spring Triticeae crops in Denmark (>45%). In addition, the rotation patterns of the three crops other than rapeseed in the Netherlands were consistent, that is, dominantly rotated with other land covers (>77%).

Advantages of the TTL-Based Method for Large-Area Crop Mapping
Although a large number of satellites imageries are increasingly accessible with unprecedented temporal, spatial, and spectral frequency, the availability of ground data did not keep up, hampering the accurate training and validation of commonly used supervised classification methods [22,24,26,55]. Since large-area crop mapping is paramount for ensuring food security and supervised classification methods are predominant tools for realizing such a mission, it is necessary to investigate the feasibility of transferring a supervised model trained in limited regions with abundant ground data to large spatial extents with scarce data. The 10 m crop maps for the 10 EU countries exhibited high accuracy and clearly identified the field boundaries across various agricultural landscapes, suggesting that our TTL-based method is a reliable and promising method to apply to crop classification over large areas with limited ground data. Nevertheless, some points need to be noted: (1) the training samples should cover different growing conditions as reasonably as possible because the phenological differences are relatively large across different regions, and (2) the training samples are required to include plentiful land cover types, especially those that were absent in sampling regions but widespread in other areas. This is well supported by the findings of Kluger et al. (2021) [55], in that it is better to collect training samples used for TTL-based crop classification in a region with diverse crop-type composition than to collect training samples in a region with low crop-type diversity. Previous studies used growing degree days (GDD) to represent the differences in crop phenology and indicated that lower similarity of GDD in two regions could result in poorer performance of an RF model transferred from one to the other [24]. Given our study area, maize phenology is advanced with latitude decreasing ( Figure S5). Thus, how to reasonably cope with such intra-class phenology differences is essential to achieving sound TTL-based crop mapping over large areas. To avoid omission of maize fields that differed from training regions in phenology, we randomly sampled a small number of maize samples from land cover products from France as France was characterized by diversified climatic (e.g., temperate, continental, and polar and alpine climates) characteristics, which to a large extent contain the climate types of the rest of study area (except for dry climate) [41]. This procedure allowed for the capture of phenology differences in other regions such as Italy and Germany, whereas merely using training samples generated in England cannot work. Moreover, the CROME data for England did not contain some types of land cover that frequently occur in other countries such as soybean, sunflower, rice, and mineral surfaces, which are readily misclassified as desired crop types. Therefore, we also collected samples for these classes to augment the diversity of the training sets.
In short, the TTL-based method makes ample use of ground data acquired from limited regions to classify crop types over large spatial extents where labeled data are not available. Compared with previous studies, we first attempt to reasonably enrich the diversity of the training sets to overcome the disadvantages of the TTL method that performs poorly in settings where the differences in crop composition and phenology are relatively large.

Detailed Investigation of Crop Rotation Patterns
The rigorous validation by both point-based validation and subnational-scale area comparison demonstrated high accuracy and reliability of the resultant 10 m crop maps during 2018 and2019. Therefore, it is possible to further investigate the rotation patterns of the four crops, including the rotation interval and land cover types. We found that the crop rotation interval of these EU countries was characterized by at least 1 year, which was consistent with actual situations [56]. Crops are traditionally grown in rotation in a large number of European countries, and it is recognized that (1) most crops benefit from the nutrients released by the residues of the preceding crop [57] and (2) crop rotation contributes to boosting yield, ameliorating soil fertility, and reducing pests and diseases, and consequently improving the resilience of the crop production system under climate change [58][59][60]. In addition, we revealed that maize was generally rotated with winter Triticeae crops or converted to other land covers. Rapeseed was dominantly grown in rotation with winter Triticeae crops, whereas the rotation patterns of winter and spring Triticeae crops were diversified. These findings were supported by some field experiments in Europe [56,[61][62][63]. We provided the refined spatial characteristics of crop rotation patterns across the EU at the parcel level. Such information is novel and essential for many purposes such as agro-ecosystem modeling, land use, and cropping system optimization.

Uncertainty
Although the resultant crop maps have been proven to be reliable, there are still some uncertainties. First, the largest constraint is the availability of S2 data which is primarily restricted by clouds. Even though the moving median composite method partly fills the gaps due to cloud contamination, it continues to malfunction, especially in regions (e.g., Northern Europe) with dense cloud coverage, resulting in a high omission error. This issue is expected to be solved using Harmonized Landsat and Sentinel-2 (HLS) data, which have been applied to crop mapping [26,37,64]. However, such data have not been made available in the GEE platform, hampering rapid and highly efficient crop mapping over large areas. Second, the TOA reflectance product was used rather than surface reflectance (SR) owing to the limited temporal availability of the latter on the GEE platform for our study area. Although previous studies have substantiated that TOA reflectance products are reliable for crop classification [18,49], it might dampen the classification accuracy because of the pronounced scattering and absorption of visible radiation by atmospheric contamination [65]. Here we not only excluded shorter-wavelength bands (e.g., red and blue bands) that were most susceptible to atmospheric effects from feature selection but also used indices based on spectral band ratios to reduce such effects [50,66]. With more S2 SR data available for the study area since 2019, it is promising to improve the classification accuracy in subsequent years. Third, we did not collect training samples for winter and spring Triticeae crops covering different environmental conditions and such ground data were lacking in regions other than England. It not only led to an underestimation of area compared with the official statistics but also restricted the scalability of our TTL-based method. For instance, winter wheat phenology is advanced with latitude decreasing, especially in Spain where the anthesis date can be approximately 40 days earlier than that in England [67,68], hampering the application of the TTL method in such regions with crops intensively cultivated due to the scarcity of corresponding ground data. Fourth, we grouped winter and spring Triticeae crops as single classes, respectively. Nevertheless, there are slight differences in phenology among wheat, barley, rye, and oat. Further researches should attempt to apply the TTL method to refine crop classification (both in crop category and spatial resolution) over larger spatial extents using increasing and precise ground data.

Conclusions
We produced the crop maps at 10 m spatial resolution for four main crops across 10 EU countries from 2018 to 2019 through the GEE platform, by transferring an RF model trained in limited areas with plentiful ground data to other larger extents with scarce ground samples. Compared with previous efforts to classify crops using the TTL method, we reasonably augmented the diversity of training data by additional sampling covering different growing conditions and land cover types, addressing the poor performance of the TTL method when the differences in crop composition and phenology are relatively large. We also demonstrated the important roles of REP and NDYI in crop mapping across the EU. The overall accuracy of the final maps was generally greater than 0.89, with PA and UA ranging from 0.72 to 0.98. In comparison with the NUTS-2 level official statistics, the resultant maps proved to be highly accurate with R 2 consistently greater than 0.9. These rigorous validations indicated that the TTL method was robust in classifying crops over large spatial extents where ground data are limited. The crop rotation analyses showed that the rotation interval across these EU countries was predominantly at least one year. We found that maize was generally rotated with winter Triticeae crops or converted to other land covers in the following year. Rapeseed was grown in rotation with winter Triticeae crops, whereas the rotation patterns of winter and spring Triticeae crops were complicated. The TTL method used in this study can be readily replicated to classify crop types in other large regions, given the strong computing power and ever-increasing satellite observations in the GEE platform. The datasets (available at https://doi.org/10.17632 /4frtdtxsk9.2, accessed on 1 March 2022) are crucial for many purposes including agro-ecosystem modeling, crop growth monitoring, land use and cropping system optimization, and policy making.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs14081809/s1, Figure S1: Crop maps at 10 m resolution of ten EU countries for 2018. Figure S2: Feature importance values for the tope of 50 variables from random forest models in 2018 (a) and 2019 (b). Figure S3: Visual comparison between our final crop maps (a1-a8) and existing reference datasets in 2018 for France (b1-b6) and Netherlands (b7,b8). For France, winter and spring Triticeae crops was concluded as a single category (green color). Figure S4: Distribution of crop area along latitude gradients (a-d) and longitude gradients (e-h). Figure S5: Temporal profiles of NDVI for maize of different countries. Lines depicts the mean values. Shaded area depicts error bars with one positive/negative standard deviation of maize. Table S1: The number of Sentinel-2 images processed in this study for each country in 2018 and 2019.

Data Availability Statement:
The data presented in this study are openly available at https://doi. org/10.17632/4frtdtxsk9.2.

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