Land Cover Classification

Land cover classification turns satellite and aerial images into maps that say what covers the ground: trees, crops, grass, water, buildings, bare soil, snow. Land use is the related question of what people do with the ground, such as pasture, an airport or a park. The two are often mapped together, and many class schemes mix them. For computer vision the task is image classification when one label goes to a whole image patch, and semantic segmentation when every pixel gets a class. What makes it a subject of its own is everything around the classifier. The classes come from a mapping standard with rules about minimum areas and percentage cover. The labels usually come from other maps or from interpreters looking at sharper images, and both carry errors. A map’s accuracy is a statistical estimate made from a probability sample, not a score on a held-out split. Neighbouring pixels are so alike that a random split of training pixels can overstate accuracy by a wide margin.

The sensors, missions, resolutions, vegetation indices, mixed pixels and cloud masks behind the task are described under remote sensing. The worked example below trains a pixel classifier on a synthetic multispectral scene, compares random and spatially blocked cross-validation with the true accuracy of the map, and estimates class areas from a stratified sample.

Problem Definition

Given co-registered images of an area, predict for each mapping unit one class from a fixed legend. The mapping unit is a pixel, a patch or an object such as a field. Three choices define the problem before any model is trained.

• The legend. A list of classes with definitions. Definitions are often thresholds on cover: ESA WorldCover’s tree cover class is any area “dominated by trees with a cover of 10% or more” (WorldCover Product User Manual, v2.0). Two maps with different legends cannot be compared pixel by pixel without a crosswalk between classes.
• The minimum mapping unit. The smallest area a map records as its own patch. CORINE Land Cover records areas of 25 ha and linear features 100 m wide (Copernicus Land Monitoring Service), while the Dynamic World annotators drew polygons with a minimum of 50 × 50 m, five by five Sentinel-2 pixels (Brown et al., 2022). A pixel classifier has no minimum mapping unit unless one is imposed afterwards.
• The reference period. Land cover changes with seasons and years. A map is valid for a date or a year, and its reference data must come from the same period.

Input and Output

The input is usually a stack of bands from one sensor at one date, a time series of such stacks, or several sensors together. Optical bands are often combined with synthetic aperture radar backscatter, elevation and climate layers. ESA WorldCover 2021, for instance, uses Sentinel-2 surface reflectance, Sentinel-1 backscatter, a digital elevation model, position and climate variables (WorldCover Product User Manual, v2.0).

The output is a raster of class codes in a map projection, often with per-class probabilities. Dynamic World publishes nine probability bands, one per class, that sum to one, plus a label band holding the most probable class (Brown et al., 2022). Probabilities let users pick their own thresholds or compute composites over a chosen date range. A scene classifier instead outputs one label, or several labels for multi-label patches, per image chip.

Variants

• Scene classification. One or more labels per patch, as in the single-label EuroSAT and the multi-label BigEarthNet benchmarks described under remote sensing.
• Pixel-level mapping. Dense prediction over large areas, the form of operational products such as WorldCover and Dynamic World.
• Annual against near-real-time maps. Most products are annual. Dynamic World instead classifies every Sentinel-2 image whose metadata gives a cloudy-pixel percentage of at most 35% as it arrives, from 27 June 2015 onwards; its authors describe this as a new kind of product (Brown et al., 2022). The Earth Engine catalog lists the collection as still updated, as of October 2026.
• Time-series and parcel classification. Crop types are separated by how they develop over a season more than by any single image. Sainte Fare Garnot et al. (2020) classify agricultural parcels from 24 Sentinel-2 dates between January and October 2017, with labels from the French Land Parcel Identification System.
• Land use against land cover. Land use often cannot be read from spectra alone: a grassed area may be a pasture, a park or a sports field. CORINE’s legend contains land use classes such as airports and sport and leisure facilities (Copernicus Land Monitoring Service), which depend on spatial context and on what the area is used for.

Class Schemes

Each major map defines its legend for a purpose, and a model trained on one map’s labels learns that map’s definitions. The schemes below are described from their official pages and documents, as of October 2026.

Scheme	Coverage and resolution	Classes	Source of the map
CORINE Land Cover	Europe (39 countries for 2000–2018), vector and 100 m raster	44 classes in three levels (5, 15, 44)	Interpretation of satellite images by national teams, updated every six years
ESA WorldCover	Global, 10 m, maps for 2020 (v100) and 2021 (v200)	11 classes defined with the FAO Land Cover Classification System	Gradient-boosted trees on Sentinel-1 and Sentinel-2 features, with expert rules
Dynamic World	Global, 10 m, one map per Sentinel-2 image	9 classes	Fully convolutional network on single Sentinel-2 images

