Earth observations for official crop statistics in the context of scarcity of in-situ data

2.1In-situ data and sampling

The ground-truth data are available in the format of geo-referenced polygons drawn around the fields that were visited during the field survey in June 2018 [3]. A single crop class is assigned to each polygon. In order to obtain the training and validation samples, we automatically generated points inside the polygons within a minimum distance of 30 m from the border of the polygons and to the closest point, respectively. For crop classes with more than 50 polygons only one point per polygon was generated (Fig. 2). We then randomly selected 50 samples for training of the crop type classifiers. The remaining samples were used for validation (Table 1). Training and validation data points were reviewed for quality control. Each point was checked in the EOSTAT CropMapper administrator tool for consistency with the mean Normalized Difference Vegetation Index (NDVI) signature of a given crop category. If the NDVI signature of a given point was not consistent with the mean signal of a given crop class, then the point was removed from the dataset. RF classifiers have the advantage that they are not sensitive to outliers in the training data set. Furthermore, data points were individually reviewed through visual inspection using as very high resolution images as background in Google Earth to verify positional correctness.

Figure 2.

Detail of the CAWa ground-truth reference data available for Kashkadarya. The map shows the polygons as provided by the CAWa dataset, as well as one point per polygon generated inside each shape.

Two supervised classifiers, the RF and the DTW were trained with monthly composites of Sentinel-1 and Sentinel-2 satellite images available from the period March to October 2018. Several rounds of training were performed using different numbers of randomly sampled training data points (5–10 points for DTW, 10–50 points for RF). We then applied the trained algorithms to classify all pixels of Kashkadarya region at a spatial resolution of 10 m, except for urban areas (based on the Copernicus Global Land Cover Layers1), or pixels above a certain elevation (2500 m asl) or with steep slopes (> 5∘). For validation purposes, we looked at the classification accuracy with respect to all available validation points, as well as with respect to a subset of 30 random validation points per crop class. The subset was generated to account for the class imbalance in the validation data set (Table 1).

2.2Temporal composites

Sentinel-1 and Sentinel-2 time series data are harmonized and reduced in size by implementing a temporal aggregation. Temporal aggregation requires a temporal interval parameter to be chosen, to generate (almost) cloud-free composites from all images within the chosen temporal interval. This results in a harmonized time series of equal-size temporal intervals if the same temporal interval is applied for the entire year. With knowledge of the crop calendar, different temporal intervals for aggregation can be chosen for different periods of the year, to get a denser profile around key crop growth stages between onset and senescence. We have tried both approaches:

The geometric median was used to generate the temporal composites, as this preserves the spectral relationship between the composited bands, smoothens local spectral artefacts, and ensures consistency across scene boundaries.

As an alternative to temporal aggregation, data gaps could be filled spatially using a technique like the spectral angle mapper based spatio-temporal similarity (SAMSTS) developed by Yan and Roy in 2018 [11]. However, temporal aggregation was preferred to a gap-filling approach because data reduction was a desired outcome, considering the number of pheno-spectral input features used in this study. Based on the resulting temporal composites, the full extent of Kashkadarya region was deemed to have sufficient cloud-free observations in the span of a year, and therefore than no additional spatial-temporal gap filling was necessary.

2.3Harmonics fitting

After temporally aggregating data using a given temporal interval, data gaps and undesirable artefacts may remain in the data. While data gaps need to be filled, and artefacts need to be smoothened, the same method of time series harmonic fitting is suitable to address both issues. Also known as Fourier analysis, this method allows a complex curve to be expressed as the sum of a series of sine and cosine waves. Each wave is defined by a unique amplitude and a phase angle where the amplitude value is half the height of a wave, and the phase angle (or simply, phase) defines the offset between the origin and the peak of the wave over the range 0–2 π. Therefore, high seasonal variation in NDVI or EVI of crop pixels will be characterized by high amplitude values and phase angles [12]. Although both linear and non-linear harmonic fitting approaches have been considered to model time series for cropland mapping [13] it is not clear whether non-linear fitting improves the classification results. Therefore, this study has only applied linear fitting to model the time series, using a harmonic parameter of 4, which matches well the number of composite observations (16 for most pixel locations in the study area).

