The RapeseedMap10 database: annual maps of rapeseed at a spatial resolution of 10 m based on multi-source data

Large-scale, high-resolution maps of rapeseed (Brassica napus L.), a major oilseed crop, are critical for predicting annual production and ensuring global energy security, but such maps are still not freely available for many areas. In this study, we developed a new pixel- and phenology-based algorithm and produced a new data product for rapeseed planting areas (2017–2019) in 33 countries at 10 m spatial resolution based on multiple data. Our product is strongly consistent at the national level with official statistics of the Food and Agricultural Organization of the United Nations. Our rapeseed maps achieved F1 spatial consistency scores of at least 0.81 when compared with the Cropland Data Layer in the United States, the Annual Crop Inventory in Canada, the Crop Map of England, and the Land Cover Map of France. Moreover, F1 scores based on independent validation samples ranged from 0.84 to 0.91, implying a good consistency with ground truth. In almost all countries covered in this study, the rapeseed crop rotation interval was at least 2 years. Our derived maps suggest, with reasonable accuracy, the robustness of the algorithm in identifying rapeseed over large regions with various climates and landscapes. Scientists and local growers can use the freely downloadable derived rapeseed planting areas to help predict rapeseed production and optimize planting structures. The product is publicly available at https://doi.org/10.17632/ydf3m7pd4j.3 (Han et al., 2021).

1 Introduction

Although fossil fuels are currently the main source of energy (Fang et al., 2016; Shafiee and Topal, 2009), their overexploitation is increasing various threats to human survival, such as greenhouse gas emission and environmental pollution (Fang et al., 2016; Höök and Tang, 2013). Biofuel energy seems to be a promising alternative energy source (Hassan and Kalam, 2013). Rapeseed (Brassica napus L.) is an important source of biofuels, edible oil, animal feed, and plant protein powder (Firrisa et al., 2014; Malça and Freire, 2009; Sulik and Long, 2016). Data products on rapeseed planting densities, growth conditions, and productivity are dependent on precise and accurate planting area maps (Zhang et al., 2019), but such maps are still unavailable.

Global agricultural statistics on rapeseed in many regions are derived from field surveys, field sampling, and producer reports (Arata et al., 2020; Fuglie, 2010). Ground-based methods, which are time-consuming and labor-intensive, fail to provide detailed spatial information on rapeseed fields (J. Wang et al., 2020). In contrast, remote sensing technology plays an important role in agricultural monitoring and yields accurate, objective spatial–temporal crop information (Dong et al., 2016; Salmon et al., 2015).

Many current land cover products obtained by remote sensing have a publicly available cropland layer. Examples include the Fine Resolution Observation and Monitoring of Global Land Cover project (Gong et al., 2013), the Global Land Cover 2000 (GLC2000) map (Bartholomé and Belward, 2005), ChinaCropPhen1km (Luo et al., 2020), and Global Food Security-support data at 30 m (GFSAD30) (Phalke et al., 2020; Xiong et al., 2017). Nevertheless, cropland identified by these products is either undifferentiated as to crop type, has a coarse spatiotemporal resolution (Teluguntla et al., 2018), or does not include rapeseed. Few large-scale rapeseed maps, especially at 10 m resolution, are currently available. A decision tree classification method based on a large number of training samples has been used to classify various crops for the 30 m resolution Cropland Data Layer (Boryan et al., 2011) in the USA and the Annual Crop Inventory in Canada (Fisette et al., 2013), but this method is hard to apply to other developing regions because ground training samples are lacking (Xiong et al., 2017). A new method to map large-scale annual maps with a high spatial resolution that would be widely applicable to regions with few ground training samples is thus strongly needed.