CORINE Land Cover. The European Commission began the CORINE programme in the 1980s to harmonize national land cover maps, and the first CORINE Land Cover dataset was produced in 1990. The inventory has 44 thematic classes and was updated for 2000, 2006, 2012 and 2018 (Copernicus Land Monitoring Service). The nomenclature has three levels: five headings at the first (artificial surfaces, agricultural areas, forest and semi-natural areas, wetlands, water bodies), 15 at the second and 44 at the third. Its definitions mix cover and use. Forest is defined by trees normally higher than 5 m with a canopy closure of at least 30%, while classes such as port areas or green urban areas are uses. Change between two dates is mapped as its own layer with a 5 ha minimum unit, and the service warns that differencing two status layers will not reproduce the change layer, because the units differ. The service’s technical summary gives a thematic accuracy target of at least 85% for each edition, marked as “probably not achieved” for 1990 and “achieved” for 2000. As of October 2026, the 2024 edition is scheduled for publication in the first quarter of 2027.

ESA WorldCover. WorldCover’s 2020 map was released in October 2021 with a global overall accuracy of 74.4%, and the 2021 map on 28 October 2022 with 76.7% (WorldCover Product User Manual, v2.0). Its 11 classes are tree cover, shrubland, grassland, cropland, built-up, bare or sparse vegetation, snow and ice, permanent water bodies, herbaceous wetland, mangroves, and moss and lichen. They are defined with the Land Cover Classification System of the UN Food and Agriculture Organization, a hierarchical system that lets a legend’s detail be adjusted to the information available. The two maps were made with different versions of the algorithm, and the user manual says that most changes between them are due to the algorithm rather than to the ground. The product is free under the Creative Commons Attribution 4.0 licence (WorldCover website).

Dynamic World. The nine classes are water, trees, grass, flooded vegetation, crops, shrub and scrub, built area, bare ground, and snow and ice. Brown et al. (2022) chose them after reviewing earlier legends, including the USGS Anderson system, and kept them close to the land use categories of the IPCC good practice guidance (forest land, grassland, cropland, wetland, settlement and other) so that the data could be used for carbon accounting. The classes were applied through dense polygon annotations of homogeneous areas, not through percentage-cover rules. The dataset is licensed under CC BY 4.0 with a required attribution (Earth Engine Data Catalog).

Other schemes appear with benchmarks. DeepGlobe uses seven classes following the Anderson classification, among them rangeland and forest land defined by at least 20% tree crown density (Demir et al., 2018). LoveDA uses seven classes from a Chinese national survey standard (Wang et al., 2021). Brown et al. (2022) also list older global products: the NASA MCD12Q1 500 m dataset, the ESA Climate Change Initiative 300 m dataset and the Copernicus Global Land Service 100 m dataset. WorldCover’s algorithm is based on the last of these (WorldCover Product User Manual, v2.0).

Challenges

Label sources and their errors

Labels for land cover come from four main sources, and each brings its own errors.

• Existing maps. BigEarthNet’s labels come from CORINE Land Cover 2018: 43 of its 44 third-level classes appear in the archive (our count of their Table 2). The inventory was produced by national reference centres under the European Environment Agency, which interpreted the classes from the texture, pattern and density of objects in satellite images (Sumbul et al., 2019). A map with a 25 ha minimum unit cannot describe smaller patches, so a 1.2 km BigEarthNet patch, 144 ha, holds only a few map units; small features inside it are absorbed into their surroundings (our reasoning from the two sizes). The map used was also a preliminary version. Clasen et al. (2025) report that later updates of CORINE 2018 corrected wrong and missing labels, and they show a patch of sea that BigEarthNet had labelled urban fabric. EuroSAT used the European Urban Atlas to find its patches and then checked them by hand (Helber et al., 2019).
• Interpretation by people. Annotators label images, often helped by sharper imagery. Dynamic World’s labels were drawn by 25 annotators with photo-interpretation or remote sensing experience on about 4,000 tiles and by 45 non-experts on about 20,000 more (Brown et al., 2022). On the 409 held-out tiles, the non-experts agreed with the expert consensus on 77.8% of pixels. The non-experts’ producer’s accuracy for grass was 22.4% against the consensus, and their shrub and scrub labels often covered pixels the experts had labelled trees, crops or bare ground (Brown et al., 2022, Table 5). The authors therefore smoothed each one-hot label towards a uniform distribution, by 0.2 for experts and 0.3 for non-experts, so that the training targets expressed their confidence in each group.
• Administrative records. Farmers’ yearly declarations label about 10 million parcels each year in France alone, as Sainte Fare Garnot et al. (2020) note; their own crop labels came from the French Land Parcel Identification System records.
• Field and inventory data. Karasiak et al. (2022) took forest types from the French national forest inventory’s stand map. Its dates ran from 2006 to 2015 while their images were from 2016, so they removed stands whose cover had changed with an NDVI threshold and eroded each polygon by 20 m to avoid mixed pixels at the edges.

