Cropland Mapping over Sahelian and Sudanian Agrosystems: a Knowledge-based Approach Using Proba-v Time Series at 100-m

Early warning systems for food security require accurate and up-to-date information on the location of major crops in order to prevent hazards. A recent systematic analysis of existing cropland maps identified priority areas for cropland mapping and highlighted a major need for the Sahelian and Sudanian agrosystems. This paper proposes a knowledge-based approach to map cropland in the Sahelian and Sudanian agrosystems that benefits from the 100-m spatial resolution of the recent PROBA-V sensor. The methodology uses five temporal features characterizing crop development throughout the vegetative season to optimize cropland discrimination. A feature importance analysis validates the efficiency of using a diversity of temporal features. The fully-automated method offers the first cropland map at 100-m using the PROBA-V sensor with an overall accuracy of 84% and an F-score for the cropland class of 74%. The improvements observed compared to existing cropland products are related to the hectometric resolution, to the methodology and to the quality of the labeling layer from which reliable training samples were automatically extracted. Classification errors are mainly explained by data availability and landscape fragmentation. Further improvements are expected with the upcoming enhanced cloud screening of the PROBA-V sensor.


Introduction
In West Africa, about 70% of the population depends on cropping systems to meet their basic needs [1].Rainfed and smallholder agriculture are the two major components of the food production system [2][3][4].Nevertheless, inter-annual crop production variability, mainly due to the variability of the rainfall distribution and poor soil fertility, has led to chronic food insecurity [5].Recurrent severe regional droughts occurred in 1972-1973 and 1983-1984 and led to famines that affected millions of people [6].Even though no regional drought has occurred since 1984, localized famines recently hit some Sahelian countries [7].In addition to natural variability, socio-economic and political factors, such as population growth, social conflict and low economic development, put pressure on cropping systems [8,9].It is of primary importance to monitor land use change [10][11][12] and to develop an in-depth understanding of the detailed spatial patterns and the temporal dynamics of cropland [13].Early warning systems e.g., the Famine Early Warning Systems Network (FEWS-NET) and the Global Information and Early Warning System (GIEWS), for food security require accurate and up-to-date spatial information about the cropland to monitor food production [7].Nevertheless, mapping cropland remains challenging in this region, especially because of the agricultural landscape fragmentation, the spatial heterogeneity of the cropland, the diversity of the cropping systems and the mosaic of cropland, fallow and natural grassland [14,15].Many global land cover products delivered at different resolutions ranging from 30-m to 10-km contain cropland classes.High resolution Landsat-based (30-m) cropland maps, such as Cropland-Use Intensity (CUI), Africover and GlobeLand 30, are consistent, but the methodology prevents regular updates [16,17].Medium resolution land cover maps, such as MODIS product, ESA Climate Change Initiative (CCI) Land cover, GlobCover, Global Land Cover 2000 (GLC 2000), Global Land Cover by National Mapping Organizations (GLCnmo), Global Food Security Analysis-support Data (GFSAD), International Institute for Applied Systems Analysis (IIASA) cropland, EcoClimap II and GLC-Share [15,[18][19][20][21][22][23][24][25], offer current cropland information.Yet, most of these global land cover maps do not focus on agriculture and tend to neglect the regional specificities of Western Africa.Specific cropland maps do exist, such as the global map of rainfed cropland (GMRCA) and the global irrigated area map (GIAM), but their resolutions of 10-km are not suitable for regional applications [26,27] and are too coarse for monitoring and planning purposes.
Several studies have highlighted discrepancies in the extent and the spatial distribution of cropland within global products [16,28,29].For Africa, the MODIS product contains 30% more crop area than the GlobCover product [30].The global cropland extent fluctuates between 1.11 and 3.62 billion hectares according to different global products [26].The complexity of the Sahelian and Sudanian agrosystems increases the inconsistencies between global croplands (Figure 1).Some of the differences are explained by the spatial resolution, the year of production, as well as differences in the thematic legend classes and the cropland definition (Table 1).Table 1.Diversity of agriculture-related classes for the 10 previous global land cover maps.GlobeLand 30 is based on Landsat images (30-m); GlobCover and CCI are based on MERIS images (300-m); MODIS and GLCnmo are based on MODIS images (500-m); GLC 2000 is based on SPOT-Vegetation images (1 km); the JRC MARS map hybrid product at 250-m; and IIASA cropland, GLC Share and Global GFSAD hybrid products at 1-km.
In this comparison, the mosaic classes of cropland and natural vegetation are considered as cropland.The GlobeLand 30 cropland class considers land used for agriculture, horticulture and gardens, including paddy fields, irrigated and dry farmland, vegetation and fruit gardens [17], while MODIS cropland is defined by land covered with temporary crops followed by harvest and a bare soil period [18].Bontemps et al. [31] highlighted high year-to-year instabilities in labels, especially between cropland and natural vegetation for GlobCover and MODIS products.
Many classification algorithms have been developed to discriminate land cover from satellite data.Parametric classifiers, such as the maximum likelihood classifier, and non-parametric classifiers, such as the neural network classifier, decision tree classifiers and machine learning algorithms, are among the most popular [32,33].Due to their discriminating power, many types of neural networks have been developed [34].Decision tree classifiers transform very complex classification problems into simpler decision making processes [35].The use of machine learning algorithms is becoming increasingly important, yielding high accuracies for land cover classification and often outperforming other algorithms [36][37][38][39][40].
Taking advantage of the spatial resolution of PROBA-V, this research aims at developing a new method to map cropland in the Sahelian and Sudanian agrosystems where global products remain inconsistent and cropland mapping challenging [29].The fully-automated methodology limits the use of field data and could produce regional cropland on a regular basis for crop monitoring.The relative importance of spectral-temporal features used in the method is studied.Then, the resulting 2014 Sudano-Sahelian Cropland map is validated with photo interpretation on very high resolution imagery and compared to other global products.The spatial distribution of errors is assessed, and potential explanatory variables are studied.Finally, an error analysis focuses on challenging areas where most discrepancies between global products occurred.

Study Site
This research focuses on the Sahelian and Sudanian agrosystems in West Africa with the area of interest ranging from 17°W-23°E to 9°N-18°N, comprising a region more than 4300-km long and 1000-km wide (Figure 2).This region is known for a strong north-south gradient linked with the movements of the Intertropical Convergence Zone (ITCZ).Climatic conditions range from arid and desert regions in the north to sub-humid areas in the south through semi-arid areas.The study area is characterized by two clearly distinct seasons: a rainy season and a dry season.A flat topography with low plateaus and valleys characterizes the landscape.The Sahelian region is the area between isohyets 250-mm and 500-mm.Sowing occurs at the beginning of the rainy season.The growing period is short (60-90 days).The start and the length of the growing period are positively correlated with the onset of precipitation, which is more variable than the ending date of the rainy season [5,41].The Sudanian region is considered in this study as bounded by isohyets 500-mm and 1000-mm, although there is no broad agreement on its definition.The length of the rainy season in the Sudanian region equals the length of the dry season, resulting in growing period of more than 150 days between May and November [42] (Figure 3).The fields are usually very small and of various shapes, with many fields with an area less than 1 ha.

Data
PROBA-V is the successor satellite of SPOT-Vegetation and was launched in May 2013 in order to fill the gap between SPOT-Vegetation and Sentinel-3.While it was not in the initial specification of the mission, this sensor currently offers global coverage every five days at a 100-m spatial resolution in four wavebands, blue, red, near-infrared (NIR) and mid-infrared (SWIR), thanks to its central camera.PROBA-V 100-m daily surface reflectance has been composited on a weekly basis using a mean composite approach.The time series of eleven months (48 composites) from 19 March 2014 to 12 February 2015 was prepared and includes the growing season.The data availability mainly depends on cloud/shadow contamination and the non-five-day global coverage at 100-m in the lower latitudes due to the 517-km swath [43].False cloud detections are systematically observed over some high brightness areas, such as desert (Figure 4).However, these areas do not contain crops and will not degrade the classification results.In the study area, 37 out of the 48 seven-day composites were on average usable (77%), making it possible to understand the temporal behavior of the main crop types.However, clouds are more frequent in lower latitudes and during the vegetation growth season,.These factors could limit classification quality in those regions (Figure 5).
A labeling layer map was built from GlobeLand 30 [17] with the specific contribution of ESA CCI Land Cover [19] for the irrigated cropland.The labeling layer is the reference map from which the training dataset is extracted for the supervised cropland classification.GlobeLand 30 offers a global land cover of 10 classes at 30-m based on Landsat images for 2010.The accuracy of this global dataset was first computed with an independent validation dataset of Zhao et al. [44] and results in an OA of 88%.GlobeLand 30 was then validated with a set of 7 high resolution images.A Pareto boundary analysis showed the performance of the 30-m spatial resolution and the better balance between omission errors and commission errors for GlobeLand 30.The area under the Pareto boundary defines the performance of the spatial resolution, as the smaller the area, the more performant the resolution (for more details, see [45]).This area is two-fold smaller for GlobeLand 30 than for the JRC product and three-fold smaller than for the MODIS and the CCI products.The CCI Land Cover project delivers three global land cover products, the most recent relating to the 2008-2012 period.In the study area, agriculture is mainly rainfed, while irrigation occurs at specific locations, e.g., the Niger basin.The semi-arid conditions lead to strong spectral contrast between rainfed cropland and irrigated cropland.Therefore, a separate class for irrigated cropland has been created with all single pixels belonging to the crop class in GlobeLand 30 and irrigated crop in the ESA CCI Land Cover Map.The final labeling layer contains eleven classes: the ten of GlobeLand 30 plus the irrigated crop class from the CCI Land Cover.

Methodology
The challenge of cropland mapping in West Africa is dealing with the very strong latitudinal gradient in terms of seasonal dynamics, as well as land cover.A knowledge-based approach identifies relevant temporal features throughout the growing period to optimize cropland capture.A support vector machine (SVM) classifier was selected for its good performances compared to other classification algorithms' classifiers [36,38].First, the 2014 Sudano-Sahelian Cropland map was compared to existing global products to study discrepancies.Then, classification errors were related to potential explanatory variables, such as location, valid cloud-free observation and landscape characteristics.Finally, classification errors were also studied in regions where global products disagreed on cropland.The method developed is fully automated, i.e., meaning no human action during the whole processing from images to cropland map.

Cropland Classification
The development of cropland mapping in Sahelian and Sudanian regions proceeds in a three-step method: (i) extraction of the temporal features from filtered PROBA-V time series; (ii) efficient local training based on trimmed data; and (iii) classification using SVM (Figure 6).

Extracting Temporal Features from PROBA-V
A knowledge-based approach permits a better classification of cropland based on relevant temporal and spectral features.The features were surface reflectance composited at specific events of the growing season when crops are expected to behave differently than other land covers do.Among the cropland discrimination studies [46][47][48][49], five temporal features were selected: the maximum of the red band (max red), the minimum and maximum of the Normalized Difference Vegetation Index (NDVI) and the increasing and decreasing slopes of the NDVI profile (Figure 7b).Soil preparation practices, such as tillage and sowing, clear the land surface contrasting with naturally-vegetated areas and resulting in higher reflectance in the red band.This timing corresponds to the maximum of red reflectance.During the rainy season, the vegetation development and the growth rate are expected to be superior in cropland than in non-cropland areas.This was captured by the NDVI maximum and minimum and the increasing slope of the NDVI.At the end of the season, crops are harvested, leading to a sharp decrease in vegetation, which is captured by the decreasing slope of NDVI.These features are specific to the vegetation cycle and do not relate to a specific time, thus making it possible to take into account the high latitudinal gradient and the diversity of agriculture practices (crop type, management).PROBA-V reflectances were smoothed with a Whittaker filter to fill the gaps and remove the residual noise [50].Some studies confirmed the adequacy of this smoother for time series filtering [51,52].These smoothed reflectance values were extracted at a specific stage of temporal features for each pixel of the study area.The max red feature for instance was derived from the reflectances (4 bands) of the date at which the red reflectance reached a maximum for each pixel.Hence, neighboring pixels can belong to different dates, but correspond to the same phenological event.This led to twenty inputs for classification (five temporal features and four spectral bands) called spectral-temporal features (Figure 7).

Trimming and Local Training
Due to the ability of SVM to handle small training sets [53], a trimming procedure was performed to remove misleading data due to the absence of an error-free and up-to-date labeling layer and recent land cover change.The trimming used spectral values of all classes to identify mislabeled pixels in a specific class.This was done by clustering so as to select pixels with the highest probability of belonging to a specific cluster.A random sample was used to select pixels from two sets of pixels: (i) pixels of class i; (ii) pixels from all classes except i.The clustering classification identified clusters of class i that belonged to a different class.Following the clustering, the probability of each cluster to be of class i was assessed.The probability of each cluster to be of class i corresponds to the purity of the cluster, i.e., the number of pixels of class i divided by the total number of pixels for a given cluster.Samples of group i belonging to the cluster with higher than a 75% probability of being classified as i were selected as the training dataset for this specific class i.This approach was repeated for each class.The class proportion in the training set was maintained similar to that initially found in the labeling layer.Additional cleaning was applied to discard pixels belonging to a class that deviated strongly from the distribution (95% of the confidence interval) of that class in at least one of the four bands of the five features.

SVM Classification
The trimmed data training set was used to calibrate the SVM classifier (for more details, see [54]).SVM searches the optimal separating hyperplane in a multidimensional space that maximizes the margin between the defined plan and the data.The one-against-one method was used to perform multiclass-classification with i classes on the SVM.Therefore, i(i − 1)/2 binary classifiers were trained, and the final class was obtained through a voting procedure.The SVM classifier identifies the hyper-plan using support vectors and margins.The cost-support vector classification was trained with a Gaussian radial basis kernel function to ensure a high level of automation.The γ parameter of the SVM was set at ten to optimize performance and computation time [40,55].

Handling the Spatial Gradient and the Landscape Diversity
The combination trimming-training-classification was applied on a moving window of 1 × 3 • to handle the latitudinal gradient.In fact, the vegetation dynamic in cropland and natural vegetation can fluctuate strongly with latitudes [56][57][58][59].To limit the boundary effect between grid cells, the training was performed on a larger grid cell.A buffer of 1 • in longitudinal direction and of 0.3 • in latitudinal direction was chosen as a compromise between processing time, high latitudinal gradient and recovery rate (40%).To mosaic smoothly the moving window results, five staggered cropland classifications were compiled by majority voting to obtain regional consistent cropland.

Relative Importance of Spectral-Temporal Features
A preliminary analysis of the respective spectral-temporal features' contribution was completed to identify the subset of the twenty spectral-temporal features that maximizes the classification accuracy.Twenty spectral-temporal features could increase both noise and processing time.A random forest algorithm was run to calculate the spectral-temporal feature importance by turning off one feature (while others remained constant) and computing the resulting loss in accuracy through the mean decrease in the Gini index (GI) [60].This analysis carried out for four classes of crop proportions calculated on quantiles: very low crop proportion (<1.4%), low crop proportion (1.4%-4.4%),medium crop proportion (4.4%-17.8%)and high crop proportion (17.8%-72.2%)(Figure 8).The distribution of crops varies across a latitudinal gradient with a higher proportion of crops between 10 • and 13 • of latitude.The feature importance analysis was computed for each grid cell within a crop proportion class.Then, a random sample of five grid cells per crop proportion class was used to perform the classification with the best selected features (Figure 8).The best features were selected according to a decision rule; the mean decrease in GI must exceed the mean of the mean GI decrease for a given grid cell.The accuracy of SVM classification with the best selected features was compared to SVM classification with all features using a Student test computed on the F-score.

Validation
A two-stage sampling was designed to collect the validation sample.The first sampling units were 40 grid cells within the 126 grid cells of the study area.The second sampling unit was a stratified random sampling, as recommended by Wagner and Stehman [61], within the selected grid cells.The strata corresponded to the eleven land cover classes of the labeling layer.The sample size allocated to each stratum was twenty pixels.All of the random 100-m by 100-m pixels automatically generated were visually assessed with 2014 images in Google Earth.Each random pixel was labeled by visual interpretation in one of the two classes: crop and non-crop, depending on whether more than 50% of the pixel was covered by the specific class.The availability of 2014 images was the main limiting factor in building this dataset.A set of 2315 points with 1616 non-crop and 699 crop validation samples had finally been collected.In total, 30% of the validation points belong to the crop class.This proportion (30%) matches the magnitude of the crop proportion in this region observed on the global product (crop proportion ranging from 13% to 42%).
Typically, the accuracy of a map is assessed by measuring the degree of agreement between the classification output and the validation dataset.A confusion matrix was used to derive the metrics of classification accuracy [62].Liu et al. [63] recommended amongst fourteen category-level and twenty map-level accuracy measures the user's accuracy (UA), producer's accuracy (PA) and overall accuracy (OA) as primary accuracy indices.The OA was computed as the ratio of the sum of all correctly-classified pixels on the sum of all validation pixels.The F-score was calculated as a combination of PA and UA for a specific land cover class (Equation ( 1)).The UA for a specific class is the ratio between correctly-classified pixels and all pixels classified as the specific class.The PA for the specific class is the fraction of correctly-classified pixels and all ground truth-specific class pixels.
Because of this research focus on cropland, the F-score was only computed for crop class and is called the F-score in the following sections.

Error Analysis
The accuracy assessment described in Section 3.4 provides a global performance evaluation, although it is well established that classification accuracy varies across space and that errors are not equally distributed spatially [59,[64][65][66].Eight potential explanatory variables were proposed to explain the classification accuracy computed by OA and the F-score: latitude and longitude of the grid cell center, availability of cloud-free data and five landscape metrics indices.Within the five landscape metrics, some are specific to crops, such as crop proportion, Matheron index and crop fragmentation, while others are related to the overall landscape, such as fragmentation and entropy.These landscape features were measured on the labeling layer in order to explain the spatial distribution of errors.The correlation and the importance rank computed with the random forest algorithm were calculated for all of the potential explanatory variables.
The perimeter area ratio for all classes of land cover is used as a proxy of landscape fragmentation: where i represents a specific class of land cover and n the total number of land cover classes.The Shannon entropy index measures the diversity of the landscape due to the number of land cover classes and the proportions of those classes (Equation ( 3)) [67].In general, this index increases when the number of classes is higher and when the proportions of all existing classes are equal.
where P i is the proportion of a class i.
The Matheron index measures the ratio between the total outer perimeter of crop cells edges (sum for all patches) and the product between the area of crop and the total area [68]:

Spectral-Temporal Feature Importance
First, the diversity of the features seems to complement them according to the cropland proportion.The maximum in red reflectance and the minimum NDVI seem the two most discriminant features in higher crop proportion regions (Figure 9).These features refer to the start of the growing period when differences between cropland and natural vegetation are high due to land preparation.In fact, high crop proportion areas correspond to regions with higher amounts of rainfall, which makes natural vegetation greener at the beginning of the season when crops are considered as bare soils due to land preparation.All temporal metrics seem important in one or another crop proportion class without a specific distinction for any one of them.The mean decrease in GI was computed for the 20 spectral-temporal features for each of the four crop proportion classes.When a given spectral-temporal feature is left out, the higher mean decrease in GI compared to the GI obtaining using all 20 features leads to the higher importance of this spectral-temporal feature.Regardless of the crop proportion classes, as expected, the blue band contributes the least in cropland discrimination due to the impacts of aerosols and atmospheric scattering [69] (Figure 9).The SWIR band is of higher importance than the NIR band, while the red and the SWIR bands are the two most important bands for the classification.Following the feature importance analysis, the best features were selected (SVM select) to compare classification accuracy with regards to the classification using all of the features (SVM all).Due to the absence of validation points in the very low crop proportion area, only three crop proportion classes were considered, and five grid cells were randomly selected in each crop proportion stratum.OA and F-score were computed for the two classifications: SVM select and SVM, both by crop proportion classes.SVM select classification slightly outperforms SVM for low and medium crop proportion strata, while SVM excels in the high crop proportion landscape (Figure 10).Standard deviations are rather larger for SVM select than for SVM.The two classifications were compared using a Student test on the F-score.The assumptions of the test were verified: the normality of the variable and the equality of the variances of both samples were confirmed by a Shapiro-Wilk test (p-value = 0.735) and a Fisher test (p-value = 0.872), respectively.The two F-score values are considered equal (p-value = 0.873), which avoids rejecting the null hypothesis of the equality of samples.Based on these results and to make the method as generic as possible for all agricultural landscapes, the SVM classification was applied at the regional scale using all twenty spectral-temporal features.

Qualitative Analysis of 2014 Sudano-Sahelian Cropland map
The visual comparison of the 2014 Sudano-Sahelian Cropland map and the labeling layer, mostly derived from GlobeLand 30, highlights the actual contribution of the PROBA-V classification.The images used for this analysis are a combination of Digital Globe images acquired during 2014, Spot 5 and RapidEye imagery for the years 2012 and 2013.Images 1 and 2 are respectively a World View 2 of 26 June 2014 and a GeoEye of 14 June 2014.Images 3-7 are Spot 5 of the 2013 season, while Image 8 is a RapidEye image acquired on 29 January 2012 (Figure 11).
Cropland boundaries are more precise in the labeling layer, probably due to higher spatial resolution, while the non-cultivated lithosols are better delineated by the 2014 Sudano-Sahelian Cropland map (Image 1, in Figure 11).Roads are better defined by the labeling layer, which, however, includes more noise and incorrect non-crop area.The village patches are larger in the 2014 Sudano-Sahelian Cropland map, but unlike the labeling layer, it does display the full extent of the village (Image 2 in Figure 11).The contrasted landscape shown in image 3 with cultivated valleys and rocky plateaus yields an accurate classification in both cases.The labeling layer with the higher spatial resolution even captures the forest in the very bottom of the valley as non-crop (Image 3 in Figure 11).In a more complicated landscape, the 2014 Sudano-Sahelian Cropland map failed to discriminate bare soil areas (in pink in Image 4) and dark surfaces (in black), while some of the cropland is mapped as non-crop in the labeling layer (Image 4 in Figure 11).In arid areas, errors observed in both datasets are probably due to the labeling layer, which drives the classification training.Crops in these arid areas are very sparse, and discrimination between bare soil and crop is challenging (Image 5, in Figure 11).This is confirmed by Image 7 in the north of Burkina Faso, where crops alternate with bare soil.The labeling layer missed considerable cropland areas, while these are correctly mapped by the 2014 Sudano-Sahelian Cropland map.Hence, in this area, spectral-temporal features allow the better capture of crop areas thanks to the temporal dimension of the PROBA-V time series (Image 7, in Figure 11).In irrigated regions, the labeling layer failed to capture all cropland areas.The 2014 Sudano-Sahelian Cropland map captures more crops, but also more noise (Image 6, in Figure 11).Finally, Image 8 sums up what is often found in these illustrations.The 2014 Sudano-Sahelian Cropland map seems to better capture all of the cropland areas, but has more difficulty in drawing exact boundaries and in capturing small landscape features, such as forest in the bottom valley.Conversely, the labeling layer is sometimes too restrictive and misses some of the cropland (Image 8, Figure 11).To conclude, difficulties in differentiating between cropland and fallows, as well as misclassification of shrubland, bare soil and desert area as cropland remain in the final product.
Supervised classification with visual determination of the region of interest (ROI) was performed on each of these sites.Wall to wall validation between high resolution images and GlobeLand 30 or 2014 Sudano-Sahelian Cropland map was assessed.Proximity between the commission-omission point of the 2014 Sudano-Sahelian Cropland map and the Pareto boundary implies a high performance of the algorithm classification (Figure 12).To better quantify the visual interpretation of the 2014 Sudano-Sahelian Cropland map and the labeling layer, quantitative analysis of the accuracy of the 2014 Sudano-Sahelian Cropland map map was assessed.The produced cropland was also compared to other existing global products.

Accuracy of the Cropland Map and Comparison with Existing Global Products
The overall accuracy computed from the 2315 validation points reaches 84% and the F-score for crop class 74%.The confusion matrix shows better PA and UA for the non-crop class: 89% compared to 74% for crop classes.UA and PA are identical for both classes (Table 2).The same validation dataset was applied on the crop/non-crop maps derived from the 10 previous global products.All products were resampled to the 100-m resolution using the majority method.A high variability of accuracies was observed regardless of the production year of the map (Figure 13).
Table 2. Confusion matrix for the regional cropland map using the 2315 validation points.Overall accuracy (OA), producer's accuracy (PA) and user's accuracy (UA) are computed in percent.The 2014 Sudano-Sahelian Cropland map and the GlobeLand 30 product almost reach the same performance levels (OA: 84% vs. 86%; and F-score: 74% vs. 75%, respectively).Other products studied show substantially lower accuracy indices.This confirms the adequacy of the labeling layer and the efficiency of the methodology in exploiting recent annual remote sensing datasets.Similarities of the accuracy indices of GlobeLand 30 and our cropland product must be balanced against disagreement between those products (Figure 14).Eighty five percent of pixels are labeled identically by GlobeLand 30 and the 2014 Sudano-Sahelian Cropland map, while 41% of pixels classified as crops in one or another are mapped as crops in both products.As expected, because of the higher OA and F-score of GlobeLand 30 for the entire validation dataset (Figure 13), the OA in the conflict region (414 validation pixels) reaches 54% for GlobeLand 30 and 46% for the 2014 Sudano-Sahelian Cropland map.
However, one might argue that these measures are excessively harsh, as the size of the validation pixels for the >100-m products is smaller than the products' pixel sizes.They might be revised when considering a validation pixel size that matches the products.Nonetheless, the variation in OA when validating a product at 300-m in place of one at 100-m with the 100 squared meter validation pixels is only of 2%.This number rise to 6% when a 1-km product is used, which still cannot explain the total discrepancies between global land cover products.It was computed by aggregating the labeling layer at 100-m, 300-m and 1-km, defining the label of each pixel as crop or non-crop when more than 50% is covered by the specific class and, finally, computing accuracy.

Spatial Distribution of Errors
Grid cells including more than twenty validation points, i.e., 38 grid cells, were selected to estimate locally the OA and the F-score.The eight possible explanatory variables were also computed for each grid cell.In total, 41% of the variance of the OA is explained by the eight potential explanatory indices for only 21% of the F-score (Table 3).OA and F-score are strongly linked to fragmentation; the lower the fragmentation, the more accurate the classification.Entropy is negatively correlated with both F-score and OA.This means that the lower the number of classes and the higher the proportion in the classes, then the higher is the accuracy.As expected, the accuracy is highest when few classes are present with small fragmentation and a high proportion in each class.Spatial heterogeneity is confirmed as an important factor in driving land cover map accuracy [70,71].
Looking at the spatial distribution of latitude impacts both the F-score and the OA: it increases the accuracy in northern regions where less clouds are present.Conversely, longitude is correlated with OA, so that accuracy rises when going west (Table 3).This might be due to the fact that isohyets move south when going east, resulting in greener and more vegetated areas in the west of the study area.Hence, misclassification issues in less vegetated areas, such as desert or bare soil, could lead to lower accuracy indices.Specific cropland characteristics do not affect the F-score more than OA in terms of correlation, which was unexpected.Nevertheless, the importance rank computed by random forest highlighted the relevance of crop-specific indices for cropland classification (Table 3).The Matheron index and crop proportion are negatively related to accuracy because a higher fragmentation and a higher number of crop patches make classification difficult due to the number of mixed pixels at the edges [64].Data availability affects OA and F-score.The encountered errors were partially explained by the high fragmentation of the landscape, and the low latitude linked with lower data availability.
Very few validation points were available in high latitudes due to the absence of high resolution imagery in these regions.15).This regression (R 2 = 0.4124) was then applied to the entire study area to predict the OA (Figure 16).According to the model, a degradation of OA is observed in the southeast of the study region.

OA and F-Score in the Disagreement Region of Global Products
To further identify the accuracy of the products, it is of paramount interest to focus the analysis on challenging areas where most discrepancies between maps are observed.A global agreement map built from the global products (see Table 1) ranges from zero, where all products agree on the non-crop class, to ten, where all products agree on the crop class.A value of five represents the highest disagreement between products as five products classified the pixel as crop and five classified the pixel as non-crop.OA and F-score were computed for each of the eleven classes of agreement using the validation points belonging to that class.The OAs of all global products, except GlobeLand 30 and the 2014 Sudano-Sahelian Cropland map, have a similar behavior (Figure 17).The highest accuracy of the majority of products is met in regions labeled as non-crop by all products, and accuracy then drops until values around 0.45, when most products disagree.Finally, OA increases to yield an accuracy of 0.8 when all products agree on crop class.GlobeLand 30 and the cropland map produced in this study have a nearly identical profile, which could be summarized as a slight decrease in accuracy when agreement on crop arises.F-score increases when agreement on crop occurs.GlobeLand 30 and the 2014 Sudano-Sahelian Cropland map behave very differently to others.The F-score of the new product and GlobeLand 30 starts with a higher value in regions where many products agree on non-crop.This means that these two products correctly label crops in regions where almost all global products mislabeled crop as non-crop.q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.4 0.6 0.8 OA q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00

Fragmentation of the Landscape
The proportion error defined by [72] represents the proportion by which individual classes are over-or under-estimated at a lower resolution than the initial one (Equation ( 5)): where p 0 is the original proportion of the crop class at 30-m and p r is the estimated crop proportion at resolution r (90-m or 300-m in this study).The 90-m cropland map is obtained by aggregating the 30-m for three pixels and the 300-m for 10-pixels.The label of this obtained pixel is crop when more than 50% of the area is classified as crop.The 90-m products are consider to be very close to the 100-m resolution, as available for the PROBA-V sensor.In low and medium crop proportion classes, the underestimation of crops is of 5% for the 90-m products and more than 25% for the 300-m products (Table 4).The advantage of 100-m in a highly fragmented landscape is useful and can fill the gap between high resolution at medium resolution.

Discussion
The proposed classification method of PROBA-V 100-m time series yields an overall accuracy of 84% similar to that of GlobeLand 30 from which training data were derived.Both products overpass all other global products in terms of accuracy.The accuracy of GlobeLand 30 ought to be linked to its spatial resolution of 30-m and a highly interactive verification procedure [17].GlobeLand 30 is produced with a POK-based (pixel object knowledge) methodology, a pixel and object-based method with knowledge.The knowledge-based interactive verification procedure of GlobeLand 30 was developed through a web service and relies on intensive manual corrections.The fully-automated method presented aims at overcoming this limitation by enabling the production of cropland maps on an annual basis.Three main factors can explain the accuracy of the 2014 Sudano-Sahelian Cropland map: the 100-m spatial resolution of the time series, the methodology itself and the quality of the labeling layer.The hectometric spatial resolution makes it possible to resolve fine landscape patterns, especially in those regions where field size often does not exceed 1 ha.Similarly, a recent study has highlighted significant better performance using 100-m PROBA-V data than 300-m data for crop type classification [73].
Misclassifications are mainly explained by low data availability due to high cloud cover during vegetation growth (25% of OA explained by data availability) and by landscape fragmentation (39% of OA explained by fragmentation).Remaining errors in the training dataset are minimized by the stringent trimming procedure limits, from which 30% of the best points were retained.Highly fragmented landscape is better captured by high resolution images as demonstrated with Pareto boundaries.The Sentinel 2, 5-day coverage at 10-m in 13 spectral bands, will probably offer a good combination of spatial-spectral and temporal resolution.Landsat 8, 16-day coverage at 30-m in 11 spectral bands, offers a spatial resolution able to catch the landscape fragmentation of the study region.Nevertheless, the temporal coverage might be too small to catch vegetation behavior during the vegetation growth.Remaining errors in the training dataset are minimized by the strong trimming procedure limits.Discrepancies between the GlobeLand 30 product and the 2014 Sudano-Sahelian Cropland map can partly be linked to the time production of the map, the spatial resolution, the input data and the interactive verification for GlobeLand 30.
Spectral-temporal features were used to fully exploit the PROBA-V time series.Nevertheless, remaining noise in the PROBA-V surface reflectance image could partly explain some misclassification.The noise might be related to undetected clouds and shadows.PROBA-V cloud screening uses both blue and SWIR bands.Due to the time lag between bands (12 s between NIR and SWIR), two separate masks are built and then merged [74].Threshold rules are applied on the SWIR (blue) band to detect cloud with an additional check on the blue (SWIR) band.These threshold values are static and applied globally, leading to misclassifications of bright surfaces as clouds [43].Furthermore, borders of large clouds are poorly detected, and the remaining haze effect affects the reflectances.Some cases of cloud omission, others of cloud commission, are also present.In addition, part of the noise is related to the spectral temporal features themselves.Spectral-temporal features are based on extreme values and are thus more sensitive to noise, as noise itself is characterized by extreme values [66].Region-dependent and dynamic cloud screening are currently being invested by the PROBA-Vegetation team and could significantly improve the cloud mask [75].This could improve feature extraction and thus classification by significantly decreasing noise in the time series.If no significant improvements are made to the cloud product delivered, we believe that a stronger temporal filter should be applied prior to performing the Whittaker filter in order to identify undetected clouds.

Conclusions
This fully-automated method has delivered the first cropland at 100-m for the Sahelian and Sudanian region using the PROBA-V sensor.The product reaches an overall accuracy of 84% and an F-score of 74% for the crop class.This accuracy indices are significantly larger than major global products (except GlobeLand 30 which was used as labeling layer) thanks to the spatial resolution of 100-m, the crop-specific knowledge-based methodology and the adequate choice of the labeling layer.The diversity of the spectral-temporal features chosen to take advantage of the entire PROBA-V time series and its full spectrum seems to adequately complement each others for cropland discrimination in heterogeneous crop proportion landscapes.Misclassifications observed in the 2014 Sudano-Sahelian Cropland map are partially explained by undetected clouds, haze and shadows, data availability and fragmentation of the landscape.
Further improvements are expected with the upcoming enhanced cloud screening.The recent launch of the European Sentinel-2 sensor offers new possibilities for monitoring land and vegetation with its spatial resolution of 10-m and its 13 spectral bands.The methodology is scale independent and could be applied on Sentinel-2 time series to better capture the landscape specificities in the Sudano-Sahelian regions.Annual cropland maps could be updated on a regular basis using this fully automated method in order to monitor cultivated area extension and abandonment over a large area.

Figure 1 .
Figure 1.Cropland over Sahelian and Sudanian agrosystems on ten global products showing the high variability of cropland extent.

Figure 2 .
Figure 2. Study site with the validation sample (blue points) and regional grid (in grey) with isohyets of 250-mm, 500-mm and 1000-mm.The background image corresponds to the surface reflectances in infrared, mid-infrared and red bands for the dates corresponding to the maximum of red reflectance (maximum red feature).

Figure 3 .
Figure 3. Calendar of seasons and agricultural management.The study region contains two clearly distinct seasons with crop cultivation during the rainy season and land preparation at the end of the dry season.

Figure 4 .
Figure 4. Percentage of PROBA-V valid cloud-free observation for the study area during the whole time series.False detection of clouds is systematically observed over the desert area.The valid cloud-free observation percentage decreases with latitude.

Figure 5 .
Figure 5. Percentage of data available for each degree of latitude ranging from 9°to 18°.Low data availability is observed in low latitudes and during vegetation growth (from July to October).

Figure 6 .
Figure 6.Three-step methodology for cropland classification: (i) extraction of the temporal features; (ii) local training based on trimmed data; (iii) classification using SVM.

Figure 7 .
Figure 7. (a) Spectral response of the four bands (blue, red, NIR, SWIR) of the PROBA-Vegetation instrument; (b) Representation of the five temporal features (minimum NDVI, maximum NDVI, increasing slope, decreasing slope and maximum red).

Figure 8 .
Figure 8. Crop proportion of each grid cell and randomly-sampled grid cells for SVM classification with the best selected features.A higher crop proportion is observed in lower latitudes.

Figure 9 .
Figure 9. Respective contribution of each spectral-temporal feature for crop discrimination for the different cropland density classes.The dotted red line corresponds to the decision rule for selecting the best features.

Figure 10 .
Figure 10.Comparison of the accuracy for SVM classification with all features (SVM) and selected features (SVM select).The whiskers represent the standard deviation of the five selected grid cells for each crop proportion class.

Figure 12 .
Figure 12.Three examples of Pareto boundary, commission and omission errors for GlobeLand 30 and 2014 Sudano-Sahelian Cropland map.Pareto boundaries were also computed for a 10-m product and a 300-m product.(a) Site in center Mali ; (b) Site in Sikasso region (Mali South East); (c) Site in South Mali.

Figure 13 .
Figure 13.Comparison of the 2014 Sudano-Sahelian Cropland map with the ten previous global products showing high variability in the accuracy indices and better performance for the 2014 Sudano-Sahelian Cropland map and GlobeLand 30 used as the labeling layer.

Figure 14 .
Figure 14.The 2014 Sudano-Sahelian Cropland map derived from PROBA-V at 100-m and the comparison with GlobeLand 30 derived from Landsat imagery acquired around 2010.

Figure 15 .
Figure 15.OA and F-score predicted by the multiple regression compared to the OA and F-score observed.

Figure 16 .
Figure 16.Prediction of OA over the study area based on a multiple regression using eight explanatory variables related to spatial localization, data availability and landscape characteristics.

Figure 17 .
Figure 17.Evolution of OA and F-score with agreement computed on ten global products.The X axis corresponds to agreement between global products with 0 = all products classify pixels as non-crop and 10 = all products classify pixels as crop; N corresponds to the number of validation points corresponding to a given class of agreement.

Table 3 .
Correlation and importance rank for Overall accuracy (OA) and F-score computed on eight potential explanatory variables.Multiple Linear Regression to Explain OA and F-Score A multiple linear regression based on eight explanatory variables explains 41.24% of OA variance and 21.05% of F-score variance (Figure

Table 4 .
Proportion error for two cases (initial resolution = 30-m; and lower resolution = 90-m or 300-m) within the four crop proportion classes.