How to Download Land Cover: ESA WorldCover 10 m Tiles for Your Study Area

A land cover map says what covers the ground at each place: trees, crops, buildings, water. Making one yourself, as in the classification tutorials, takes labelled samples and care. ESA WorldCover is a free 10 m land cover map of the whole world that you can use as it is, or as a reference to compare your own map against.

WorldCover is made from Sentinel-1 radar and Sentinel-2 optical images. It exists for 2020, as version 1.0, and for 2021, as version 2.0. Each year is one map for the whole year, with eleven classes.

The classes

The file holds a colour for each code, so QGIS shows the map with these colours as it opens. A code without a class, such as 0, means no data.

Read a tile name

WorldCover is cut into tiles of 3 degrees by 3 degrees. A tile name looks like N30E075. N30 is the latitude of the south-west corner, 30 degrees north, and E075 is its longitude, 75 degrees east. The tile therefore covers 30 to 33 N and 75 to 78 E. Jagraon, at 30.79 N and 75.47 E, is in this tile.

To find the tile for your place, round the latitude and the longitude down to a multiple of 3. For a place at 30.79 N, the multiple of 3 below it is 30, and for 75.47 E it is 75.

Route A: download a tile with a link

The tiles are stored in a public bucket. A link downloads the file, with no account.

https://esa-worldcover.s3.eu-central-1.amazonaws.com/v200/2021/map/ESA_WorldCover_10m_2021_v200_N30E075_Map.tif

Change N30E075 to your tile. The file is about 100 MB, because it is a map of 3 by 3 degrees at 10 m. For the 2020 version, the path is v100/2020/map/ and the file name has 2020_v100. Open the file in QGIS and the colours appear.

Route B: Python

The script reads only your window of the tile. You need Python with numpy, pandas and rasterio.

import math
import numpy as np
import pandas as pd
import rasterio
from rasterio.windows import from_bounds

# Step 1: your area, as west, south, east, north in longitude and latitude
study_area = [75.40, 30.74, 75.55, 30.84]        # a box of about 14 km by 11 km near Jagraon, Punjab

# Step 2: the tile of the centre of the box. Round the latitude and longitude down to a multiple of 3.
def tile_name(lon, lat):
    lat_part = f"{'N' if lat >= 0 else 'S'}{int(math.floor(lat / 3) * 3):02d}"
    lon_part = f"{'E' if lon >= 0 else 'W'}{abs(int(math.floor(lon / 3) * 3)):03d}"
    return lat_part + lon_part

centre_lon = (study_area[0] + study_area[2]) / 2
centre_lat = (study_area[1] + study_area[3]) / 2
tile = tile_name(centre_lon, centre_lat)

# Step 3: the address of that tile for 2021. Nothing is downloaded yet.
address = f"https://esa-worldcover.s3.eu-central-1.amazonaws.com/v200/2021/map/ESA_WorldCover_10m_2021_v200_{tile}_Map.tif"
print("tile:", tile)

tile: N30E075

Check that the box lies inside the tile. A box that crosses the 3 degree line needs two tiles, which the guide on mosaics explains.

# Step 1: open the tile online and look at what it is
with rasterio.open(address) as src:
    print("tile size:", src.width, "x", src.height, "cells | cell size:", round(src.res[0], 7), "degrees | CRS: EPSG", src.crs.to_epsg())
    print("bounds:", [round(v, 1) for v in src.bounds], "| overviews:", src.overviews(1))

    # Step 2: read only the window of your box
    window = from_bounds(*study_area, src.transform).round_offsets().round_lengths()
    land_cover = src.read(1, window=window)

    # Step 3: the colour of each class and the credit line are stored in the file
    colours = src.colormap(1)
    credit = src.tags().get("copyright")
    licence = src.tags().get("license")

print("window:", land_cover.shape[1], "x", land_cover.shape[0], "cells")
print("copyright:", credit)
print("licence:", licence)