Three kinds of error recur across these sources (our summary). The legend of the label source may not match the target legend. Geometry may not match, through generalized polygons, small misregistration or shifts in the sharper imagery used for interpretation. And the date may not match. The WorldCover validation team found visible shifts between Google and Bing basemaps at its sample sites and handled them by also accepting the majority label of each sample’s 3 × 3 neighbourhood (WorldCover Product Validation Report, v2.0).

Other challenges

• Ambiguous classes. Shrubland, grassland and wetlands are hard for people and models alike. WorldCover 2021’s user’s and producer’s accuracies are 49.1% and 46.9% for shrubland and 30.5% and 44.6% for herbaceous wetland, against 80.0% and 91.9% for tree cover (WorldCover Product Validation Report, v2.0). Dynamic World’s model also did worst on the classes its annotators were least sure of (Brown et al., 2022).
• Mixed pixels and scale. At 10 to 30 m many pixels mix classes, which biases hard classification and area estimates.
• Geographic shift. Spectra of a class vary with climate, soil, season and practice, and a model trained in one region can degrade in another. LoveDA was built to measure this between urban and rural scenes (Wang et al., 2021). It is dataset bias in geographic form.
• Spatial dependence. Nearby pixels are more alike than distant ones, which inflates accuracy when training and test pixels are close (below).
• Clouds and season. Optical classifiers see the ground only on clear days. Operational products composite many dates or mask clouds per image.

Approaches

Spectral, temporal and spatial features

Classical land cover mapping classified each pixel from its band values. Ratios and normalized differences such as NDVI, described in the remote sensing article, were added as features. Because classes such as cropland and grassland differ more in how they change over a year than at any date, operational systems summarize time series. WorldCover v200 builds 10-day median composites from cloud-masked Sentinel-2 images, fills gaps by linear interpolation, and computes the 10th, 50th and 90th percentiles and the interquartile range of each band and vegetation index. It does the same for Sentinel-1 VV and VH backscatter after a multitemporal speckle filter, and adds altitude, slope, position, biome and yearly climate statistics: 115 features per pixel (WorldCover Product User Manual, v2.0).

Spatial information enters in three ways. Neighbourhood statistics such as texture describe the pixels around a target. Object-based approaches first segment the image into regions, for instance with superpixels, and classify the regions. Convolutional networks learn spatial features directly. Karasiak et al. (2022) note that such contextual features can improve classification, but that the same spatial dependence also makes accuracy estimates optimistic when it is ignored in evaluation.

Classical classifiers

Before deep learning, practice relied on decision trees, random forests and support vector machines on hand-crafted features, as Zhu et al. (2017) summarize for hyperspectral data. These methods remain in operational use. WorldCover trains gradient-boosted decision trees (CatBoost) on its 115 features: one model on all features, one without Sentinel-1 and one without Sentinel-2. Expert rules then combine the three and use auxiliary layers such as OpenStreetMap and settlement maps to correct known confusions, such as bare ground against built-up areas and mangroves against trees (WorldCover Product User Manual, v2.0). Sainte Fare Garnot et al. (2020) explain the lasting popularity of random forests among practitioners by their speed and small memory use.

Deep learning

Patch classifiers adapted from image classification label scene chips. Dense mapping uses fully convolutional networks and encoder–decoders such as U-Net and DeepLab; Zhu et al. (2017) review fully convolutional and deconvolutional networks on aerial benchmarks such as ISPRS Vaihingen and Potsdam. Dynamic World is an example at global scale (Brown et al., 2022). It applies one fully convolutional network to each Sentinel-2 Level-1C image, with all bands except B1, B8A, B9 and B10 resampled to 10 m. The reflectances are log-transformed and mapped through a sigmoid to spread their skewed distribution. The authors describe the network as almost 100 times smaller than U-Net and DeepLab v3+. A second, mirrored network learns to synthesize the input bands back from the class probabilities, a weak supervision pathway that the authors use to compensate for deficiencies in the human labels. Training adds per-pixel weights for class balance and down-weights the centres of labelled polygons, where annotators often missed small details.

