Mapping rapeseed planting areas using an automatic phenology-and pixel-based algorithm (APPA) in Google Earth Engine

Jichong Han,Zhao Zhang,Juan Cao,Yuchuan Luo

Academy of Disaster Reduction and Emergency Management Ministry of Emergency Management & Ministry of Education,School of National Safety and Emergency Management,Beijing Normal University,Beijing 100875,China

Keywords:Automatic mapping Spectral indices Polarization Phenology Rapeseed

ABSTRACT The timely and rapid mapping of rapeseed planting areas is desirable for national food security.Most current rapeseed mapping methods depend strongly on images with good observations obtained during the flowering stages.Although vegetation indices have been proposed to identify the rapeseed flowering stage in some areas,automatically mapping rapeseed planting areas in large regions is still challenging.We developed an automatic phenology-and pixel-based algorithm (APPA) by integrating Landsat 8 and Sentinel-1 satellite data.We found that the Normalized Rapeseed Flowering Index shows unique spectral characteristics during the flowering and post-flowering periods,which distinguish rapeseed parcels from other land-use types(urban,water,forest,grass,maize,wheat,barley,and soybean).To verify the robustness of APPA,we applied APPA to seven areas in five rapeseed-producing countries with flowering images unavailable.The rapeseed maps by APPA showed consistently high accuracies with producer accuracies of (0.87-0.93 and F-scores of 0.92-0.95 based on 4503 verification samples.They showed high spatial consistency at the pixel level with the land cover Scientific Expertise Centres (SEC) map in France,Crop Map of England in United Kingdom,national-scale crop-and land-cover map of Germany,and Annual Crop Inventory in Canada at the pixel level.We propose APPA as a highly promising method for automatically and efficiently mapping rapeseed areas.

1.Introduction

Obtaining spatial and temporal information on the growth status and yield of crops is the first aim of local and national food security[1-6,7].Crop maps are essential for various studies in agriculture,including analyses of yield predictions [6,8],phenology inversions[9,10],and spatiotemporal cropland patterns [11].However,the mask layers of crops show lower accuracies in some global and regional yield prediction systems or Earth observation systems owing to their coarse spatial resolutions [7] and lack of available annual crop maps.The unavailability of high-quality (e.g.,timely and accurate)crop maps hinders agricultural research [12].

Recently,new remote-sensing technology and platforms,among them Google Earth Engine (GEE),have provided good opportunities for mapping crops [13].Diverse classification methods have been proposed,including machine learning methods[14-16],object-oriented methods [17],phenology-based methods[18-20],and combinations of these [21].Several long-term time series of farmland mask datasets are popularly used,including the Cropland Data Layer in the U.S.and the Annual Crop Inventory in Canada from Agriculture and Agri-Food Canada [22,23].These products were generated by supervised classification methods,which depend on field surveys and large numbers of real samples[24].In developing countries,crop mask datasets are rare.

Various satellite data are used to map crops including optical and microwave measurements.Optical data provide spectral information about the reflective and emissive features of surfaces.Microwave data provide structural information about surface roughness and texture.Landsat and Sentinel-1/2 are among the most widely used optical data for crop mapping because they both have high spatial resolution and are publicly available [12,20,25].Both Landsat and Sentinel-2 optical images are affected by clouds.Sentinel-2 has higher spatial and temporal resolution (10 m) than Landsat(30 m).There are large differences between some spectral indices calculated from calibrated top-of-atmosphere(TOA)reflectance data and atmosphere-corrected surface reflectance(SR)data.Theoretically,SR data can better reflect the true spectral characteristics of crops than TOA data.Level-2A is the SR product of Sentinel-2,and global coverage of Level-2A started in December 2018 (https://sentinel.esa.int/web/sentinel/user-guides/sentinel-2-msi/product-types/level-2a,last access: September 22,2021).Thus,Landsat stores more available SR images of historical years than Sentinel-2 on a global scale.The Sentinel-1 satellite with 10-m spatial resolution is not affected by clouds.However,owing to the complexity of the land surface,it is difficult to identify a crop using using Sentinel-1 data.Sentinel-1 and Landsat data contain different but complementary spectral and spatial information[26-28].

Crop maps of rapeseed,a vegetable oil crop that accounts for much of the world’s vegetable oil supply,are needed.Previous studies that focused on mapping rapeseed can be categorized into five groups.(1) Field surveys or visual interpretations of remote sensing images during the flowering phase [3,4,29,30].This method is strongly dependent on the quality of the images obtained during the peak flowering dates.Consequently,it is difficult to apply the method owing to its time-consuming and labor-intensive properties,and using it in areas frequently disturbed by clouds and rain is also difficult,especially in tropical regions [31].(2) Machine learning methods[2,24,32].In these methods,classifications are mainly based on many training samples,which are difficult to obtain [16,24,32].Classification results vary among samples selected from different rapeseed stages or areas.(3) A threshold-segmentation algorithm based on continuous time-series data(such as Normalized Difference Vegetation Index (NDVI)and Enhanced Vegetation Index (EVI)) [20].Ashourloo et al.[12] identified rapeseed using the feature of increased reflectance in the red and green bands during the flowering stage.However,this method is applicable only if a long-term,good observations time series of images is available.(4) Extraction algorithms based on colorimetric transformations and spectral features.Wang et al.[4] proposed a method for mapping rapeseed in which the spectral and color information of images obtained during the flowering stages were combined.The change in color of the rapeseed during the flowering stage causes visual differences.They accordingly used a color space technology to extract rapeseed.(5) Extraction algorithms based on the incorporation of hyperspectral remote sensing data[29,33]by multirange spectral feature fitting.However,the remote-sensing images used in the fourth and fifth methods must be obtained during the rapeseed flowering stages,and the latter method is limited to small coverage areas and is unsuitable for largescale mapping owing to the higher costs of hyperspectral images.

All five groups of methods depend strongly on either images with good observations or images obtained during flowering stages or require many training samples.The flowering stages of rapeseed usually last less than one month.There are no strong color features at the non-flowering stage,making it difficult for the method to accurately identify rapeseed.For example,the long revisit cycle of Landsat and the fact that the images are frequently disturbed by clouds make it difficult to obtain good observations in some areas during the rapeseed flowering stages [12,24,32,34].A new method to map rapeseed automatically that does not rely on peak flowering images is needed.

Our objectives were (a) to identify the spectrum and polarization characteristics of rapeseed,(b) to develop an automatic phenology-and pixel-based algorithm (APPA),(c) to apply the APPA to seven regions in diverse geographical environments to test its robustness,and (d) to compare the performance of the APPA with existing data products.

2.Materials and methods

2.1.Study area

We selected seven sites in five rapeseed-producing countries:China,France,United Kingdom,Germany,and Canada (Fig.S1).The selected study areas cover diverse climate zones,topographic conditions,and crop calendars.In most study areas,there were not enough Landsat images available during the rapeseed flowering period owing to cloud interference (Fig.S2).

China (Qinghai) zone I is located in the northeastern Tibetan Plateau (101°17′48′′E-101°57′16′′E,37°15′4′′N(xiāo)-37°39′8′′N(xiāo)).The lowest temperatures in Qinghai occur in January and December,and the highest temperatures in July and August.Spring rapeseed is the main crop in this study area.Rapeseed is usually sown in May and harvested in September,and the flowering stage occurs during July,generally lasting for one month [35].

China (Hebei) zone II is situated in the southern North China Plain (114°42′27′′E-114°47′27′′E,36°8′54′′N(xiāo)-36°12′12′′N(xiāo)).Rapeseed grown in zone II is the main source of oil for local farmers.The climate is arid,steppe-like,and cold [36].Winter wheat and winter rapeseed are the main crops in this region.Rapeseed is usually sown in November and harvested in the following June,and the flowering stage is from mid-March to mid-April.

China (Hubei) zone III is located in the Yangtze River Basin where the climate is temperate,without a dry season,and the region experiences hot summers (115°41′7′′E-116°7′39′′E,29°41′57′′N(xiāo)-29°59′24′′N(xiāo)) [36].The Yangtze River Basin is the largest major rapeseed plantation area in China.Winter rapeseed is usually planted in October and harvested in May.The flowering stage is in March [35].

Zones IV and V are located in northern France (0°48′19′′E-2°2 9′49′′E,47°51′55′′N(xiāo)-49°30′29′′N(xiāo)) and the south of England in United Kingdom (0°4′23′′W-2°51′48′′W,50°22′19′′N(xiāo)-52°7′31′′N(xiāo)),respectively.Rapeseed is a main oil crop in both zones.The climate in these two regions is temperate,with no dry seasons and warm summers [36].Winter rapeseed is usually planted in September and harvested in June,and the flowering stage spans from April to May [10,27].

Zone VI covers Brandenburg and Sachsen in Germany(12°38′23′′E-14°17′15′′E,50°44′39′′N(xiāo)-52°5′1′′N(xiāo)) and has a climate characterized by cold temperatures,without a dry season and with warm summers [36].Winter rapeseed,maize,and other winter cereals are the main crops in this study area.Rapeseed is an economic oil crop in Brandenburg and Sachsen.Winter rapeseed is usually planted in September and harvested in July.The flowering stage is in May [37].

Canada zone VII is located in Saskatchewan,Canada(107°11′37′′W-109°7′29′′W,51°53′43′′N(xiāo)-53°14′26′′N(xiāo)),where the climate characteristics include cold temperatures,no dry seasons,and warm summers [36].Saskatchewan is one of the largest rapeseed-growing regions in Canada.Winter rapeseed,barley,corn,soybean,and winter wheat are the main crops in this region.Winter rapeseed is usually planted in May and harvested in September,and the flowering stage is in July [38].See Fig.S3 for more details.

2.2.Data

2.2.1.Satellite data

(1) Landsat 8 Surface Reflectance data and processing

All Landsat 8 OLI Surface Reflectance (SR) images used in our study were obtained from GEE (Tables S1,S2),a platform with high-performance computing capability [13].The reason we used Landsat 8 SR data instead of Sentinel-2′s Level-2A SR data is that the global coverage of Level-2A SR started in December 2018.In addition,some zones in our work were studied before 2018.Landsat 8 stores more available SR images globally for historical years than Sentinel-2.We used the Landsat Quality Assessment (QA)band instead of the blue band to remove clouds and cloud shadow pixels from all SR images[39].In comparison with the QA band,it is difficult to remove the shadow area in the blue band (Fig.S4).Only images that were unaffected by clouds or shadows were used for rapeseed mapping.We analyzed the spectral information of the green,red,near-infrared (NIR),and swir2 bands,with a 30-m resolution.We also used Sentinel-2 satellite images (at 10-m resolution) to collect reference samples and estimate accuracy [40](Fig.S5).

(2) Sentinel-1 backscattering coefficient data and processing

We also collected Sentinel-1 (S1) C-band synthetic aperture radar(SAR)imagery(from the Level-1 Ground Range Detected product) [41].The interferometric wide swath(IW) of the S1 imagery includes multiple polarizations,including dual-band crosspolarization (VV) and vertical transmit/horizontal receive (VH)polarization [41].The S1 imagery in the GEE data pool was preprocessed.More details can be found at https://developers.google.com/earth-engine/guides/sentinel1 (last access: June 22,2021).We removed speckle noise with a morphological mean filter (a circle kernel type).See https://developers.google.com/earth-engine/guides/image_morph (last access: June 22,2021) for more information about morphological operations.The spatial resolution of the VH time-series data from GEE was unified to 30 m in this study.

2.2.2.Digital elevation model

The 30-m-resolution elevation data from the Shuttle Radar Topography Mission (SRTM) were used to generate slope maps in the study areas[42].Areas with steep slopes(≥30°)that were unlikely to grow the crop were then masked out on the GEE platform.

2.2.3.Cropland data

We collected existing rapeseed map products to verify the accuracies of the rapeseed maps generated by our method for different regions (zones IV-VII): (1) land cover SEC data collected in France with a 10-m spatial resolution.The overall accuracy of the crop in SEC is about 0.90 [43];(2) a Crop Map of England (CROME) in United Kingdom with a hexagonal-cell(0.41-ha)spatial resolution.Both the producer accuracy (PA) and user accuracy (UA) of rapeseed in CROME exceed 0.85;(3) a national-scale crop-and landcover map of Germany (NCLM) with a 30-m spatial resolution[24].The PA and UA of rapeseed in NLCM are 0.87 and 0.97,respectively[24];and(4)and the Annual Crop Inventory(ACI)in Canada with 30-m spatial resolution.The cropland in ACI meets the overall target accuracy of at least 85% (https://developers.google.com/earth-engine/datasets/catalog/AAFC_ACI#description).All datasets can be found in Table S1.To unify the spatial resolutions of the rapeseed maps,the products were resampled to 30-m resolution.Note that the products used in this study are not taken as ground truth.Instead,the comparisons offer an overall assessment of agreement by treating these products as baselines in the absence of field sample points.

2.2.4.Field surveys

We collected field survey data in two counties (Wei and Linzhang County)in Handan city,Hebei province,China(zone II),during the rapeseed growing period in 2020.We continued to observe the rapeseed sample parcels from February to June with a revisit interval of 2-5 days.Finally,we collected time series pictures of rapeseed growth and 414 rapeseed sample points.The observation dates and geographic coordinate information were recorded (e.g.,Figs.1A,S5).We applied the pixels extracted from rapeseed parcels for threshold determination and statistical analysis.The field survey data were used to identify the spectral (reflectance and vegetation index) and polarization characteristics (backscatter coefficients) of rapeseed parcels and to test the accuracy of the rapeseed maps generated by the proposed method in zone II.

Fig.1.Field photos and index for rapeseed in China in 2020.(A)Field photos of growth stages of rapeseed in zone II(114°45′48′′E,36°10′22′′N(xiāo)).Differences in saturation of the photos are caused by the time and angle of observation.(B)Time-series profiles of the green,red,NIR,and swir2 bands for rapeseed.Contours indicate the mean values of the rapeseed reference samples obtained from the field investigations.Dark gray rectangles indicate the flowering stages.Llight gray rectangles indicate the month after the flowering stages (pod-filling stages).

2.2.5.Collecting reference samples

Using images with higher resolutions is one of the most extensive methods used to select reference samples [44-47].We collected reference samples by integrating Sentinel-2 TOA images of flowering rapeseed in other zones.The reason we used Sentinel-2 TOA data instead of Sentinel-2′s Level-2A SR data is that there are no Sentinel-2 Level-2A SR images available during the study periods of zones I and III.The visual difference caused by the change in the color of rapeseed during the flowering stages is a feature used to identify rapeseed [3,4,29,30,34,38,48].The rapeseed parcels in different waveband combinations of Sentinel-2 TOA imagery all show unique visual characteristics(Fig.S5).For example,rapeseed is yellow-green in the truecolor (R: G: B=band4: band3: band 2) images during the flowering stages (Fig.S5).Pure pixel reference samples were collected from the centers of the parcels to avoid mixed-edge pixels [49].Correspondingly,the reference samples were used to identify the differences between rapeseed and other land types (urban,water,forest,grass,maize,wheat,barley,and soybean).The pixels extracted from reference samples were used for determining thresholds.

2.3.Calculation of vegetation index

Yellow flowers are one of the most specific characteristics of rapeseed at peak greenness;the inflorescences of these yellow flowers appear in clusters at the ends of the stems [12,29,48].In comparison with other land cover types,this obvious visual difference in rapeseed can generally be observed during the whole flowering stage.Correspondingly,rapeseed also has unique spectral phenological characteristics [3,50].The green,red,and NIR bands showed continuous increases during the flowering stage (Fig.S6);these increases differed from the properties of other land cover types [3,12,29,50].Moreover,the reflectance values of swir2 remained stable throughout the whole flowering stage and the following month (pod-filling stage),with lower values than those of the green band (Fig.1B).We calculated the normalized rapeseed flowering index (NRFI) [27].Because the NRFI includes green and swir2 bands (Eq.(1)),it has a higher potential than other indices for differentiating rapeseed from other land cover types.The NRFI values were greater than 0 during the rapeseed flowering and the following month (pod-filling stage) (Figs.1D,S7).NRFI can be calculated as follows:

where ρgreenand ρswir2are respectively the green and shortwave infrared-2 bands of the Landsat 8 OLI imagery.

2.4.Feature development

The spectral and polarization characteristics provided valuable information for mapping rapeseed.We characterized the growth dynamics of NRFI and VH for the main land-use types,including urban,water,forest,grass,maize,winter wheat,barley,and soybean(Fig.S8).The reference samples were collected from field surveys and available land cover dataset(Table S1).We found that the NRFI values of all land use types except water bodies were less than 0.05.The NRFI of rapeseed (in the flowering and following month) and water were both greater than 0.The pod stage refers to the pod-filling stage in this study.However,the backscattering coefficient of water bodies was lower than that of rapeseed.Considering that the NIR of rapeseed reached maximum values (NIRmax) during the flowering and pod stages,NIRmax was used as an indicator to identify rapeseed.VH reached a maximum at the pod stage owing to the intertwined pods and stems of rapeseed increasing the backscattering coefficient.Accordingly,the maximum value of VH (VHmax) during rapeseed growth is considered a valid indicator to distinguish rapeseed from other crops [25].Given that rapeseed has a large range of VHs before and after flowering,we calculated the coefficient of variation of VH (VHcv) over the period from two months before flowering to two months after flowering.Thus,four indicators(NRFI,NIRmax,VHmax,and VHcv)were combined to identify rapeseed.

2.5.Generation of rapeseed map

Owing to variation in the timing of flowering and the pod stage among different regions,the rapeseed classification results obtained from one image or image composite cannot capture this dynamic [12].Additionally,some regions may not have good observations during the classification periods.To reduce the effects of phenology and the lack of images on mapping rapeseed,we used a combined approach from time series Landsat images to map rapeseed.The analysis processes of the APPA are shown in Fig.2B.The rapeseed classifier rules were as follows.

Step 1:The NRFI is calculated based on the red and green bands of each Landsat 8 SR image during the classification periods.Then,the areas where the NRFI values fall within the threshold range on multiple dates are extracted.Finally,all results extracted above are combined (P=max(NRFIt1,NRFIt2,...,NRFItn).With reference to a previous study [51],the threshold was obtained from the histogram of the rapeseed reference samples.The histograms can be found in Fig.S9.

Step 2: The maximum NIR value in the classification period(NIRmax) is calculated.Then,the region (F=(P and NIRmax)) is obtained by combining the NIRmax value with P.

Step 3: The maximum value (VHmax) and the coefficient of variation of VH(VHcv)from two months before the flowering stage to two months after the flowering stage are calculated.The rapeseed map (R=(F and Slope and VHmax and VHcv)) comprises the area where F,the slope map,VHmax,and VHcv intersect.This step integrates Landsat 8,Sentinel-1,and topographic slope data.

Step 4:The speckle noise is removed by connecting the domain and filling the small gaps inside the parcels [52].The thresholds obtained for the different indicators in different regions can be found in Table S3.

To test the robustness of the proposed method,we extracted rapeseed images in three scenarios based on Landsat images obtained in the flowering stage,pod stage (pod-filling),and flowering-to-pod stage for each zone.

2.6.Validation of the rapeseed map derived from APPA

Verification points and existing crop products were used as reference data for accuracy verification.We used verification points to evaluate the accuracies of the rapeseed maps in China (zones I,II,and III).The validation reference samples were randomly generated.We checked and labeled each sample and omitted points not labeled with a distinct land cover type based on field surveys and Sentinel-2 images during the rapeseed flowering stage(Fig.S10).Finally,a total of 397 rapeseed samples (1203 nonrapeseed samples),414 rapeseed samples (1187 non-rapeseed samples),and 334 rapeseed samples (968 non-rapeseed samples)were collected in zones I,II,and III,respectively,for validation(Fig.S10).The rapeseed maps of zones IV,V,VI,and VII were compared with the land cover SEC maps in France,the CROME in United Kingdom,the NCLM in Germany,and the ACI in Canada at the pixel level.The confusion matrix(Table S4),user accuracy(UA,Eq.(2)),producer accuracy (PA,Eq.(3)),overall accuracy (OA,Eq.(4)),and F-score(Eq.(5))were used to evaluate the map accuracies.TP and TN are the number of samples classified as rapeseed and nonrapeseed,respectively.FP and FN are the number of samples misclassified as rapeseed and non-rapeseed,respectively.

Table 1 Confusion matrix and evaluation metrics of rapeseed maps in zones I-III based on the validation points at several stages of growth.

Fig.2.Flowchart of this study.The main steps include(A)collecting satellite and topographic data,pre-processing data(excluding bad observations),and developing multitemporal features,(B) mapping rapeseed by the APPA method,and (C) assessing accuracy based on validation points and inter-comparison with existing products.

Fig.3.Rapeseed maps were generated based on the APPA in diverse regions and at diverse stages of growth.(A)Rapeseed maps in different zones.(B)Regional-scale views of rapeseed maps.The subfigures depict the red-outlined regions in panel A.The data for each subgraph can be accessed publicly at https://doi.org/10.5281/zenodo.6451104.

3.Results

3.1.The spectral and polarization characteristics of rapeseed

According to field surveys in zone II(Fig.1),the peak flowering dates occur on approximately April 3 based on the International Biologische Bundesanstalt,Bundessortenamt und Chemische Industrie (BBCH) scale [2,53].The green,red,and NIR reflectance values increased during the flowering stage due to yellow rapeseed petals [3,12,33].The reflectance values of swir2 were lower than those of the green band throughout the flowering period (March 20,2020 to April 20,2020) and the following month (pod-filling,April 20,2020 to May 20,2020)(Fig.1B).This spectral feature disappears in the late pod stage of rapeseed.Although the reflectancevalues of swir2 were also lower than that of the red band during the flowering stage,the difference between the red band and the swir2 band was smaller during the following month.

Fig.4.Comparison of the verification results of rapeseed maps generated by the APPA in three stages of growth and different zones (zones I-VII).PA,Producer accuracy;UA,user accuracy.

The results showed that the VH time series had a local minimum (-17.49,April 3,2020) during the flowering stages because the rapeseed petals reduced the backscattering coefficient(Fig.1E).In addition,VH reached a maximum value (-10.38,April 27,2020)during the pod stage.The intertwining of pods and stems increased the VH values.During the growing stage,the VH time series was consistent with those reported in previous studies[3,25,27,34].Note that the backscattering coefficient is affected by the incidence angle.We extracted the time series profiles of VH under different incident angles (31.68-43.10°) in the sample rapeseed parcels (Fig.S11).We found that the incidence angle had a slight effect on the backscattering coefficient during the growth of rapeseed.Furthermore,although the climatic and other natural conditions vary among different regions,the NRFI and VH time series of rapeseed both presented the same characteristics among the various regions (Fig.S7).

3.2.Rapeseed map from APPA

The spatial distributions of rapeseed generated by the APPA in multipl regions are shown in Fig.3.Many rapeseed areas were omitted from the classification results at the flowering stage owing to the lack of available images.The map of rapeseed at the pod and flowering-to-pod stage identified more rapeseed areas than the results for the flowering stage.

3.3.Accuracy assessment

We assessed the accuracies of the rapeseed maps generated by the APPA in three stages of growth based on verification points(zones I-III)and existing land cover datasets(zones IV-VII).Table 1 and Table S5 show the confusion matrix.The omission error of each zone was high during the rapeseed flowering stage,while during the pod stage there were fewer omission errors and more correctly identified rapeseed pixels.The combination rapeseed map generated from the flowering and pod stages showed the lowest omission and commission error for each zone.

Table 2 The numbers of all images obtained in multiple zones (I-VII) and at diverse stages of growth.

Figs.4 and 5 show that the classification results obtained for different stages in all regions had high user accuracies(UA ≥ 0.9).However,the producer accuracies of the Landsat images based on the flowering stages (PA ≤0.5) in all study areas except zone VI were much lower than those based on the other two stages(PA ≥0.7).The results of the F-score were similar to those of PA.The results show that the APPA method can be used to identify rapeseed from Landsat images at the pod stages even if no images are available at the flowering stages.

Fig.5.Accuracy assessment based on verification points.(A)Verification results of rapeseed maps in three stages of growth and multiple regions.(B)Regional-scale views of the verification results.The subfigures depict the red-outlined regions in panel A.

We compared the accuracy results of the APPA with those of existing land cover datasets at the pixel level.We found a marked difference between the rapeseed map obtained with the APPA method and the existing crop products in zones IV,V,and VII.The rapeseed map obtained at the flowering stage is highly consistent with the NCLM of Germany in zone VI because images were available during the flowering stages.Moreover,the rapeseed pixels obtained in the images taken during the pod stages based on the APPA show a very consistent spatial pattern with those of the existing crop products (Fig.6).

Fig.6.Spatial comparison between the rapeseed maps generated by the APPA and existing crop products.(A)Inter-comparison results at three different stages of growth in multiple zones.(B)Regional-scale views of results for inter-comparison.The subfiguresdepict the red-outlined regions in panel A.The red pixels show that the rapeseed maps and existing products do not overlap.Yellow pixels show overlap between rapeseed maps and existing products.The data for each subgraph can be accessed at https://doi.org/10.5281/zenodo.6451104.

4.Discussion

4.1.Advantages of the APPA for classifying rapeseed

We proposed the APPA to identify rapeseed based on the unique spectral characteristics of Landsat images.The APPA is quite different from the methods developed in previous studies[2,4,12,20,29,30,32,34].(1)Most of these studies used a single type of optical satellite imagery.Optical images reflect the spectral features of rapeseed but not the structural information such as geometric features and texture of rapeseed canopy.The integrated spectrum and polarization indices incorporated in the APPA are effective for distinguishing rapeseed from other land cover types.(2) Most of the methods developed in previous studies almost always rely on good-observed images obtained during the flowering period [4,29,30,35].Compared to previous studies,the APPA avoided the effects on classification accuracy of diverse flowering times and the lack of good observations at the flowering stage.The NRFI used in APPA has unique spectral characteristics during the flowering and post-flowering stages (Figs.1D,S7).That the reflectance values of swir2 were lower than those of the green band throughout the flowering stage and the following one month(pod-filling stage) (Fig.1B) offers a new opportunity to identify rapeseed in the pod stage.Thus,using the APPA method,we might incorporate all possible images that are available during the classification period,which extends approximately 60 days.Even if some images were disturbed by clouds on certain days of the flowering phase,the rapeseed area could still be identified based on good observations obtained within the flowering month (Table 2).Although good observational images were unavailable during the flowering stages for most cultivation area (Figs.7,S12),we still mapped the main planting areas using the APPA.

We compared APPA with the Colorimetric transformation and Spectral features-based oilseed Rape extraction Algorithm (CSRA)method proposed by Wang et al.[4] and the NGVI method proposed by Fang et al [54].The parameter settings in the CSRA and NGVI came from reference [4].Rapeseed was identified based on Landsat imagery using CSRA and NGVI methods.Landsat images and the validation sample were the same as those used by APPA.Non-vegetation and non-crop pixels were first masked using NDVI ≥0.3 and NIR ≥0.23 in the CSRA and NGVI methods.The accuracy and spatial details of these methods are listed in Tables 3,S6 and Fig.S13,respectively.The results showed that CSRA had high producer and low user accuracy.APPA had high user accuracy and overall accuracy.The main reason was that CSRA overestimated the area of rapeseed.Among the three methods,NGVI had the lowest accuracy.

Fig.7.Numbers of good observations obtained in different zones and under different stages of growth by excluding clouds,snow or ice,and cirrus.

We used all good observations obtained during the classification period to reduce the impacts of phenological differences on the accuracies of the resulting maps.Compared with the machine learning classifier used in previous studies [2,24,30,34],machine learning requires a large number of training samples,while no such demand exists for the APPA,and the APPA can be maximized using available imagery.Additionally,some areas in the rapeseed maps generated by the APPA have less‘‘salt and pepper”noise than that seen in the land cover SEC map in France(Fig.8A)or the ACI in Canada(Figs.8C,S15).Compared with CROME,our rapeseed maps have more detailed spatial information (Fig.8E,F).

4.2.Uncertainty analysis and implications for future studies

The accuracy and stability of the APPA maps may be affected by several factors.The accuracy of the APPA algorithm is influenced by information on rapeseed growth phenology.Accurate phenology information will reduce omission and commission error in rapeseed mapping.The phenology of rapeseed varies in different regions with climate and other factors [25,35].Some countries and regions have accurate crop phenology observations,such as data from agrometeorological stations in China [9] and Deutscher Wetterdienst(DWD)in Germany[37].However,many regions still lack continuous crop phenology records,posing a challenge for the global mapping of rapeseed.Several studies have tried to simulate the phenology of rapeseed [2,27,35].d’Andrimont et al.[2] simulated the expected phenology of rapeseed in Germany using growth degree days calculated from meteorological data,but the method they used has not been validated in other countries.Combining remote sensing,topographic and meteorological data to simulate rapeseed growth phenology may prove an effective method or improving the accuracy of APPA rapeseed global mapping.

Crop structure will affect the accuracy of crop mapping because it is associated with the complexity of the mixed pixels in the landscape[55].There are multiple types of crop structures in the study area.APPA showed better performance in areas with more homogeneous pixels.Rapeseed and wheat are the main crops in zones III.The main crops in zone III are rapeseed and rice.However,the farmland in the hilly areas of southern China is fragmented,a possible reason for the lower producer accuracy of rapeseed maps in zone III than in zones I-II.Figs.S16-18 show the crop and land cover of zones IV,V,and VII based on existing land cover products(Table S1).In addition,rapeseed is often rotated with other crops(e.g.rice and wheat) to reduce damage from pests and diseases[25,56,57].Our previous study[25]found that rapeseed crop rotation breaks are at least two years in most European countries.There may be large differences in fertility stages in some areas owing to different management practices,increasing mapping errors.

Table 3 Comparing the accuracy of rapeseed maps generated by APPA and two other methods.

It is still a challenge to identify rapeseed in areas where there are no available good observations at flowering and pod stages.The spatial resolution of the remote sensing images being used is another source of uncertainty [58-60].Mixed pixels are widespread in mountainous and hilly areas,where cropland fields are fractional and irregular [4,19].Consequently,some rapeseed parcels are too small to be recognized.More recently,as more Sentinel-2 Level-2A SR images with higher resolutions are avail-able and the Planet satellite has begun providing images with higher temporal resolutions;these images have been successfully used to map crop areas [12,61,62].The increasing number of good-quality images could further improve the robustness of the APPA and allow the mapping of rapeseed areas with higher accuracy.

Fig.8.Spatial detail comparison of rapeseed maps obtained from existing crop products and the APPA.The panels show the red-outlined regions in Fig.S14.

5.Conclusions

A new method is demonstrated for automated rapeseed mapping that exploits the specific spectral and polarization characteristics of rapeseed during the growing season and does not rely on peak flowering images.The user accuracy and overall accuracy of the resulting rapeseed map exceeded 0.9 and the producer accuracy exceeded 0.87 in three diverse regions of China.The spatial distribution of rapeseed maps in regions of France,United Kingdom,Germany,and Canada showed high consistency with existing crop maps.It provides a new opportunity for mapping long-term crop areas at large spatial scales.

Data availability

The rapeseed maps generated by the APPA algorithm in this study are available at https://doi.org/10.5281/zenodo.6451104.

CRediT authorship contribution statement

Jichong Han:Methodology,Investigation,Validation,Formal analysis,Writing -original draft,Visualization.Zhao Zhang:Conceptualization,Writing -review &editing,Project administration,Funding acquisition.Juan Cao:Formal analysis,Writing-review&editing.Yuchuan Luo:Formal analysis,Writing-review&editing.

Declaration of competing interest

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

Acknowledgments

We thank the editors and anonymous reviewers for their valuable comments.This research was funded by the National Natural Science Foundation of China (42061144003).

Appendix A.Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2022.04.013.