Land Cover Products Explained

Problem statement

A global land cover map looks like a photograph with the colours tidied up. It is a model's answer to a classification question, and the question differs between products: which class dominates this 10 m pixel in this year, according to training labels drawn by a particular team, from a particular satellite, with a particular class list. ESA WorldCover, Impact Observatory's annual maps, Google's Dynamic World and the Copernicus Global Land Service all answer versions of that question, and they disagree most where the answer is most interesting — edges, mosaics, regrowth and change.

Choosing a product is choosing a set of definitions and a set of errors. Knowing how they are built is how to choose.

Quick answer

import numpy as np, rasterio, pystac_client, planetary_computer
from rasterio.windows import from_bounds
from rasterio.warp import transform_bounds

W, S, E, N = -62.95, -10.75, -62.80, -10.60
cat = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1")

def read_class(collection, dt):
    item = next(cat.search(collections=[collection], bbox=[W, S, E, N], datetime=dt).items())
    href = planetary_computer.sign(item.assets["map" if "map" in item.assets else "data"].href)
    with rasterio.open(href) as src:
        win = from_bounds(*transform_bounds("EPSG:4326", src.crs, W, S, E, N), src.transform)
        return src.read(1, window=win), item.id

wc, wc_id = read_class("esa-worldcover", "2021-01-01/2021-12-31")
io, io_id = read_class("io-lulc-annual-v02", "2021-01-01/2021-12-31")
WC = {10: "tree cover", 20: "shrubland", 30: "grassland", 40: "cropland", 50: "built-up", 80: "water"}
IO = {1: "water", 2: "trees", 4: "flooded veg.", 5: "crops", 7: "built area", 11: "rangeland"}
for name, arr, labels in (("WorldCover 2021", wc, WC), ("Impact Observatory 2021", io, IO)):
    vals, counts = np.unique(arr, return_counts=True)
    top = sorted(zip(counts, vals), reverse=True)[:4]
    print(f"{name}: " + ", ".join(f"{labels.get(int(v), v)} {c / arr.size:.0%}" for c, v in top))

WorldCover 2021: grassland 60%, tree cover 40%, cropland 0%, built-up 0%
Impact Observatory 2021: rangeland 70%, trees 30%, crops 0%, water 0%

The same 16 × 16 km of Rondônia in the same year: both products agree it is mostly trees and grassland/rangeland, and disagree on how much of each — because their "grassland" and "rangeland" classes are not the same class.

Class definitions, year, resolution, stated accuracy

1. Read the class definitions

"Tree cover" in WorldCover means at least 10% tree canopy; "trees" in Impact Observatory is a learned label for areas dominated by trees. "Rangeland" includes grassland and shrubland. Classes with the same name are not the same class.

2. Check the year and version

WorldCover's 2020 and 2021 maps used different algorithms: differences between them are partly method. Annual products are more consistent across years but still change with model versions.

3. Match the resolution to the features

At 10 m, field boundaries, riparian strips and small clearings are resolved; at 30 m or coarser they blend into mixed pixels. The product's resolution sets the smallest thing it can say anything about.

4. Look at the stated accuracy by class

Overall accuracies of 75–85% hide much lower accuracy for rare or transitional classes. Check the class-level figures for the classes you care about.

5. Do not measure change by differencing two products

Differences between two maps include every classification error in both. Change needs a product designed for change, or a time series from one consistent method.

Crosswalking classes and comparing resolutions

The data used on this page

The Rondônia window used throughout the forest guides (63.0°–62.5°W, 10.8°–10.3°S): Hansen Global Forest Change v1.12 — tree cover in 2000, the year of loss, and the data mask — read from Google's storage into a three-band rondonia_gfc.tif, as in the Hansen reading guide.

import os

import rasterio
from rasterio.windows import from_bounds

W, S, E, N = -63.0, -10.8, -62.5, -10.3
V = "GFC-2024-v1.12"
if not os.path.exists("rondonia_gfc.tif"):
    bands = []
    for name in ("treecover2000", "lossyear", "datamask"):
        url = f"https://storage.googleapis.com/earthenginepartners-hansen/{V}/Hansen_{V}_{name}_10S_070W.tif"
        with rasterio.open(url) as src:
            window = from_bounds(W, S, E, N, src.transform)
            bands.append(src.read(1, window=window))
            profile = src.profile | {"count": 3, "width": bands[0].shape[1], "height": bands[0].shape[0],
                                     "transform": src.window_transform(window), "compress": "deflate",
                                     "driver": "GTiff", "tiled": False}
    profile.pop("blockxsize", None)
    profile.pop("blockysize", None)
    with rasterio.open("rondonia_gfc.tif", "w", **profile) as dst:
        for i, band in enumerate(bands, start=1):
            dst.write(band, i)
with rasterio.open("rondonia_gfc.tif") as src:
    print(src.count, src.shape, src.crs, src.res)

3 (2000, 2000) EPSG:4326 (0.00025, 0.00025)

Example 1 — the class crosswalk in practice