For time series, Sainte Fare Garnot et al. (2020) replace the convolutions of earlier convolutional–recurrent models with a pixel-set encoder that treats the pixels of a parcel as an unordered set, followed by temporal self-attention. On their dataset of 191,703 parcels in one Sentinel-2 tile in southern France, they reported 94.2% overall accuracy and 50.9 mean per-class IoU, against 92.5–93.8% and 42.1–48.1 for the other deep models they re-implemented, over 5-fold cross-validation (their Table 2; their text gives the margins as 0.4 and 3 to 9 points). Self-supervised pretraining on unlabelled archives, the topic of the remote sensing article’s section on foundation models, is a further route when labels are scarce (transfer learning and self-supervised learning).

Datasets and Benchmarks

The table summarizes benchmarks for land cover with the protocols their papers define. The remote sensing article gives EuroSAT’s and BigEarthNet’s sizes, sources, licences and baselines; they are repeated here only as far as the protocol needs.

Benchmark	Imagery	Labels	Split protocol	Metric
EuroSAT (2017/2019)	Sentinel-2, 13 bands, 64 × 64 patches	10 classes, located with Urban Atlas, checked by hand	Random 80/20, class by class	Overall accuracy
BigEarthNet (2019)	Sentinel-2, 120 × 120 patches at 10 m	Multi-labels from preliminary CORINE 2018 (43 classes)	Random 60/20/20 in the baselines	Precision, recall, F1, F2
reBEN, BigEarthNet v2.0 (2024/2025)	Sentinel-1 and Sentinel-2 pairs, 549,488 patches	Updated CORINE 2018, 19-class nomenclature, pixel maps and multi-labels	Geographic split within each tile	Average precision, F1 (macro and micro)
DeepGlobe land cover (2018)	50 cm RGB, 1,146 images of 2448 × 2448 px	7 classes, professional annotators	803/171/172 images	Mean IoU over classes
LoveDA (2021)	0.3 m Google Earth images, 5,987 tiles of 1024 × 1024 px	7 classes	Spatially separate districts; urban and rural domains	Mean IoU

BigEarthNet and reBEN. Among the issues Clasen et al. (2025) list with the original BigEarthNet, three concern labels and evaluation: labels from a preliminary version of CORINE 2018, patches processed with an older atmospheric correction, and a recommended split, assigned on a grid, that placed training, validation and test patches next to each other. reBEN keeps the same patch locations but labels them from the updated CORINE 2018 map with the 19-class nomenclature of BigEarthNet’s later version, adds pixel-level reference maps, and drops patches with less than 75% of pixels labelled. Its split divides each Sentinel-2 tile into an outer frame (training), an inner frame (validation) and a central square (test) in a 2:1:1 ratio. This maximizes the distance between training and test patches and keeps overlapping patches of neighbouring tiles in the training set. With Sentinel-2 alone, a ResNet-50 reached a macro average precision of 70.72% and with both sensors 70.76%, averaged over three training runs (Clasen et al., 2025, Table I). These numbers are not comparable with results on the original split (our note).

DeepGlobe. The land cover track of the DeepGlobe 2018 challenge used DigitalGlobe Vivid+ images of 50 cm pixels over 1,716.9 km². Agriculture covers 56.76% of the labelled pixels and water 3.74%. Annotators were asked to label any instance of a class larger than about 20 × 20 m. The metric is the IoU of each class averaged over the six classes; a seventh class, unknown (clouds and others), is excluded from the score. The organizers’ baseline, a modified DeepLab with a ResNet-18 backbone and class weights, reached an IoU score of 0.433 (Demir et al., 2018).

LoveDA. The images cover 536.15 km² of Nanjing, Changzhou and Wuhan in 18 administrative districts, nine urban and nine rural (Wang et al., 2021). The authors cite Tobler’s first law of geography, that near things are more related than distant ones, as the reason they split training, validation and test sets so that they are spatially independent. The dataset also defines domain adaptation tasks between urban and rural districts. In the rural-to-urban task, a model trained on rural districts reached a mean IoU of 31.27% on urban test districts, against 46.79% for the same kind of model trained on urban districts (Wang et al., 2021, Table 6). The README states that LoveDA may be used for academic purposes only, under CC BY-NC-SA 4.0, with Google Earth’s terms applying to the images (as of October 2026).

Product validation sets. Operational maps are evaluated on reference samples of their own. Dynamic World’s 409 test tiles, each labelled by three experts and one non-expert, were drawn as ten tiles per biome per hemisphere and are archived on Zenodo; the training labels are on PANGAEA (Brown et al., 2022). WorldCover’s validation dataset, built for the Copernicus Global Land Service, is not released to the public (WorldCover Product Validation Report, v2.0).

Evaluation Metrics

The error matrix and its accuracies