Five remote-sensing-based methods for rapeseed mapping have been developed in recent decades: (1) machine learning methods (Griffiths et al., 2019; Preidl et al., 2020; She et al., 2015; Tao et al., 2020), (2) classification based on time series data (Ashourloo et al., 2019), (3) threshold segmentation based on phenology (Tian et al., 2019), (4) multi-range spectral feature fitting (Pan et al., 2013), and (5) mapping based on HSV (hue, saturation, and value) transformation and spectral features (Wang et al., 2018). Most of these methods, however, only produce rapeseed maps for a small area using very limited imageries taken on rapeseed peak flowering dates (Ashourloo et al., 2019; She et al., 2015). Rapeseed peak flowering dates vary by area and cultivar because of differences in natural conditions and cultivation habits, especially over large regions (d'Andrimont et al., 2020; Ashourloo et al., 2019; McNairn et al., 2018). Using the above methods to automatically map rapeseed areas with a finer resolution over a large region is thus still a huge challenge.

Taking into consideration the unique phenological characteristics of different crops, various researchers have developed potentially useful phenology-based methods for crop identification over large areas (Ashourloo et al., 2019; Dong et al., 2016; J. Liu et al., 2018; Zhang et al., 2020). These algorithms, which generate classification rules by analyzing the unique characteristics of the studied crop, have been successfully applied for mapping rice (Dong et al., 2016), soybean (Zhong et al., 2014), corn (Zhong et al., 2016), and sugarcane (J. Wang et al., 2020) but have rarely been applied to rapeseed. Rapeseed has unique reflectance and scattering characteristics (Ashourloo et al., 2019; McNairn et al., 2018; Sulik and Long, 2015, 2016) and has three canopy morphologies based on leaves, yellow petals, and pods/branches (Ashourloo et al., 2019; Rondanini et al., 2014). Each canopy shape strongly influences how solar radiation is intercepted (Sulik and Long, 2016). The specific features of reflectance values and scattering coefficients of rapeseed from S-1/2 data can thus provide information for the automatic mapping of rapeseed over larger areas and with a finer resolution.

Another relevant aspect of rapeseed imaging concerns crop rotation, which is beneficial for pest and disease management in crop production (Harker et al., 2015; J. Liu et al., 2018) and a major factor in rapeseed yield (Harker et al., 2015; Ren et al., 2015). The physical and chemical properties of soil are altered during crop rotation, and these changes affect rapeseed growth (Ren et al., 2015). Most current studies have been limited to field observations (Peng et al., 2015). The spatial distribution of rapeseed rotation in different regions is still unclear because high-resolution rapeseed maps are lacking. To aid cultivation and management, the characteristics of rapeseed rotation need to be explored.

Taking into account the above-mentioned issues, we integrated multi-source data to develop a new method for identifying rapeseed. We then applied our new method to generate rapeseed maps with a spatial resolution of 10 m across the main planting areas of 33 countries from 2017 to 2019 and analyzed the geographical characteristics of rapeseed cultivation and crop rotation.

2 Materials and methods

2.1 Study area

We identified rapeseed planting areas in 33 countries, the world's main rapeseed producers, on three continents: North America, South America, and Europe (Fig. 1). The largest areas of rapeseed cultivation are located in Canada and the European Union (Carré and Pouzet, 2014; van Duren et al., 2015; Rondanini et al., 2012). In 2008, 79 % of biodiesel feedstock crops in Europe, which produces a large amount of biodiesel for export every year, were rapeseed (van Duren et al., 2015). In Chile, the main rapeseed producer in South America, the yield of rapeseed in 2018 was 3887.7 kg ha⁻¹. Rapeseed cultivation in these countries is important for food and energy security (Carré and Pouzet, 2014). The climates of the three continents are different because of factors such as latitude and topography (Peel et al., 2007). The rapeseed planting season varies among countries because of these differences in climate and other natural conditions (Singha et al., 2019; Wang et al., 2018), thus posing great challenges to the mapping of rapeseed.

Figure 1 Locations of 10 km radius sample blocks for phenological monitoring in the 33 countries in this study. The 33 countries include Canada (CAN), United States of America (USA), Chile (CHL), Ireland (IRL), United Kingdom of Great Britain and Northern Ireland (GBR), France (FRA), Spain (ESP), Netherlands (NLD), Belgium (BEL), Luxembourg (LUX), Germany (DEU), Switzerland (CHE), Denmark (DNK), Sweden (SWE), Poland (POL), Czechia (CZE), Austria (AUT), Slovenia (SVN), Croatia (HRV), Slovakia (SVK), Hungary (HUN), Estonia (EST), Latvia (LVA), Lithuania (LTU), Belarus (BLR), Ukraine (UKR), Republic of Moldova (MDA), Romania (ROU), Bulgaria (BGR), Serbia (SRB), North Macedonia (MKD), Greece (GRC), and Turkey (TUR). Country names and codes are the same as those used by the Statistics Division of the United Nations Secretariat. The three-digit alphabetical codes assigned by the International Organization for Standardization can be found at https://unstats.un.org/unsd/methodology/m49/ (last access: 8 June 2021).

2.2 Data

2.2.1 Remote sensing data

We collected imagery from Sentinel-2 (S2) and Sentinel-1 (S1) satellites (Table 1) launched by the European Space Agency (ESA) (Drusch et al., 2012; Torres et al., 2012). We used 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 poor-quality images (e.g., cloud-obscured information) in this study. Removal of all such information was difficult, however, because of the quality of the QA band (J. Wang et al., 2020; Zhu et al., 2015). We used the interferometric wide-swath mode of S1, which provides dual-band cross-polarization (VV) and vertical transmit–horizontal receive (VH) with a 12 or 6 d repeat cycle and 10 m spatial resolution (Torres et al., 2012). The S-1/2 images were obtained using the Google Earth Engine (GEE) (Gorelick et al., 2017). In addition, we used QA bands to remove most of the poor-quality images on GEE. (Sample code can be found at https://code.earthengine.google.com/?scriptPath=Examples:Datasets/COPERNICUS_S2, last access: 8 June 2021) Further details are provided in Table 1.

2.2.2 Digital elevation model

We used a spatial resolution of 1 arcsec (approximately 30 m) elevation data from the Space Shuttle Radar Terrain Mission (Table 1) (Farr et al., 2007). We then calculated the spatial distribution of slope using GEE (sample code at https://code.earthengine.google.com/?scriptPath=Examples:Datasets/USGS_SRTMGL1_003, last access: 8 June 2021). Finally, we extracted areas with a slope of less than 10^∘ to mask hilly terrain (Jarasiunas, 2016).

2.2.3 Cropland and agricultural statistical 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). Existing crop data products containing rapeseed information comprise four datasets: (1) the 30 m Annual Crop Inventory (ACI) in Canada (Fisette et al., 2013) and (2) the 30 m Cropland Data Layer (CDL) in the USA (Boryan et al., 2011), both of which were downloaded from GEE; (3) the Crop Map of England (CROME) generated in GBR; and (4) the 10 m Land Cover Map of France (LCMF) in France (Inglada et al., 2017). These four crop layer products were generated from satellite images and a large number of training sample collections. In this study, rapeseed maps in ACI, CDL, CROME, and LCMF were used for accuracy verification at the pixel level. For accuracy verification, we selected statistics on major crop areas in different countries and regions released annually by the Food and Agricultural Organization of the United Nations (FAO). Details are provided in Table 1.

2.2.4 Crop calendars

We used two crop phenological datasets to assist in the extraction of rapeseed phenological parameters: crop calendars in different countries (https://ipad.fas.usda.gov/ogamaps/cropcalendar.aspx, last access: 8 June 2021) and field records of crop phenology in Germany. The crop calendars originated from the United States Department of Agriculture, which only records rapeseed planting and harvest times for selected countries. The crop phenology field records in Germany were in situ observations from crop phenological records shared by the Deutscher Wetterdienst (DWD) in Germany (Kaspar et al., 2014). The DWD provides field observations of crop phenological periods throughout Germany following the Biologische Bundesanstalt, Bundessortenamt, and Chemical Industry (BBCH) scale (Table 1). DWD records include the start and end dates of rapeseed flowering (d'Andrimont et al., 2020; Kaspar et al., 2014). Neither the two crop calendars nor the DWD records contain information on rapeseed peak flowering dates. To extract rapeseed phenological parameters, we used all stations that fully recorded start and end dates of the flowering period from 2017 to 2019, namely, 281, 269, and 253 stations in 2017, 2018, and 2019, respectively.

Table 1 Detailed information about the data collected in this work.

2.3 Methods

2.3.1 Optical and SAR characteristics during the rapeseed growing period

We selected available rapeseed parcels and in situ observations of the DWD from different climate regions and different years to analyze the optical (reflectance and vegetation index) and SAR (VV and VH) characteristics of rapeseed over time. As an example, Fig. 2 shows the time series of one rapeseed parcel around a DWD station (station ID: 13126) in 2018. This rapeseed parcel exhibited unique visual characteristics during the flowering period (Fig. 2e). When rapeseed approached peak flowering, the flowers became yellow (d'Andrimont et al., 2020; Pan et al., 2013; Tao et al., 2020; Wang et al., 2018). Rapeseed was yellow-green on the true-color images of S2 and Google Earth during the flowering period (Fig. S1). The reflectance of the green and red bands separately increased – from 0.09 and 0.06, respectively, before flowering (17 April 2018) to 0.16 and 0.14 at peak flowering (7 May 2018) – and then decreased (Fig. 2a). The reflectance of the blue band was lower than that of the red and green bands during flowering. This outcome is similar to the results of previous research (Ashourloo et al., 2019; Sulik and Long, 2015). We also calculated the normalized difference yellow index (NDYI), which can capture increasing yellowness in a time series (d'Andrimont et al., 2020; Sulik and Long, 2016), as follows:

$$\begin{array}{}\text{(1)}& \mathrm{NDYI}={\displaystyle \frac{\mathrm{green}-\mathrm{blue}}{\mathrm{green}+\mathrm{blue}}},\end{array}$$

where green is the TOA reflectance of the green band (b3) of the S2 imagery, and blue is the blue band (b2) reflectance. NDYI increased from −0.03 on 17 April to 0.21 on 7 May (Fig. 2b) and reached a peak during rapeseed flowering. This unique spectral feature of rapeseed during the flowering period was due to the yellow petals.

Figure 2 Time-series profiles of four features of rapeseed pixels around one DWD station (ID = 13126; longitude: 11.333268424^∘ E, latitude: 52.200000463^∘ N) in Germany in 2018. (a) Mean reflectance values (red, green, and blue). (b) Mean NDYI. (c) Mean VH. (d) Mean VV. The light-shaded areas indicate the standard deviation. The BBCH scale was used for in situ observations of rapeseed phenology, with BBCH61 and BBCH69 respectively corresponding to the start and end of flowering. (e) The rapeseed parcel around the DWD station is bounded in red (image source: Copernicus Sentinel-2 data 2018).

S1 backscattering changes with rapeseed growth. We used VV and VH time series smoothed by the Savitzky–Golay (SG) filter (window size 3) (Chen et al., 2004) as inputs to identify the phenological parameters of rapeseed parcels. We used the SG filter algorithm in MATLAB 2020b, which uncovered local minima in both the VV (−11.20, 8 May) and VH (−15.60, 5 May) time series during rapeseed flowering (Fig. 2c, d). Furthermore, VH reached a maximum (−9.64, 1 June) during the pod period. Unlike other crops, rapeseed has two distinct green-up phases: the flowering period and the pod period (Ashourloo et al., 2019; Bargiel, 2017; Mercier et al., 2020; Veloso et al., 2017). The petals of rapeseed decrease the scattering of VV and VH, while the pods increase the scattering intensity of VH (d'Andrimont et al., 2020; Bargiel, 2017; McNairn et al., 2009; Mercier et al., 2020). The NDYI and backscattering (VV and VH) time series of rapeseed in different climate regions (Fig. S1) exhibited the same characteristics. Therefore, we used both optical and SAR features to identify rapeseed flowering and pod periods in this study. Because of differences in the revisit times of the S1/2 satellites, rapeseed peak flowering dates were defined as the median dates extracted using optical and SAR indicators.

2.3.2 Sample blocks collected for phenological monitoring

As a prerequisite to large-scale mapping, the phenology of rapeseed in different countries must be identified and delineated (Dong et al., 2016; Zhang et al., 2020), but not enough observational records of rapeseed phenology are available on a large scale. In accordance with the DWD method of phenological observation (Kaspar et al., 2014), we created sample blocks with a radius of 10 km over rapeseed-producing areas of different countries and randomly sampled 10 rapeseed parcels per block. The rapeseed plots were identified by their phenological characteristics, which were obtained by visual interpretation and analysis of reference data, including high-resolution images available in S2 and Google Earth as well as spectral reflectance (red and green bands), spectral index (NDYI), and scattering coefficient profiles (VV and VH) from the S1/2 time series. Google Earth images taken during rapeseed flowering were used to assist with the visual interpretation of rapeseed parcels. Rapeseed parcels with no available high-quality time-series imagery were omitted. Finally, 75 sample blocks in 2017, 84 sample blocks in 2018, and 84 sample blocks in 2019 were uniformly and randomly collected (Fig. 1).

2.3.3 Detection of flowering and pod phases in different countries

To determine the flowering dates of rapeseed in different countries, we evaluated each phenological sample block from 2017 to 2019 (Fig. 3). First, we calculated the average values of all pixels in the 10 previously selected rapeseed parcels in each block during the rapeseed growth period in conjunction with the crop calendar. VV and VH time series for each sampled rapeseed parcel were smoothed using the SG filter. Second, S1/2 peak flowering dates and pod dates were derived for all sample blocks according to the method in Sect. 2.3.1. We found that the peak flowering dates of rapeseed, especially in Europe, followed an obvious latitudinal gradient (Fig. 3j).

Figure 3 The spatial distribution of rapeseed flowering dates. (a)–(i) Flowering dates (Julian day) in different sample blocks in 2017, 2018, and 2019. (j) Characteristics of the latitudinal gradient in Europe. The peak flowering date for each latitudinal interval is the mean of the flowering dates of all sample blocks within that interval.

We also observed that the signal with the maximum VH occurred within 45 d of the peak flowering date (Fig. S2). We then calculated the difference in the peak flowering date of each sample block between different years, which revealed that the flowering peak dates of most sample parcels were advanced or delayed by only 10 d (Fig. 4d). Using the same period for different years in a given area was thus considered to be reasonable for rapeseed identification in this study. Previous studies and field observation records have indicated that the flowering period of rapeseed is approximately 30 d (d'Andrimont et al., 2020; Chen et al., 2019; Kaspar et al., 2014; She et al., 2015). Therefore, we divided each month into two time periods, with the 15th day serving as the dividing line. Two consecutive half-months were defined as a suitable period for classifying flowering dates (Fig. 4a–c). Finally, we designated the flowering period for each administrative unit based on the sample blocks.

Figure 4 Flowering phenology of rapeseed. (a)–(c) The spatial distribution of rapeseed flowering periods among sample blocks. (d) Boxplot showing changes in peak flowering dates of sample blocks over 3 years.

Figure 5 “Rainbow” cloud effect origins, examples, and solutions. (a) Staggered detector configuration of S2 (ESA, 2015a). (b)–(f) Examples of spectral misregistration effects and the performance of cloud-masking methods. Each image was masked with the quality assurance band (QA60) for the Sentinel-2 TOA image. The red arrows indicate the “rainbow” appearing around clouds at high altitudes in the S2 image (image source: Copernicus Sentinel-2 data). (g) Sentinel-2 TOA image of rapeseed at the flowering stage. The yellow arrow indicates rapeseed fields. (h) Scatter plot of NDYI vs. NRGBI of rapeseed field samples and samples with a “rainbow” around clouds in the S2 images. Relative pixel density is indicated by the color scale on the right.

2.3.4 Development of a phenology- and pixel-based algorithm for mapping rapeseed

Our temporal profiling of rapeseed parcels along with the results of many previous studies indicated that the spectrum at the flowering stage and the scattering signal at the pod stage are key features for identifying rapeseed (Ashourloo et al., 2019; Bargiel, 2017; Han et al., 2020; Mercier et al., 2020; Sulik and Long, 2015; Veloso et al., 2017). We developed a single phenology- and pixel-based rapeseed mapping algorithm that relies on four features: spectral bands (red and green), spectral indices (NDYI), polarization bands (VH), and terrain (slope). Four primary steps were used to map annual planting areas (Fig. 6).

Figure 6 Workflow for mapping rapeseed areas using the proposed phenology- and pixel-based algorithm. GFSAD30, Global Food Security-Support Analysis Data at 30 m; NDYI, normalized difference yellowness index; NRGBI, the new spectral index; DWD, Deutscher Wetterdienst; FAO, Food and Agriculture Organization of the United Nations; RMSE, root-mean-square error; MAE, mean absolute error; R², R-squared; CDL, Cropland Data Layer; ACI, Annual Crop Inventory; CROME, Crop Map of England; LCMF, Land Cover Map of France; UA, user's accuracy; PA, producer's accuracy; F1, F1 score.

In the first step, we determined the threshold of the feature indicators. We analyzed the histograms of random samples selected from different countries as suggested by previous study (Zou et al., 2018). Green and blue bands and NDYI (Fig. S3) were similar during the flowering stage in all samples from the different regions. Most rapeseed pixels (98 %) had the following values: red > 0.07, green > 0.11, and NDYI > 0.05. We observed some pixels, however, with a relatively high NDYI due to contamination by a cloud with a “rainbow” appearance, which would cause them to be misclassified as rapeseed (Fig. 5). Because of the limited quality of the QA band and the simple cloud score algorithm, such misclassifications arising from poor-quality observations from the S2 image cannot be removed (X. Wang et al., 2020b; Zhu et al., 2015). The rainbow in the cloud is the result of the push-broom design of S2 (Fig. 5a) and spectral misregistration (for more details, see ESA, 2015a, and b). Taking into account the relative displacement of each spectral channel sensor in the S2 push-broom design (Frantz et al., 2018; X. Liu et al., 2020; Zhao et al., 2018), we developed a new spectral index (NRGBI) to reduce the influence of the rainbow (Eq. 2):

$$\begin{array}{}\text{(2)}& \mathrm{NRGBI}={\displaystyle \frac{\mathrm{red}-\mathrm{blue}}{\mathrm{red}+\mathrm{blue}}}-{\displaystyle \frac{\mathrm{green}-\mathrm{blue}}{\mathrm{green}+\mathrm{blue}}},\end{array}$$

where red, green, and blue are the TOA reflectance values of the red (b4), green (b3), and blue (b2) bands of the S2 imagery, respectively. A scatter plot of NDYI vs. NRGBI of rapeseed parcel samples and “rainbow” samples around clouds (visual interpretation) demonstrated that the NRGBI (threshold = −0.05) can effectively distinguish rapeseed from the rainbow (Fig. 5h). The GEE code for NRGBI index calculations can be found at https://code.earthengine.google.com/a39fc699a276d018778d59c5b085d960 (last access: 8 June 2021). In addition, NRGBI can be calculated using Eq. (2) in other GIS software programs (e.g., QGIS and ArcGIS) on a local computer.

The second step was the identification of all rapeseed pixels from different images during the flowering period and their subsequent aggregation into annual rapeseed planting areas (Fig. 6). Because peak flowering dates and the number of available images of rapeseed fields vary within a region (Fig. S4), rapeseed classifications based on a single image may fail to capture rapeseed flowering dynamics (Ashourloo et al., 2019). To avoid misclassification due to poor-quality observations during the rapeseed flowering stage, we aggregated all results classified from available S2 images during this period. The use of a larger number of images resulted in better performance (Fig. S4).

In the third step, we combined optical data with SAR images to ensure the accuracy of the rapeseed maps. High VH values during the pod stage are another distinct feature that can distinguish rapeseed from other crops (Mercier et al., 2020; Tian et al., 2019; Van Tricht et al., 2018; Veloso et al., 2017). Taking into consideration the variability of flowering in different fields and the duration of the pod stage (Sect. 2.3.2), we calculated the maximum VH between the second half of the flowering stage and the next 30 d (ca. 45 d; Fig. 6). Within this 45 d interval, at least three S1 satellite images were available in the study area. In addition, areas with a slope ≥10^∘ were removed (Jarasiunas, 2016). All pixels meeting these requirements were defined as rapeseed.

In the fourth step, we removed “salt and pepper” noise by applying a threshold based on the number of connected components (objects), that is, the size of the neighborhood in pixels, and then filling the gaps inside the parcels (Hirayama et al., 2019). We used an eight-connected rule, which means that the edges or corners of the pixels were connected. If two adjacent pixels were connected, they were considered part of the same object (https://www.mathworks.com/help/images/ref/bwareaopen.html, last access: 8 June 2021). The bwareaopen function in MATLAB (version 2020b) was used to remove objects not meeting a given threshold. The thresholds of different indicators in different regions are given in Table S1.

2.4 Accuracy assessment

To test the accuracy of our proposed algorithm, we first compared rapeseed areas retrieved using the new method with FAO statistics. Our rapeseed data constituted a binary (0 or 1) map with a spatial resolution of 10 m. We calculated the total area of rapeseed maps in each country and compared these numbers with FAO national rapeseed statistics. To verify the accuracy of rapeseed mapping, we used the RMSE (Eq. 3), MAE (Eq. 4), and coefficient of determination (R², Eq. 5), which were calculated as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathrm{RMSE}=\sqrt{{\sum }_{i=\mathrm{1}}^n{\displaystyle \frac{{\left(y_i-f_i\right)}^{\mathrm{2}}}{n}}},\text{(4)}& {\displaystyle }\mathrm{MAE}={\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n\left|y_i-f_i\right|,\text{(5)}& {\displaystyle }{R}^{\mathrm{2}}={\displaystyle \frac{{\left({\sum }_{i=\mathrm{1}}^n(y_i-\stackrel{\mathrm{‾}}{y_i})(f_i-\stackrel{\mathrm{‾}}{f_i})\right)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n{\left(y_i-\stackrel{\mathrm{‾}}{y_i}\right)}^{\mathrm{2}}{\sum }_{i=\mathrm{1}}^n{\left(f_i-\stackrel{\mathrm{‾}}{f_i}\right)}^{\mathrm{2}}}},\end{array}$$

where n is the total number of countries, y_i is the mapped rapeseed planting area for country i, $\stackrel{\mathrm{‾}}{y_i}$ is the corresponding mean value, f_i is the rapeseed planting area recorded by the FAO for country i, and $\stackrel{\mathrm{‾}}{f_i}$ is the corresponding mean value.

We also compared our rapeseed maps with four open-access datasets that include rapeseed layers at the pixel level: ACI, CDL, CROME, and LCMF in Canada, the USA, GBR, and France, respectively. We used the data from 2018 and 2019 in these datasets as a reference (Boryan et al., 2011; Fisette et al., 2013). To unify the spatial resolution of our rapeseed maps, we resampled CDL, ACI, and CROME to 10 m resolution to allow comparison. To check the accuracy of our classification, we calculated UA (Eq. 6), PA (Eq. 7), and F1 (Eq. 8) based on confusion matrices (Table S2).

We also randomly selected verification samples based on previous studies (Pekel et al., 2016; X. Wang et al., 2020a) to validate our rapeseed maps. A 0.2^∘ × 0.2^∘ latitude–longitude grid was superimposed on our 2018 rapeseed map (Fig. S5). Two points – one rapeseed and the other non-rapeseed – were randomly generated in each grid by visually interpreting images available from S2 and Google Earth along with spectral reflectance (red and green bands), spectral index (NDYI), and scattering coefficient (VV and VH) profiles from the S1/2 time series. Confusion matrices were similarly used to assess accuracy according to Eqs. (6)–(8):

$$\begin{array}{}\text{(6)}& {\displaystyle }\mathrm{UA}={\displaystyle \frac{x_{ij}}{x_j}},\text{(7)}& {\displaystyle }\mathrm{PA}={\displaystyle \frac{x_{ij}}{x_i}},\text{(8)}& {\displaystyle }\mathrm{F}\mathrm{1}=\mathrm{2}\times {\displaystyle \frac{\mathrm{UA}\times \mathrm{PA}}{\mathrm{UA}+\mathrm{PA}}}.\end{array}$$

In the above equations, x_ij is the value of the ith row and jth column, x_i is the sum of the ith row, and x_j is the sum of the jth column. Although the statistical data and existing products did not completely reflect the actual areas and locations of cultivated rapeseed, these datasets were still beneficial for validating the accuracy of our rapeseed maps at national and pixel scales.

3 Results

3.1 Accuracy assessment

We compared our derived rapeseed areas with those from FAO statistics. The total planting areas of rapeseed exhibited good consistency with the agricultural statistics at the national level, with a RMSE of 1459.64 km², a MAE of 785.25 km², and an R² of 0.88 (Fig. 7). We found that the derived areas in 2018 were larger than those in 2017 and 2019, especially in countries with relatively small rapeseed areas. The greater availability of S2 images and higher-quality data in 2018 may have contributed to the derivation of the larger areas by our new method (X. Liu et al., 2020).

Figure 7 Comparison of rapeseed areas with FAO statistics at the national level.

Figure 8 Classification validation results. (a) Percentage of rapeseed areas based on FAO statistics classified as such in existing products and our rapeseed map database. (b) Accuracy of our classifications in four countries (Canada, USA, GBR, and France) using existing products as a reference. UA, user's accuracy; PA, producer's accuracy; F1, F1 score.

Figure 9 Spatially explicit details of rapeseed maps in eight countries with diverse crop structures in different years. The names of climate zones are given in yellow. RGB composite images comprise red (b4), green (b3), and blue (b2) bands from Sentinel-2 good-quality observations during the rapeseed flowering period (image source: Copernicus Sentinel-2 data). Climate zone data are from the Food Insecurity, Poverty and Environment Global GIS Database.

As indicated by their accuracy based on confusion matrix values, our rapeseed maps were consistent at the pixel level with maps of the American CDL in 2018 and 2019 and the Canadian ACI, British CROME, and French LCMF in 2018 (Table S3). As shown in Fig. 8a, the rapeseed areas calculated from our maps were consistently more comparable to FAO statistics than those from existing products. UA, PA, and F1, which varied by country, ranged from 0.93–0.97, 0.70–0.80, and 0.81–0.86, respectively (Fig. 8b). The rapeseed areas determined using our algorithm accounted for approximately 71 % of the 2018 CDL, 71 % of the 2018 ACI, 80 % of the 2018 CROME, 70 % of the 2018 LCMF, and 79 % of the 2019 CDL. In addition, the distributions on our rapeseed maps were consistent with those of existing products at the pixel level (Figs. S7 and S8). The differences in accuracy may have been due to the varied number of high-quality images available in different regions (Dong et al., 2016). Despite the different ground conditions, methods, images, and spatial resolutions among products, the results of our comparison further verify the accuracy of our rapeseed maps (Gong et al., 2020; Singha et al., 2019).

According to confusion matrix values (Table S4) based on random sampling points, the accuracy of our rapeseed maps varied in different regions. We obtained the highest accuracy (F1, 0.91) in zone II, followed by zone III (F1, 0.9), and zone I (F1, 0.84). These disparities may be due to differences in the availability of S1/2 images among the studied areas.

3.2 Additional features of rapeseed maps derived using our new method

To further characterize the rapeseed maps generated in this study, we selected various images in several areas of each country. The rapeseed maps showed good spatial consistency with the actual areas of rapeseed cultivation on the ground (Figs. 9 and S6). Fields with various planting densities, ranging from densely planted areas in Canada (Fig. 9a) to relatively sparse ones, such as in Chile (Fig. 9b) and European countries (e.g., Fig. 9c, d) (Lowder et al., 2016); various shapes, ranging from regular rectangles (e.g., Fig. 9a, h) to irregular parcels (Fig. 9c, d); and different climatic conditions, ranging from a temperate oceanic climate (Fig. 9c–e) to temperate sub-continental (Fig. 9a, f) or even subtropical (Fig. 9b) ones, were clearly and comprehensively indicated on our maps. The fragmented pattern of land in some European countries, especially that in eastern and central Europe due to land reform in 1989 (Hartvigsen, 2013, 2014), was clearly evident; Fig. 9f shows land in Estonia as an example (Jürgenson and Rasva, 2020; Looga et al., 2018). Although the algorithm was applied to different climates, terrains, and landscapes over a very large region, its classification accuracy across 33 countries was satisfactory. Our rapeseed maps can thus effectively identify fields in detail with high spatial resolution and clear field boundaries. More rapeseed classification details can be found in Fig. S6.

Figure 10 Spatial distribution of rapeseed areas at 10 m resolution along latitudinal and longitudinal gradients in 2018. (a) Europe and Turkey. (b) Canada and the USA.

3.3 Spatial patterns of rapeseed planting areas

In our maps, the largest total area of rapeseed cultivation worldwide was in Canada. Along with GBR, Poland, and Ukraine, the two leading rapeseed growing countries in Europe – France and Germany – accounted for approximately 66.3 % of European rapeseed areas. The 3-year (2017–2019) spatial patterns were consistent at the national level (Fig. S9). We also plotted the geographic distribution of rapeseed areas along latitudinal and longitudinal gradients in the study areas (Fig. 10). With the exception of steep mountainous regions and cold northern areas, rapeseed is widely planted in European countries at latitudes of 46–53^∘ N and longitudes of 2^∘ W–4^∘ E, 9–19^∘ E, and 22–27^∘ E (Fig. 10a). In Canada and the USA, areas with latitudes of 49–54 and 56–57^∘ N and longitudes of 117 to 118^∘ W and 98 to 114^∘ W have high densities of planted rapeseed (Fig. 10b).

4 Discussion

4.1 Investigation of rapeseed rotation systems

We obtained 3-year rapeseed maps at a spatial resolution of 10 m whose high accuracy was validated by annual national statistics books, open-access public products, and random sampling points. These rapeseed maps provided a new opportunity to investigate rapeseed rotation systems (X. Liu et al., 2018). Crop rotation information is an important factor in crop yield management (Harker et al., 2015; X. Liu et al., 2018; Ren et al., 2015; Rudiyanto et al., 2019; Zhou et al., 2015). To analyze rapeseed rotation patterns, we selected 25 representative areas (Fig. S10) that met the following three criteria: high image quality, high rapeseed classification accuracy, and large extent of planted rapeseed. Rapeseed rotation in these areas was calculated based on the frequency of each rapeseed pixel (Fig. 11).

Figure 11 Spatial distribution of three rotation schedules in different areas from 2017–2019.

Because only 3 years (2017–2019) of rapeseed maps were available, the longest observable rapeseed rotation break was 2 years. To more accurately discern the pattern of rapeseed rotation, we thus classified rapeseed rotation breaks in this study into three types: ≥2 years, 1 year, and 0 years. Most countries, especially European ones, were characterized by rotation breaks that were ≥2 years (mostly green areas in Fig. 12). In Canada, 70 % of fields were subjected to rotation breaks of ≥2 years, with the remainder (30 %) following a 1-year break pattern (Fig. 12a). As shown in the histogram in Fig. 12d, we identified 20 locations with ≥2-year rotation breaks, which corresponds to 90 % of planting areas. Many previous studies have found that a 2- or 3-year rotation break significantly reduces the number of fungal spores, especially those of Rhizoctonia solani and Leptosphaeria maculans, thus suggesting that a rotation system is an important component of disease control in rapeseed (Gill, 2018; Harker et al., 2015; Ren et al., 2015; Zhou et al., 2015). Rapeseed rotation also improves yield, moisture, and fertility and reduces weeds and pest insects (Bernard et al., 2012; Harker et al., 2015; Pardo et al., 2015; Peng et al., 2015; Ren et al., 2015). Additional efforts to produce longer time-series rapeseed maps and acquire detailed rotation information are thus needed.

Figure 12 Rapeseed crop rotation. (a)–(c) Proportions of rapeseed planting areas subjected to rotation breaks of 0, 1, or ≥2 years. (d) The number of areas in (a)–(c) subjected to breaks of at least 2 years.

Figure 13 Example showing the effect of low-quality observations on classification accuracy. (a) Rapeseed map of an area of France in 2017 that contains an error (longitude: 2.059824^∘ E, latitude: 46.734987^∘ N). (b) Availability of time-series Sentinel-2 images during rapeseed flowering phases. (c) Comparison of the time series of different sites showing how the peak NDYI has been missed.

4.2 Uncertainty

The generation of annual high-resolution maps for a specific crop over a large region is a major challenge (Dong et al., 2016; Luo et al., 2020; Zhang et al., 2020). Pixel-and phenological-based algorithms, multisource remote sensing data, and the GEE are useful for mapping rapeseed at high resolution and over large areas. In addition to these advantages, our proposed algorithm does not require large amounts of training sample data and reduces disturbance due to agronomic differences by combining images from multiple dates. Nevertheless, uncertainty still exists due to several issues. The first of these factors is the cropland layer. We used GFSAD30 datasets to identify cropland; however, the GFSAD30 has limitations, such as classification errors (Phalke et al., 2020). A second contributory aspect is the number of satellite images available. Although our annual rapeseed maps are consistent with FAO statistics and show higher accuracy compared with existing products, they are limited by the quality of the observations during critical growth stages. For example, Fig. 13a shows an error in an area of France in 2017 that can be attributed to the lack of clear S2 images during the rapeseed flowering period (Fig. 13b). Because the rapeseed flowering period is generally characterized by high NDYI and high red and green band reflectance, rapeseed pixels are likely to be misclassified if images are missing during the flowering stage (Fig. 13c). A third issue concerns the threshold for different indicators, which is a key factor for mapping crops (Ashourloo et al., 2019; Dong et al., 2016; X. Liu et al., 2020; J. Wang et al., 2020; Zhang et al., 2015). Although reference thresholds for the three continents in this study are provided, they should be applied with caution to other regions. Finally, the complexity of the ground environment can contribute to uncertainty. For example, landscape types might impact the accuracy of rapeseed maps (J. Wang et al., 2020).

5 Data availability

The rapeseed maps produced with 10 m resolution in this study are accessible at Mendeley Data (http://dx.doi.org/10.17632/ydf3m7pd4j.3) (Han et al., 2021). 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 independently verify the rapeseed maps. In addition, Sentinel 1/2 images, CDL, ACI, and SRTM are available on GEE (https://developers.google.com/earth-engine/datasets/, Google Earth Engine Data Catalog, 2021). For more detailed information about the data collected in this work, see Table 1.

6 Conclusions

Large-scale, high-resolution rapeseed maps are the basis for crop growth monitoring and yield prediction. We developed a new method for mapping rapeseed based on spectral and polarization features and multi-source data. We used the new algorithm to produce three annual rapeseed maps (2017–2019) at 10 m spatial resolution in 33 countries. According to the results of three different verification methods, our rapeseed maps are reasonably accurate. Compared with existing products at the pixel level in Canada, USA, GBR, and France, PA, UA, and F1 were 0.70–0.80, 0.93–0.97, and 0.81–0.86, respectively. In addition, F1 ranged from 0.84 to 0.91 based on independent validation samples. Our approach reduces misclassifications due to different planting times and low-quality observations to some degree. The 10 m rapeseed maps provide more spatial details of rapeseed. Finally, we observed that the rapeseed crop rotation interval is at least 2 years in almost all countries in this study. Our proposed rapeseed mapping method can be applied to other regions. The derived rapeseed data product is useful for many purposes, including crop growth monitoring and production and rotation system planning.