tile size: 36000 x 36000 cells | cell size: 8.33e-05 degrees | CRS: EPSG 4326
bounds: [75.0, 30.0, 78.0, 33.0] | overviews: [2, 4, 8, 16, 32, 64]
window: 1800 x 1200 cells
copyright: ESA WorldCover project 2021 / Contains modified Copernicus Sentinel data (2021) processed by ESA WorldCover consortium
licence: CC-BY 4.0 - https://creativecommons.org/licenses/by/4.0/

The tile has 36,000 cells on each side and overviews, which makes it a Cloud Optimized GeoTIFF, so only the window travels over the network. The credit line and the licence are part of the file. Copy them into your map and your report.

Count the classes

# Step 1: the names of the classes
names = {10: "Tree cover", 20: "Shrubland", 30: "Grassland", 40: "Cropland", 50: "Built-up", 60: "Bare or sparse vegetation",
         70: "Snow and ice", 80: "Permanent water", 90: "Herbaceous wetland", 95: "Mangroves", 100: "Moss and lichen"}

# Step 2: count the cells of every class
codes, counts = np.unique(land_cover, return_counts=True)

# Step 3: a cell is a rectangle of degrees. On the ground it is about 9.2 m north to south and 8.0 m east to west at this latitude.
cell_area_m2 = (0.0000833333 * 110574) * (0.0000833333 * 111320 * math.cos(math.radians(centre_lat)))
table = pd.DataFrame({"class": [names.get(int(c), "no data") for c in codes], "cells": counts})
table["share %"] = (100 * table["cells"] / counts.sum()).round(1)
table["area km2"] = (table["cells"] * cell_area_m2 / 1e6).round(1)
print(table.to_string(index=False))
print("cell area:", round(cell_area_m2, 1), "m2 | total:", round(table["area km2"].sum(), 1), "km2")

                    class   cells  share %  area km2
               Tree cover   71083      3.3       5.2
                Shrubland    7547      0.3       0.6
                Grassland   34861      1.6       2.6
                 Cropland 1825883     84.5     134.1
                 Built-up  212745      9.8      15.6
Bare or sparse vegetation    6810      0.3       0.5
          Permanent water    1058      0.0       0.1
       Herbaceous wetland      13      0.0       0.0
cell area: 73.4 m2 | total: 158.7 km2

Cropland covers most of the box, which fits a farming plain. Built-up land is the town and the villages. The tree cover is small: trees along roads and canals and around farmhouses. The areas are approximate, because the cell area was worked out for the centre latitude of the box. For an exact area, reproject the window to a UTM zone with nearest neighbour resampling, as in the reprojection guide, and count the cells.

Save the window with its colours

# Step 1: write the window to a GeoTIFF. The transform places it, and the colour table keeps the class colours.
with rasterio.open(address) as src:
    profile = src.profile.copy()
    profile.update(height=land_cover.shape[0], width=land_cover.shape[1], transform=src.window_transform(window), compress="deflate", tiled=False)
    profile.pop("blockxsize", None)
    profile.pop("blockysize", None)

with rasterio.open("worldcover_2021_jagraon.tif", "w", **profile) as dst:
    dst.write(land_cover, 1)
    dst.write_colormap(1, colours)

# Step 2: read it back to check
with rasterio.open("worldcover_2021_jagraon.tif") as check:
    print("saved:", check.width, "x", check.height, "cells | CRS EPSG", check.crs.to_epsg(), "| colours kept:", len(check.colormap(1)) > 0)

saved: 1800 x 1200 cells | CRS EPSG 4326 | colours kept: True

A land cover raster must be reprojected with nearest neighbour resampling, never bilinear, because averaging codes invents classes. The reprojection guide shows why.

Dynamic World, a land cover that updates

WorldCover is a map of a single year. Dynamic World, from Google and the World Resources Institute, is a different kind of product. It classifies every Sentinel-2 scene as it arrives, at 10 m, into nine classes: water, trees, grass, flooded vegetation, crops, shrub and scrub, built area, bare ground, and snow and ice. It has been running since mid 2015, so you can ask what a place looked like in a given month. It is published under a Creative Commons Attribution licence.