The remote sensing literature calls the confusion matrix an error matrix and orients it the other way from much of machine learning: map classes are usually the rows and reference classes the columns (Finegold et al., 2016). Three accuracies are read from it.

• Overall accuracy is the proportion of area classified correctly, the probability that a randomly chosen location on the map is correct. It is the accuracy of the map as a whole, weighted by area.
• User’s accuracy of class $i$ is the proportion of the area mapped as $i$ that is $i$ in the reference data: what a user of the map can expect where it says $i$. One minus it is the commission error.
• Producer’s accuracy of class $j$ is the proportion of the reference area of $j$ that the map labels $j$. One minus it is the omission error.

User’s and producer’s accuracy are the precision and recall of each class, and Dynamic World’s tables label them “Precision/User’s” and “Recall/Producer’s” (Brown et al., 2022). The difference from machine learning lies in the weighting. A test set in machine learning is a sample of convenience. A map’s error matrix should be built from a probability sample of the map, with each sampled unit weighted by its probability of being sampled, so that the entries estimate proportions of area (Finegold et al., 2016).

For stratified random sampling with the map classes as strata, Finegold et al. (2016) give the estimated proportion of area in map class $i$ and reference class $j$ as

$$\hat{p}_{ij} = W_i \, \frac{n_{ij}}{n_{i\cdot}},$$

where $n_{ij}$ is the number of sample units in that cell, $n_{i\cdot}$ the number sampled in map class $i$, and $W_i$ the proportion of the map’s area in class $i$. Overall, user’s and producer’s accuracies are then

$$\hat{O} = \sum_{j=1}^{q} \hat{p}_{jj}, \qquad \hat{U}_i = \frac{\hat{p}_{ii}}{\hat{p}_{i\cdot}}, \qquad \hat{P}_j = \frac{\hat{p}_{jj}}{\hat{p}_{\cdot j}},$$

with $q$ classes and row and column sums marked by dots. The guide follows Olofsson et al. (2014) for these formulas and for their variances, from which 95% confidence intervals are 1.96 standard errors. It also states that the kappa coefficient, though often reported, is not recommended, because its use has been questioned in many articles. Overall accuracy is not always representative of individual classes, so both user’s and producer’s accuracy of every class need to be considered.

WorldCover’s validation follows this design at a larger scale (WorldCover Product Validation Report, v2.0). More than 21,000 primary sampling units of 100 × 100 m were drawn by stratified random sampling over 149 strata, independently of any land cover map. Each unit was divided into 100 cells of 10 m, labelled by about 30 regional experts on the Geo-Wiki platform. The confusion matrix was weighted by inverse inclusion probabilities, and accuracies and confidence intervals were estimated for a stratified one-stage cluster design. Dynamic World’s authors explain why they did not do the same: maps generated from every image vary in quality with the image, so design-based accuracy representative of the whole continuously updated collection is hard to define. They report agreement with expert labels on their test tiles instead (Brown et al., 2022).

Area estimation

A map’s pixel counts are a biased estimate of class areas whenever omission and commission errors do not cancel. Finegold et al. (2016) recommend estimating the area of class $k$ from the reference column of the weighted error matrix,

$$\hat{p}_{\cdot k} = \sum_{i=1}^{q} \hat{p}_{ik},$$

and multiplying by the total area. (Their equation 6 prints the summation index as $j$; the sum runs over the map classes $i$.) Combining map and reference data gives a more precise estimate than the reference sample alone, so even a map of modest accuracy can support precise area estimates if the reference data are collected carefully. The guide also warns that comparing two maps does not measure change. Forest change is often less than 1% of an area, so with two maps that are each 99% accurate, the apparent change can be attributed to error between the two maps. WorldCover’s global matrix illustrates the bias: tree cover makes up 32.7% of the mapped area but an estimated 28.5% of the reference, consistent with the slight overestimation of trees that the validation report describes.

Agreement between people

Reference labels made by people are themselves uncertain, which caps the agreement any model can show. Dynamic World’s model agreed with the expert consensus on 73.8% of pixels, close to the 77.8% agreement of the non-experts with the same consensus. On the pixels where all three experts agreed, the model’s agreement rose to 88.4% (Brown et al., 2022, Tables 5–7). On the same test tiles, other products agreed less: 66.3% for the 100 m Copernicus Global Land Service map of 2019 and 61.6% for the 300 m ESA Climate Change Initiative map of 2018. Only the regional LCMAP 2017 map, at 75.0% on the 48 tiles it covers in the United States, scored higher. The authors note that these comparisons needed crosswalks between legends, that the other maps were annual, and that their dates did not match the test images.

Dense benchmarks borrowed from computer vision report intersection over union per class and its mean, as DeepGlobe and LoveDA do, and multi-label patch benchmarks report F1 and average precision.

