From parcel to continental scale -- A first European crop type map based on Sentinel-1 and LUCAS Copernicus in-situ observations

In plain words

Farm policy in the European Union needs to know which crop grows on which field. This paper builds a map of crop types for the whole EU for 2018, at a resolution of about 10 metres, using radar satellite images instead of photographs. Radar works through clouds, so the satellite can deliver usable images on a regular schedule.

The method is a teacher-and-student setup. Ground-truth points that people visited in the field (from the LUCAS Copernicus survey) say which crop grows where. A random forest, a collection of many decision trees, learns to recognise the radar signature of each crop over the growing season. The signature changes as plants grow, so the model looks at a series of ten-day summaries rather than a single image. The model is trained in two steps: first it separates broad land types, then it names the crop on the land that is farmland.

The authors then check the map three ways: against held-out field points, against farmers' declared parcels in six regions, and against official regional crop area statistics. The map is useful but imperfect, and the paper is open about where it fails. It is also an early, single-year result.

TL;DR

Two-level random forest classifiers, trained on Sentinel-1 VV and VH backscatter from January to July 2018 and on LUCAS Copernicus in-situ polygons, give a 10 m crop-type map of the EU (19 crop classes). The authors report 76.1% area-weighted overall accuracy on held-out LUCAS core points for the 19 crop classes and 80.3% for grouped crops, with weaker results for grassland, fodder crops and root crops. The code is public but has no license, so its use is limited to reading it.

Problem

Detailed, EU-wide crop type maps at parcel scale were needed for evaluating agricultural policy, but optical imagery (cloud-limited) and national declaration data (mostly not open) did not cover the whole EU consistently §1, §1.1, §1.2. The paper argues that Sentinel-1 radar is unaffected by atmospheric conditions and so is a consistent input for machine learning over a large area §1.1. It also argues that the LUCAS Copernicus module (2018) provides the first polygon-level in-situ training data with crop labels at EO-compatible detail §1.2.

The approach assumes that the 2018 season is representative: the authors note that 2018 had droughts and heatwaves and that a model trained on one year may not transfer to others §4, perspectives.

Method

Overview: the workflow is a sequence of pre-processing, two-phase classification per stratum, a post-hoc mask, and three validations fig. 1. The study area is EU-28 with a 10 km buffer [§2, study area; fig. 2].

Sentinel-1 input (§2.1, §2.1.1). S1 IW-mode GRD products are used, in VV and VH, from both ascending and descending orbits. The GEE processing produces sigma-0 at 10 m. Terrain flattening and speckle filtering are not applied §2.1.1; §4.

Compositing steps as described in the paper §2.1.1:

1. Edge masking: groups of neighbouring pixels with VV below -25 dB are removed before compositing.
2. Averaging over 10-day periods: sigma-0 is averaged in each period across all ascending and descending acquisitions, separately for VV and VH, and the average is converted to dB. The average is taken on natural (linear) values, as the order of the steps implies.
3. The VH/VV cross-polarisation ratio is also computed and averaged per 10-day period §2.1.1. It is not used in the final map §3.1.

The result is 36 composites per polarisation for 2018, independent of the number of acquisitions in each period §2.1.1.

# Build 10-day composites for one pixel  [§2.1.1] [code: S1/GEE_S1_compositing.js#L20-L93]
# Input:  scenes of (VV, VH) in dB, one per acquisition
# Output: 36 composites of (VV, VH) in dB for 2018
def composite_pixel(scenes, bin_edges):
    # Remove edge pixels with VV below -25 dB  [§2.1.1, step 1]
    scenes = mask_edges(scenes, vv_threshold_db=-25)
    # Average in linear units, not dB  [§2.1.1, step 2]
    linear = to_linear(scenes)
    composites = []
    for start, end in bin_edges:
        # All ascending and descending passes in the 10-day window  [§2.1.1]
        window = linear[(linear.date >= start) & (linear.date < end)]
        mean_linear = window.mean(axis=0)
        # Convert the mean back to dB  [§2.1.1, step 2]
        composites.append(to_db(mean_linear))
    return composites

LUCAS in-situ data (§2.2). The LUCAS 2018 core survey and the LUCAS Copernicus module are both used. The Copernicus module records land cover polygons up to 51 m in four cardinal directions around a point, with a minimum mapping unit given as 78.53 m² §2.2.2. Copernicus polygons with a level-3 legend code are used; the paper gives different totals for these (see Pitfalls).

Legend (§2.3, Table 1). LUCAS classes are regrouped into two levels. Level 1 has five land-cover classes: artificial land (100), arable land (200), woodland and shrubland (300), grassland (500) and bare land (600). Level 2 has 19 arable crop classes tab. 1. The caption says seven level-1 classes, but the table and the text list five.

Level-2 classes tab. 1, with the LUCAS codes that map to them:

*F40 with agricultural land use (U111/112/113) tab. 1.

Stratification (§2.4.2). The EU-28 is split into two strata by biome §2.4.2, fig. 2. Stratum 1 holds temperate broadleaf and mixed forests, temperate grasslands and savannas, boreal forests/taiga, tundra and temperate conifer forests. Stratum 2 holds Mediterranean forests, woodlands and shrub.

Two-phase classification (§2.4.3). Each stratum has two classifiers: level 1 separates the five land-cover classes, and level 2 names the crop for pixels of the arable class. In total four random forests are trained (two strata times two levels) §2.4.3. Random forests are chosen for robustness and low tuning sensitivity, as argued from cited work §2.4.3.

Training data (§2.4.1). The classifiers are trained on the LUCAS Copernicus polygons whose labels are mapped to the legend §2.4.1.

Feature selection (§2.4.4, §3.1). Two choices are made. The time window uses an 80/20 polygon split of the training data §2.4.4; the split for the index comparison is not stated:

• The time window: a monthly hindcast from January to December selects January to July as a trade-off between accuracy and timeliness [§2.4.4; fig. 4; fig. 5].
• The S1 indices: VV, VH, VH/VV and their combinations are compared over January to July §2.4.4; tab. 3. VV plus VH gives the best overall accuracy in the paper's comparison tab. 3.