from rasterio.warp import reproject, Resampling
io_on_wc = np.zeros(wc.shape, dtype=io.dtype)
with rasterio.open(planetary_computer.sign(next(cat.search(
        collections=["io-lulc-annual-v02"], bbox=[W, S, E, N],
        datetime="2021-01-01/2021-12-31").items()).assets["data"].href)) as src_io, \
     rasterio.open(planetary_computer.sign(next(cat.search(
        collections=["esa-worldcover"], bbox=[W, S, E, N],
        datetime="2021-01-01/2021-12-31").items()).assets["map"].href)) as src_wc:
    win = from_bounds(*transform_bounds("EPSG:4326", src_wc.crs, W, S, E, N), src_wc.transform)
    reproject(rasterio.band(src_io, 1), io_on_wc, dst_transform=src_wc.window_transform(win),
              dst_crs=src_wc.crs, resampling=Resampling.nearest)
for wc_val in (10, 30):
    m = wc == wc_val
    vals, counts = np.unique(io_on_wc[m], return_counts=True)
    top = sorted(zip(counts, vals), reverse=True)[:3]
    print(f"WorldCover '{WC[wc_val]}' is Impact Observatory: " +
          ", ".join(f"{IO.get(int(v), v)} {c / m.sum():.0%}" for c, v in top))

WorldCover 'tree cover' is Impact Observatory: trees 73%, rangeland 27%, water 0%
WorldCover 'grassland' is Impact Observatory: rangeland 99%, trees 1%, crops 0%

WorldCover's grassland is almost all rangeland in Impact Observatory (99%), but more than a quarter of its tree cover is rangeland too: WorldCover calls a pixel tree cover from 10% canopy, Impact Observatory needs trees to dominate. A crosswalk table between products is always an approximation.

Example 2 — what 10 m sees that 30 m does not

coarse = wc[:wc.shape[0] // 3 * 3, :wc.shape[1] // 3 * 3].reshape(
    wc.shape[0] // 3, 3, wc.shape[1] // 3, 3)
mixed = np.array([[len(np.unique(coarse[i, :, j, :])) > 1 for j in range(coarse.shape[2])]
                  for i in range(coarse.shape[0])])
print(f"30 m blocks containing more than one 10 m WorldCover class: {mixed.mean():.0%}")

30 m blocks containing more than one 10 m WorldCover class: 12%

One 30 m cell in eight mixes classes at 10 m — the fraction of this landscape where a coarser product must pick one label for several land covers.

Example 3 — class shares against Hansen

with rasterio.open("rondonia_gfc.tif") as src:
    win = from_bounds(W, S, E, N, src.transform)
    tc = src.read(1, window=win)
    ly = src.read(2, window=win)
hansen_forest_2021 = ((tc > 30) & ((ly == 0) | (ly > 21))).mean()
print(f"Hansen forest 2021 (>30% canopy): {hansen_forest_2021:.0%}")
print(f"WorldCover tree cover 2021:       {(wc == 10).mean():.0%}")
print(f"Impact Observatory trees 2021:    {(io == 2).mean():.0%}")

Hansen forest 2021 (>30% canopy): 34%
WorldCover tree cover 2021:       40%
Impact Observatory trees 2021:    30%

Three honest estimates of the forest share of the same land in the same year. Which is "right" depends on the forest definition the question uses.

Why products disagree, and land cover is not land use

Why products disagree even with the same satellite

Classification depends on training labels, class definitions, the model and post-processing. Two teams using the same Sentinel-2 imagery produce different maps because they disagree about what the classes are and taught their models with different examples.

Why land cover is not land use

A pasture and a natural grassland look alike from space; so do a plantation and a forest. Land cover products describe the surface; land use — what people do with it — needs other information. Policy questions are usually about use.

Why differencing maps overstates change

If each map is 85% accurate, two independent maps disagree on a substantial share of pixels even where nothing changed. Differencing them reports that disagreement as change. Change products model time directly to avoid this.

Probabilities, tiles and licences

• Probabilities: Dynamic World provides per-class probabilities, which are more informative than the top label in mixed pixels.
• Tiles and projections: WorldCover uses 3° tiles in geographic coordinates; Impact Observatory uses UTM tiles. Reproject with nearest or mode, never bilinear.
• Accuracy reports are usually global; regional accuracy can be much lower.
• Water and wetland classes vary with season and year; comparisons across years pick up hydrology.
• Licences differ: WorldCover CC-BY 4.0, Impact Observatory CC-BY 4.0; check each.

Internal links

• How to compare land cover products for forest — the comparison in full
• Forest definitions explained — definitions behind "tree cover"
• How to classify land cover in Python — making your own map
• Global Forest Change explained — a continuous alternative
• Raster resampling explained — categorical resampling
• How to download satellite imagery from STAC — the Planetary Computer catalogue

FAQ

What is the best global land cover map?

There is no single best. Choose by year, resolution and class definitions, and check class-level accuracy for your region.

Why do WorldCover and other maps disagree?

Different class definitions, training data, models and years.

Can I detect change by comparing two land cover maps?

Not reliably; classification errors appear as change. Use products or methods designed for change.

What is the difference between land cover and land use?

Land cover is what is on the surface; land use is what people use it for.

How should I resample a land cover map?

With nearest neighbour or mode, never bilinear.

Are 10 m maps always better?

They resolve more detail but can be noisier; the best resolution depends on the features and question.