Spatial Autocorrelation and Data Splits

Remote sensing data are spatially autocorrelated: nearby pixels are more alike than distant ones. When training and test pixels are drawn from the same areas, a test pixel usually has a near-duplicate neighbour in the training set, and the score measures interpolation inside known areas rather than prediction elsewhere. This is the spatial form of the leakage that a patient-level split prevents in medical image classification.

Three studies measured the size of the effect.

• Pixel classification. Karasiak et al. (2022) mapped coniferous and broadleaf forest from Sentinel-2 data with a random forest in one French department. Non-spatial cross-validation overestimated overall accuracy by 8 to 14 points, depending on the strategy, compared with the model’s accuracy on two neighbouring departments. A spatial leave-one-out cross-validation over forest stands, which removes all training samples within a buffer around each test sample, matched the accuracy on the distant departments (82.1% in cross-validation against 80.4% and 82.2%). The buffer radius, 19.74 km in their text (19.64 km in the caption of their Fig. 9), was set from the distance at which Moran’s I, a measure of spatial autocorrelation, became negligible. Holding out whole forest stands reduced the bias but did not remove it.
• Deep segmentation. Kattenborn et al. (2022) segmented tree species in drone images from 47 sites with a convolutional network. Taking block cross-validation over whole sites as the realistic estimate, validation on random hold-out tiles overestimated performance by up to 28%, and the optimism grew as training sets became smaller.
• Regression. Ploton et al. (2020) predicted forest biomass in central Africa with a random forest on MODIS and environmental variables. Random 10-fold cross-validation gave an $R^2$ of 0.53; a spatial cross-validation over 44 clusters of reference pixels gave 0.14, and for test pixels more than about 100 km from the nearest training pixel the predictive power was virtually null.

Spatial cross-validation estimates how well a model predicts in areas it has not seen. It is not the same as the accuracy of a map, a distinction Kattenborn et al. (2022) draw with reference to Wadoux et al. (2021), who argue that spatial cross-validation is not the right way to evaluate map accuracy. In our reading, the two questions call for different tools. A method’s ability to generalize is tested with spatially separated or geographically distant data. A map’s accuracy is estimated from a probability sample of the map, independent of the training data, as in the previous section. Benchmarks have begun to build spatial separation into their splits: reBEN by distance within each tile and LoveDA by district. Grouped and blocked splits are one family of the methods described under cross-validation and the train, validation and test split.

Example: Random and Spatially Blocked Splits

The script builds a synthetic multispectral scene with known land cover and trains a pixel classifier on clustered reference data, the way labels from field visits or annotated polygons tend to come. It compares three estimates of the classifier’s accuracy with the true accuracy of the resulting map: 5-fold cross-validation with pixels assigned to folds at random, 5-fold cross-validation with folds formed by spatial blocks, and a simple random sample of 500 map pixels. A second part estimates accuracy and class areas from a stratified random sample with the formulas above.

The scene has 256 × 256 pixels, nominally 10 m each, divided into 80 regions (a Voronoi partition of random points), each given one of five classes. Four bands named blue, green, red and near infrared carry illustrative class means, not measured spectra. Each region gets its own spectral offset, standing for differences between fields, soils or crop varieties, and a smooth field adds a regional gradient; pixel noise comes last. The reference data are 30 windows of 9 × 9 pixels, at most 2,430 pixels (fewer where windows overlap). The classifier is a random forest with 100 trees on the four band values. Blocked cross-validation groups pixels by 64 × 64 pixel blocks. All parameters were fixed before the first run, and the script below is that first run unchanged. It runs in about 15 seconds with NumPy 2.3.5, SciPy 1.16.3 and scikit-learn 1.7.2, and two runs gave identical output.

import numpy as np
from scipy import ndimage
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import KFold, GroupKFold

names = ["water", "tree", "crop", "grass", "built"]
means = np.array([[0.06, 0.05, 0.03, 0.02],    # water        (blue, green,
                  [0.03, 0.05, 0.03, 0.30],    # tree cover    red, NIR;
                  [0.05, 0.08, 0.06, 0.35],    # cropland      illustrative
                  [0.06, 0.09, 0.09, 0.28],    # grassland     reflectances)
                  [0.12, 0.13, 0.14, 0.20]])   # built-up
H = W = 256