The final input is therefore 10-day VV and VH backscatter from 1 January to 31 July: 22 dekads, which gives 44 features per pixel [§3.1; code: Classification/Classification_at_scale/s1_classify_2levels.py#L34].

Hyperparameter tuning (§2.4). For each of the four models, a random search samples 100 combinations and evaluates each with 3-fold cross-validation, 300 fits in total §2.4. The seven tuned parameters are listed in tab. S3. Two are called key: the number of trees, searched from 300 to 1200 in steps of 100, and max_features, tested as the square root or the base-2 logarithm of the feature count §2.4. Implementation: scikit-learn RandomForestClassifier and RandomizedSearchCV §2.4.

# Two-phase crop classification for one stratum  [§2.4.3] [§2.5]
# Input:  stack (H, W, 44): 22 VH then 22 VV dekads, January to July
# Output: class map (H, W) with the 19 crop codes or 100/300/500/600
def classify_stratum(stack, level1_model, level2_model):
    # Flatten to one row per pixel  [§2.5] [code: s1_classify_2levels.py#L76]
    pixels = stack.reshape(-1, 44)
    # Missing values are set to 0 in the code; not stated in the paper
    pixels = replace_nan_with_zero(pixels)
    # Phase 1: land-cover class per pixel  [§2.4.3] [tab. 1]
    labels = level1_model.predict(pixels)
    # Phase 2: crop type for pixels labelled arable (200)  [§2.4.3]
    arable = labels == 200
    labels[arable] = level2_model.predict(pixels[arable])
    return labels.reshape(stack.shape[:2])

Classification at scale (§2.5). The model is applied pixel by pixel over EU-28 on the JRC Big Data Analytics Platform, using an HTCondor environment §2.5. Stratum models and tile outputs are then merged and reprojected to ETRS89-LAEA (EPSG:3035) §2.5.

Post-processing mask (§2.5, tab. S4). The map is masked for:

• Steep terrain above 1000 m, using slope above 10° from the ALOS AW3D30 DSM §2.5. The paper's wording reads as altitude above 1000 m and slope above 10°; see Pitfalls.
• Built-up areas from the GHSL 2015 settlement map at 10 m §2.5; tab. S4.
• Permanent water from JRC Global Surface Water 2019 §2.5; tab. S4.
• Inland and coastal wetlands, bare land, moors and heathlands from CLC 2018 §2.5; tab. S4. Tab. S4 also lists ports and airports.

Accuracy assessment (§2.6). Three approaches §2.6:

1. LUCAS core points from 2018 that were not used for training §2.6.1. The paper states four quality criteria: direct and in-situ observation, parcel larger than 0.1 ha, and a homogeneous land-cover class. The accuracy metrics are area-weighted, following Olofsson et al. (2014), with a 95% confidence level §2.6.1. The estimator equations are not written out in the paper.
2. GSAA farmers' declarations from six regions (Flanders, Denmark, Centre-Val de Loire, the Netherlands, North Rhine-Westphalia and Slovenia) §2.6.2, tab. 2. Only crops covering at least 1% of a region's declared area are used. Each parcel gets the modal (majority) class of the pixels inside it, which is compared to its declared crop mapped to the LUCAS legend §2.6.2, tab. S10. Producer and user accuracy are computed per class §2.6.2.
3. Area comparison with Eurostat NUTS 2 statistics (NUTS 1 for Germany and the UK), for crop areas in 1000 ha §2.6.3. Pearson correlation is computed per crop, and relative area differences are mapped [§2.6.3; fig. 10; tab. S17].

Software and data. The classification and tuning use scikit-learn; GEE does the S1 processing and the data are then downloaded §2.1.1; §4. About 10 TB of data were transferred §4.

Hyperparameters

Random forest settings of the four final models tab. S3. The Level and Stratum columns follow the paper's table. All values are as stated in the supplementary table; the paper does not give the random search results behind them.

Results

Author-reported. Setup: EU-28, 2018, Sentinel-1 VV and VH from January to July, two-level random forests trained on LUCAS Copernicus polygons, and validated against three independent references §3.

Overall accuracy and per-class producer (PA) and user (UA) accuracy against 87,853 LUCAS core points, area-weighted tab. 5; §3.3.1:

Grouped crop classes, 87,853 LUCAS core points, area-weighted tab. 4; §3.3.1: overall accuracy 80.3%.

Feature selection for the time window and indices §3.1; tab. 3:

GSAA parcel comparison, six regions, 3,149,334 parcels, 8.21 Mha [§3.3.2; tab. 2]: the paper reports common wheat, barley, maize, sunflower and rape as the most reliable classes across regions. Rape and turnip rape reaches PA above 98% and UA above 97% in the three regions where it occurs [§3.3.2; fig. 9]. Per-region confusion matrices are in tab. S11–S16.

Eurostat NUTS 2 comparison, 220 regions §3.3.3; fig. 10: Pearson correlations of 0.93 (potatoes, §3.3.3 text) to 0.99 (rape and turnip rape). Abstract and results text differ on which crop has 0.93; see Pitfalls.

Map extent §3.2.1; tab. S5: 91.74 Mha of arable land across 19 crop classes; the paper states this as 9,174 million 10 m pixels. Woodland and shrubland covers 223.4 Mha and grassland 70.9 Mha, both pixel-count estimates §3.2.1; tab. S5.

Training data §2.4.1; tab. S1: 58,178 polygons, 2,956,889 S1 pixels, 23 class rows (see Pitfalls for the count mismatches).

Paper vs code

Repository: https://github.com/ec-jrc/eucropmap, commit bb147e3027c3ddf3f071d95c7152441422e47a6d. The GitHub owner is the JRC organisation, which employs the authors, and the README names the paper and its authors. This is an institutional repository, not a third-party reimplementation, so it is used here. The repository has no LICENSE file; the code is read for comparison only and nothing is quoted from it.

Compositing matches the paper: linear averaging, then conversion to dB GEE_S1_compositing.js. The edge mask is -25 dB as in the paper. The code also calls a connected-components step with a maxSize of 100 and a 1×1 rectangle kernel, a size limit that the paper's main text does not state GEE_S1_compositing.js.

Feature window and input order match the paper: 22 decades, with VH stacked before VV before reshaping [code: Classification/Classification_at_scale/s1_classify_2levels.py#L34-L63]. The level-2 model is applied only to pixels predicted as arable (200) s1_classify_2levels.py, as in the paper's two-phase design.

Hyperparameter search differs from the paper:

• Search budget matches: RandomizedSearchCV with 100 candidates and 3-fold CV Lucas_RF_train model-crops-all-pol-jan-jul-v7-162020.ipynb.
• The n_estimators grid is [150, 300, 500, 800, 900, 1000, 1100, 1200], not a 300 to 1200 range in steps of 100 ...notebook.
• max_features includes 'auto' in addition to sqrt and log2 ...notebook.
• min_samples_split is searched over [3, 5, 7]. The paper's Table S3 has min_samples_split = 2 for three of the four models ...notebook.
• bootstrap is searched over [0, 1]. The paper's Table S3 says True for all four models. The notebook's model fits set bootstrap=0 [code: ...notebook#L606-L607; #L1212-L1213].
• The notebook's fits show several configurations, including min_samples_leaf of 4 and n_estimators of 800 ...notebook. The fitted pickled models that produced the map are not in the repository, so Table S3 cannot be checked against the code.

Training rows: the training matrix is the set of extracted S1 rows for Jan to Jul ...notebook. The notebook does not average rows to polygon means in the cells read. The paper's Discussion says polygon averages were used §4. Which one is correct is not settled by the code.

Validation filter differs from the paper. The code keeps "in situ" observations within 100 m, observations on the point, parcels not under 0.5 ha, and LC1 cover above 75% CreateLUCAS2018validationPoints.R. The paper says parcels above 0.1 ha and a homogeneous class §2.6.1, and its commented-out criteria read PARCEL_AREA_HA >1 (with a note saying 0.1 ha) and LC1 cover above 90% §2.6.1 source comments. The three versions do not agree.

Post-processing mask: the code default for maximum DEM height is 800 m, not 1000 m maskeucropmap.py. Slope is set to zero where DEM is below the threshold, so the mask applies only above the height limit maskeucropmap.py. The slope limit default is 10° maskeucropmap.py. The run arguments used for the published map are not in the repository. The code also applies a 3 by 3 median filter to slope maskeucropmap.py. The code masks the CLC codes 123, 124, 322, 331, 332, 333, 335, 411, 412, 421, 422 and 423 maskeucropmap.py. This matches the codes in tab. S4, including ports and airports, which the main text does not mention.

Not checked: the GSAA majority rule, the NUTS 2 comparison code, the per-class F-score by country code, and the pixel-accuracy module, which implements a separate IGARSS study. The GSAA scripts were not read in full.

Pitfalls

• Polygon counts disagree. The introduction gives 58,428 polygons §1.2, §2.2.2 gives 58,426 with a level-3 legend, and §2.4.1 and the results give 58,178. Use 58,178 for the training set, as the classifier section does §2.4.1. Copernicus polygons total 63,287 §2.2.2.
• Pixel counts disagree. The results give 2,296,889 pixels §3.1. Table S1 totals 2,956,889 tab. S1. The table sums match its own rows.
• Class count of the training set. The caption of Table S1 says 21 thematic classes. The text (§3.1) and the table body have 23 rows.
• Legend levels. The Table 1 caption says the level-1 legend has seven classes. The text and the table body have five §2.4.3; tab. 1.
• Rice has almost no training data. The paper states 14 polygons and reports the class as missed in the map §4; tab. S1.
• Training sample unit is unclear. The paper's Results describe pixel time series, while the Discussion says polygon averages were used for the RF training §3.1; §4. The code matches the first reading. Specify this before reimplementing.
• Cross-validation unit is not stated. The RandomizedSearchCV folds may be pixels or polygons. Pixel folds leak neighbouring samples of the same polygon into the validation split. Use polygon-grouped folds.
• Selection split is not stated for indices. The time window uses an 80/20 polygon split §2.4.4. The index comparison in Table 3 does not state its split. The reported 79.89% is therefore a selection-set number, not a held-out test accuracy.
• Overall accuracy is reported as 78.8% and 76.1%. The results give 78.8% for end of July in the monthly hindcast §3.1. The map accuracy is 76.1% §3.3.1; tab. 5. Table S2 prints 0.766 for the same overall accuracy. Use 76.1% as the map's accuracy.
• Rape accuracy in the abstract. The abstract says rape reaches more than 96% for both user and producer accuracy. That figure comes from the GSAA parcel comparison §3.3.2. The LUCAS result for rape is 82.7% PA and 80.2% UA tab. 5. Cite the source with the number.
• Eurostat difference direction is inconsistent. The Figure 10 caption defines the difference as Eurostat minus map. §3.3.3 says common wheat and maize are underestimated by the map. §4 says common wheat is overestimated in the north and maize is overestimated there too §3.3.3; §4; fig. 10. Check the sign with the figure definition before using these statements.
• Potatoes and rye correlation. §3.3.3 says Pearson r is 0.93 for potatoes and rye. The abstract says potatoes §3.3.3; abstract.
• GSAA parcel count. The text gives 401,833 parcels for Flanders §2.6.2. Table 2 gives 409,884 tab. 2.
• Terrain mask logic. The paper's wording reads as altitude above 1000 m and slope above 10° §2.5, which is an AND. The code default is 800 m maskeucropmap.py. The paper's Discussion says steep slopes above 1000 m §4. Use the paper's wording, and check the threshold used for the published map.
• Ports and airports are masked. Table S4 lists CLC 123 and 124 (port and airport areas) among the masked classes. The text does not mention them tab. S4; §2.5.
• No terrain flattening or speckle filter. Sigma-0 is used without terrain flattening, and no speckle filter is applied §2.1.1; §4. Steep terrain and incidence effects therefore remain, especially for summer crops with bare soil.
• Season ends in July. Classes that peak later, such as sugar beet, maize and potatoes, are weak in this map. The paper states this and recommends longer windows §3.3.2; §4.
• Accuracy is a single-year result. The map covers 2018 only, a drought year. The paper says a model trained on 2018 may not transfer to other years §4.
• Cost. About 10 TB of S1 composites were moved from GEE to a local facility for classification after the season §4. Downloads from GEE for large areas need a paid Google Cloud account for non-commercial use §4.
• Validation estimators are not given. The area-weighted accuracy uses estimators from Olofsson et al. (2014). The equations are not in the paper, so the 95% intervals cannot be reproduced from the text alone §2.6.1.

Lineage

• Builds on: random forest classification [Breiman 2001; §2.4.3], the LUCAS survey and its 2018 Copernicus module §1.2; §2.2.
• Compared against: no other classifier or product is compared quantitatively. The discussion compares its seasonal findings qualitatively with two cited Belgian and French studies §4.
• Superseded by: not known from content already in this wiki.