2.4Pheno-spectral feature extraction – crop reference signature

The fitted time series generated provide us with phenology of all land surfaces in the Kashkadarya region at a pixel level. To better characterize the phenology of the different surfaces, namely those of crop, phenological stages can be extracted from the time series [14]. The main advantage of extracting phenological stages is to be able to visually, as well as well as algorithmically, separate crops with distinct crop calendars. Moreover, considering the representativeness of the phenological stages extracted, the spectral values of Sentinel-1 VV channel, Sentinel-2 bands and vegetation indices corresponding to the DOY at which the phenological stages were extracted to provide additional meaningful input features to discriminate between crop types, also referred to as pheno-spectral features. The use of phenological stages alone may be subject to the inter-annual and inter-regional variation of crop calendar, but including the dynamic spectral properties of pheno-spectral features offers a more reliable identification of different crop types. The phenological stages extracted, as well as their corresponding pheno-spectral features, are summarized in Table 2. One can appreciate the resulting data reduction of using extracted phenological stages and pheno-spectral features, as opposed to raw time series (Table 3), while ensuring optimal discrimination of the targeted crop types. Considering the operational nature of the conceived system, data size and processing time are of concerns, and any data reduction that does not significantly affect the quality of the produced output is considered. The results of the extraction of the pheno-spectral signature were used to build a library of crop reference signatures.

2.5Classification algorithms

Two supervised classification algorithms were used for the generation of the crop maps: i) time-constrained Dynamic Time Warping (DTW) and ii) Random Forest (RF) supervised classification. The algorithms underwent through several rounds of training using different numbers of randomly sampled training samples (5–10 points for DTW, 10–50 points for RF). The trained algorithms were in turn used to classify all pixels of Kashkadarya region at a spatial resolution of 10 m, except for urban areas (based on the Copernicus Global Land Cover Layers), or pixels above a certain elevation (2500 m above sea level) or with steep slopes (> 5∘). For validation purposes, we looked at the classification accuracy with respect to all available validation points, as well as with respect to a subset of 30 random validation points per crop class. The subset was generated to account for the class imbalance in the test data set (Table 1). The final results obtained with the two algorithms were then compared.

2.6Time-constrained dynamic time warping

We used a DTW algorithm to simulate condition of scarce availability of in-situ data. The DTW algorithm classifies every pixel in the study area by comparing the spectral signature of that pixel to that of the reference crops contained in the crop library. The DTW measures the similarity between two temporal data sequences by identifying their optimal alignment and by producing a dissimilarity index (Fig. 3). Computing the alignment between two sequences is done recursively using the DTW matrix. The algorithm picks the smallest DTW dissimilarity value between the query pattern and the available reference patterns and attributes the corresponding reference crop label to the pixel.

Figure 3.

Computing the alignment between two sequences at hypothetical test areas (TA) TA1 andTA2. The vertical and horizontal values of the DTW matrix represent the date of an image. The alignment between two sequences is computed only for the yellow cells of the matrix, reducing the number of computations necessary (a maximum time delay, w, of 45 days is used in this example). After computing the matrix from upper left to lower right, the last element of the matrix, m[S,T], is returned, as a measure of DTW dissimilarity between the two compared sequences. Image credits: Csillik et al. [16].

DTW classification has the main advantage that only a small number of training samples are required to build the crop reference library. According to Belgiu and Csillik [15, 16] as little as 3 samples per crop type. This is a big advantage for regional and national crop type mapping, especially in countries which lack in-situ data of adequate quality. A few clean reference samples that represent the characteristic temporal pattern of the crop type are sufficient for the DTW to function. Conversely, the positional accuracy and correctness of the crop label of the training data is of paramount importance for DTW, as the algorithm is very sensitive to errors in the training data.