def make_scene(rng):
    seeds = rng.integers(0, H, (80, 2))                  # 80 Voronoi regions
    yy, xx = np.mgrid[0:H, 0:W]
    d = (yy[..., None] - seeds[:, 0]) ** 2 + (xx[..., None] - seeds[:, 1]) ** 2
    region = d.argmin(axis=-1)
    region_cls = rng.choice(5, 80, p=[0.10, 0.25, 0.30, 0.25, 0.10])
    cls = region_cls[region]
    offset = rng.normal(0, 0.02, (80, 4))                # field-to-field variation
    smooth = np.stack([ndimage.gaussian_filter(rng.normal(size=(H, W)), 20)
                       for _ in range(4)], axis=-1)
    smooth *= 0.015 / smooth.std(axis=(0, 1))            # regional gradient
    img = means[cls] + offset[region] + smooth + rng.normal(0, 0.01, (H, W, 4))
    return img, cls

def clustered_sample(rng, n_clusters=30, half=4):
    cy, cx = rng.integers(half, H - half, (2, n_clusters))
    idx = [(y + dy) * W + (x + dx) for y, x in zip(cy, cx)
           for dy in range(-half, half + 1) for dx in range(-half, half + 1)]
    return np.unique(idx)

def rf():
    return RandomForestClassifier(n_estimators=100, random_state=0, n_jobs=1)

def cv_accuracy(X, y, splits):
    correct = 0
    for tr, te in splits:
        correct += (rf().fit(X[tr], y[tr]).predict(X[te]) == y[te]).sum()
    return correct / len(y)