The catch is where it lives. Dynamic World is served through Google Earth Engine, and Earth Engine needs an account that you register for yourself, free for non-commercial use. Once you are signed in, the quickest way to use it is the Code Editor in the browser, at code.earthengine.google.com. It has a script panel in the middle, a Run button above it, and a Console tab on the right where print writes its answer.

Paste this script and press Run. It takes the Dynamic World scenes of one month over the Jagraon box used in the other guides, takes the most common class in each cell, and prints the share of each class.

var box = ee.Geometry.Rectangle([75.40, 30.74, 75.55, 30.84]);
var dwAll = ee.ImageCollection('GOOGLE/DYNAMICWORLD/V1').filterBounds(box);

function summary(start, end) {
  var dw = dwAll.filterDate(start, end);
  var label = dw.select('label').mode();          // the class that came up most often in each cell
  var h = ee.Dictionary(label.reduceRegion({reducer: ee.Reducer.frequencyHistogram(), geometry: box, scale: 10, maxPixels: 1e9}).get('label'));
  var total = ee.Number(h.values().reduce(ee.Reducer.sum()));
  return ee.Dictionary({scenes: dw.size(), cells: total.round(),
    shares: h.map(function(k, v) { return ee.Number(v).divide(total).multiply(1000).round().divide(10); })});
}

print('Feb 2023', summary('2023-02-01', '2023-03-01'));
print('Jul 2023', summary('2023-07-01', '2023-08-01'));
print('Nov 2023', summary('2023-11-01', '2023-12-01'));

I ran this script in the Code Editor, and the Console printed 1,858,688 cells for every month. The shares in percent, by class number, were:

Built area is steady, near 17 percent in all three months, as it should be for a town and its roads. Everything else moves with the seasons. In February, 82 percent of the cells are crops, the wheat of the winter. In July, crops drop to 4 percent and water rises to 74 percent. The most likely reason is that rice fields are flooded in July and Dynamic World labels flooded fields as water, together with the water of the monsoon rain, and a map of July 2023 can look like a flood map for that reason. In November, after the rice harvest, crops are back at 60 percent, and bare ground rises to 9 percent. The class snow and ice at 0.1 percent in November on the plains of Punjab is a mistake of the classifier, and a reminder to look at the result before you trust it.

The same query can be run from Python with the earthengine-api package, after a one-time earthengine authenticate on your computer. I ran this guide's numbers in the Code Editor, and I did not run the Python version.

Use WorldCover when you need one fixed, checked map of a year, and Dynamic World when you need to follow a change through time and can accept a per-scene result that is noisier. A Dynamic World image has a band called label, which holds the class number, and one band for the probability of each class.

What the map cannot tell you

• WorldCover is a global product. The project's own validation reports an overall accuracy of about three quarters, and the figure is lower for some classes and some regions. The classes of crops, grass and shrubs mix up easily. Check the validation report for your region, and check against places you know.
• The map is one picture of a year, and it does not show seasons. A field is labelled cropland whether it is green or bare in a given month.
• The 2020 and 2021 maps were made with different versions of the processing algorithm. The project advises against subtracting them to find change. Read the product user manual before you compare years.
• A cell of 10 m can hold a roof, a tree and a road. Narrow things, such as a lane or a canal, do not show up as classes of their own.
• A map from a global algorithm will not match a map that you made from labelled local samples. Compare them, and expect differences.

Common mistakes

• Using the wrong tile because the name uses the south-west corner and not the centre of the tile.
• Forgetting that a box across a 3 degree line needs two tiles.
• Counting cells and calling the number an area in square kilometres without working out the size of a cell.
• Resampling the classes with bilinear or cubic. Use nearest.
• Dropping the colour table, and getting a grey map of codes.
• Leaving out the credit line.
• Subtracting the 2020 and the 2021 versions to find change.

Key takeaways