A time-flexible method for comparing two temporal patterns by considering their temporal distortions in their alignment, the DTW has proven to achieve better results than the Euclidean distance measure for NDVI time series clustering [17]. This flexibility is desirable for crop mapping, to deal with the intra-class phenological discrepancies caused by different agricultural practices, environmental conditions, or by different weather conditions. This is beneficial also for comparisons among samples which are more geographically distant than the pixel to be classified, as they may exhibit temporal shifts in growing patterns, while still belonging to the same crop type class for computing DTW. It is apply time constraints considering the specific seasonality of different crops. For example, comparing an element of a sequence with all other elements of another sequence leads to erroneous results when aligning a winter crop with a summer crop. Applying time constraints on time warping increases the speed of processing, while providing meaningful results. A so-called time constrained DTW implementation was therefore used in our study. The elements of two time series will be compared only if the date difference is smaller or equal to ω, which is the time constraint period which we set to 30 days, following the recommendations of Csillik et al. [16].

2.7Random forest

Random forest is an ensemble learning method for classification that operates by constructing a multitude of decision trees based on available training data as shown in Fig. 4 [17].

Figure 4.

Workflow of the random forest classification.

In a multi-band RF, the number of variables available for tree construction is n times i, where n is the number of bands used (Table 2) and i is the number of times steps. The output of the random forest is the class selected by most trees (majority voting). Using multiple deep decision trees, trained on different parts of the same training set, generates a classification model with a reduced risk of overfitting the training set. This generally increases the performance of the model, especially if a large training data set is available. The RF classifier has been extensively used to map land cover mapping and crop type mapping from Landsat images [18, 19], Sentinel 1 and 2 [20]; Sentinel-2 and the Gaofen-1(GF-1, Chinese satellite data, [21]. The main strength of the RF is that it copes well with collinearity, and is not affected by outliers, high dimensionality, and noisy features. Wang et al. [22] confirm that random forests generalizes well within regions where crop compositions and phonologies remain similar.

2.8Accuracy test

In order to compare the performance of the two classifies in producing accurate crop type maps, we have initially trained the DTW using 5 and 10 training and produced the respective crop maps. Then we have trained the RF using 5, 10, 20, 30, 40 and 50 training samples and produced the respective crop maps. In line with recommendations from Foody [23] and Congalton [24, 25], we evaluated the accuracy of the DTF and the RF through the establishment of confusion matrixes using two different settings: i) we used all of the validation data (Test A), and ii) we used only 30 validation samples (Test B). The Overall accuracy (OA), the Producer Accuracy (PA), the User Accuracy (UA), were computed for each map produced by the DTW and the RF under both Test 1 and Test 2.

The PA is the map accuracy from the point of view of the map maker (the producer). It measures how often are the ground-truth features correctly shown on the classified map, or the probability that a certain land cover of an area on the ground is classified as such. The PA is the complement of the Omission Error, PA = 100%-Omission Error. It is calculated by summing the number of reference sites classified correctly divided by the total number of reference sites for that class.

The UA is the probability that a value predicted in the map to be really belonging to that class. The User’s Accuracy is the complement of the Commission Error, UA = 100%-Commission Error. The User’s Accuracy is calculating by taking the total number of correct classifications for a particular class and dividing it by the row total.

The OA is the probability that an individual will be correctly classified by a test; that is, the sum of the true positives plus true negatives divided by the total number of individuals tested.

The OA’s obtained from the validation of the maps produced by the DTW and RF classifiers were finally analyzed in function of the number of training samples used, for both Test A and Test B.

Figure 5.

Overall accuracy of the DTW and RF classifiers as a function of the number of training samples used. In a) all validation data points were used to assess the OA. In b) only 30 validation data points were used to assess the OA. In both a) and b) the DTW scores significantly higher OA than to RF, when 5 to 20 samples are used for training. Only when 30 or more training samples are used, the RF scores higher OA compared to DTW.