print("seed  random-CV  blocked-CV  map (all pixels)  random sample of 500 (95% CI)")
for seed in range(10):
    rng = np.random.default_rng(seed)
    img, cls = make_scene(rng)
    X_all, y_all = img.reshape(-1, 4), cls.ravel()
    tr = clustered_sample(rng)                           # 30 labelled 9 x 9 windows
    X, y = X_all[tr], y_all[tr]
    block = (tr // W // 64) * 4 + (tr % W) // 64         # 64 x 64 px spatial blocks
    acc_rand = cv_accuracy(X, y, KFold(5, shuffle=True, random_state=seed).split(X))
    acc_block = cv_accuracy(X, y, GroupKFold(5).split(X, y, block))
    model = rf().fit(X, y)
    pred = model.predict(X_all)
    rest = np.setdiff1d(np.arange(H * W), tr)
    acc_map = (pred[rest] == y_all[rest]).mean()
    srs = rng.choice(rest, 500, replace=False)           # probability sample
    p = (pred[srs] == y_all[srs]).mean()
    ci = 1.96 * np.sqrt(p * (1 - p) / 500)
    print(f"{seed:4d}  {acc_rand:9.3f}  {acc_block:10.3f}  {acc_map:16.3f}  "
          f"{p:.3f} +/- {ci:.3f}")
    if seed == 0:
        keep = (pred, y_all, rng)

# Part 2: accuracy and area from a stratified random sample (seed 0's map)
pred, ref, rng = keep
n_i = 50
Wt = np.bincount(pred, minlength=5) / pred.size          # mapped area proportions
nij = np.zeros((5, 5))
for i in range(5):
    s = rng.choice(np.flatnonzero(pred == i), n_i, replace=False)
    nij[i] = np.bincount(ref[s], minlength=5)
p_hat = Wt[:, None] * nij / nij.sum(axis=1, keepdims=True)   # FAO (2016) eq. 2
OA = np.trace(p_hat)
UA = np.diag(p_hat) / p_hat.sum(axis=1)
PA = np.diag(p_hat) / p_hat.sum(axis=0)
true_prop = np.bincount(ref, minlength=5) / ref.size
print(f"\nstratified sample: {n_i} px per mapped class; estimated OA {OA:.3f}, "
      f"true OA {(pred == ref).mean():.3f}")
print("class   user's  producer's   map area  estimated area  true area  (proportions)")
for k in range(5):
    print(f"{names[k]:6s}  {UA[k]:6.3f}  {PA[k]:10.3f}  {Wt[k]:9.3f}  "
          f"{p_hat[:, k].sum():14.3f}  {true_prop[k]:9.3f}")

seed  random-CV  blocked-CV  map (all pixels)  random sample of 500 (95% CI)
   0      0.971       0.853             0.902  0.908 +/- 0.025
   1      0.962       0.837             0.843  0.820 +/- 0.034
   2      0.974       0.815             0.841  0.858 +/- 0.031
   3      0.969       0.901             0.893  0.864 +/- 0.030
   4      0.965       0.853             0.866  0.858 +/- 0.031
   5      0.964       0.830             0.872  0.874 +/- 0.029
   6      0.960       0.815             0.873  0.854 +/- 0.031
   7      0.960       0.808             0.833  0.832 +/- 0.033
   8      0.965       0.852             0.863  0.866 +/- 0.030
   9      0.967       0.858             0.844  0.874 +/- 0.029

stratified sample: 50 px per mapped class; estimated OA 0.886, true OA 0.905
class   user's  producer's   map area  estimated area  true area  (proportions)
water    1.000       1.000      0.106           0.106      0.106
tree     0.980       0.846      0.136           0.158      0.171
crop     0.840       0.842      0.304           0.303      0.294
grass    0.800       0.889      0.300           0.270      0.275
built    0.980       0.926      0.154           0.163      0.154

Random cross-validation is optimistic in every scene. Across the ten scenes it reports 96 to 97% accuracy, while the maps are between 83% and 90% correct: an overstatement of 7 to 13 points. A pixel in a random test fold almost always has training pixels from the same 9 × 9 window around it, which share the region’s offset and the local gradient, so the forest can recognize the field rather than the class. The magnitude resembles the 8 to 14 points that Karasiak et al. (2022) measured, but that is a coincidence of the chosen parameters, not a prediction for real data (our reading).

Blocked cross-validation is close but noisy, and slightly pessimistic here. It differs from the map accuracy by between −5.8 and +1.4 points and is lower in eight of the ten scenes. When a whole block is held out, the model is tested on fields it never saw, which is closer to the task of mapping the rest of the scene. Each fold, though, is trained on only about four fifths of the reference pixels, while the map comes from a model trained on all of them, and the windows in one block can contain classes or fields that are rare in the others. With 30 clusters, the estimate rests on few independent units, which explains its spread (our reasoning).

A probability sample of the map estimates the map. The 500 randomly chosen pixels give estimates whose 95% intervals contain the true map accuracy in nine of the ten scenes; in scene 9 the true value lies 0.001 below the interval. This is the design-based route: the sample is independent of how the training data were collected, and its uncertainty follows from the sampling design. The interval is the normal approximation for a binomial proportion (our choice of formula).

Stratified estimates correct the map’s areas. In scene 0 the map underestimates tree cover (13.6% of pixels against 17.1% in truth) and overestimates grassland (30.0% against 27.5%). Fifty sampled pixels per mapped class are enough to move both estimates close to the truth: 15.8% and 27.0%. The accuracies say why: tree cover’s high user’s and low producer’s accuracy mean that where the map says tree it is almost always right but it misses some trees, and grassland has the lowest user’s accuracy, so other classes are mapped as grass. In this part the true values and the strata include the training pixels (at most 2,430, under 4% of the scene), so its true overall accuracy (0.905) differs slightly from the map accuracy in the first table (0.902), which excludes them.

The scene is idealized. Classes have fixed means, regions are convex, and the reference labels are perfect. Real maps add label errors, mixed pixels and geographic shift, which all lower accuracy and none of which this example includes.

Applications

Land cover maps are inputs to environmental reporting, planning and science. The CORINE service lists environmental monitoring, land use planning, climate change assessments and emergency management, and policies in agriculture, transport and spatial planning (Copernicus Land Monitoring Service). Dynamic World’s authors chose their legend to support carbon stock and greenhouse gas estimates (Brown et al., 2022), and the FAO guide addresses the forest area and forest change estimates that countries report on (Finegold et al., 2016). Crop type maps support agricultural subsidy control: Sainte Fare Garnot et al. (2020) note that the European Union needs yearly crop maps to grant its Common Agricultural Policy subsidies and that farmers currently declare their crops manually. Land cover maps are also the starting point for change detection, which compares dates; the warning above about differencing two maps applies there.

Open Problems

• Generalization across regions. Benchmarks with spatially separated or cross-domain splits, such as LoveDA and reBEN, show lower scores than random splits suggest. Models that keep their accuracy in regions, seasons and sensors outside their training data remain a research goal.
• Reference data. Labels from maps carry the maps’ generalization and dates, and labels from people disagree on ambiguous classes. Dynamic World’s model reached roughly the agreement that annotators reached with each other, which leaves the reference data, not the model, as the limit on what can be measured (our reading of Brown et al., 2022).
• Evaluating continuously updated maps. Design-based accuracy assessment assumes a fixed map. How to state the accuracy of a map that changes with every image, as Dynamic World’s authors discuss, has no settled answer.
• Consistency over time. Maps produced with different algorithm versions change for reasons other than the ground, as WorldCover’s user manual warns for its 2020 and 2021 maps. Stable, versioned time series of maps are needed before map differences can be read as change.
• Ambiguous and mixed classes. Shrubland, grassland, wetlands and mixed agricultural mosaics remain the least accurate classes in both human labels and models, and better legends may matter as much as better classifiers (our reading of the WorldCover and Dynamic World results).