--- title: "Land cover classification in Python: a Sentinel-2 random forest with scikit-learn" description: "Classify Sentinel-2 land cover in Python with a scikit-learn random forest: WorldCover labels, spatial cross-validation, F1 scores and a full map." url: https://www.technicalgis.com/tutorials/land-cover-classification-random-forest/ author: TechnicalGIS Team publisher: TechnicalGIS (https://www.technicalgis.com/about/) published: 2026-10-17 updated: 2026-10-17 language: en-GB --- # Land cover classification in Python: a Sentinel-2 random forest with scikit-learn > Classify Sentinel-2 land cover in Python with a scikit-learn random forest: WorldCover labels, spatial cross-validation, F1 scores and a full map. - Level: Intermediate - Tested with: Python 3.12 - Source code: [Download scripts](https://www.technicalgis.com/code/land-cover-classification-random-forest/) ## Key takeaways 1. Seasonal median composites (spring, summer and autumn) of ten Sentinel-2 bands plus four indices give 42 features per 20 m pixel. Spring water and red-edge indices were the most useful. 2. Sampling 2000 pure pixels per class from ESA WorldCover is a quick way to get training labels, but the model then learns WorldCover, including its errors. Hand-labelled reference data is better. 3. Hold out whole 2 km blocks, not random pixels. Here random and blocked cross-validation agreed within 0.2 percentage points (97.3% vs 97.1%) because the samples were sparse, but with dense samples the gap can be large. 4. Test-sample accuracy on pure pixels (97.7%) overstated whole-map agreement with WorldCover (87.0%) by more than 10 points, because the hard pixels on edges and in mixed areas were never tested. ## Before you start - Comfortable with NumPy and the basics of xarray; the NDVI time series tutorial introduces odc-stac and cloud masking in more detail - Some familiarity with scikit-learn's fit and predict pattern helps but is not essential - Python 3.10+ with scikit-learn, odc-stac, dask, xarray, pystac-client, planetary-computer, rasterio, scipy and matplotlib - An internet connection: imagery and labels are streamed from the Microsoft Planetary Computer and need no account **Short answer:** for supervised land cover classification in Python, build a cloud-free feature stack (for example seasonal Sentinel-2 median composites of the bands plus NDVI and other indices), collect labelled pixels for each class, train a **scikit-learn `RandomForestClassifier`**, and assess it on test data that is **spatially separate** from the training data, using a confusion matrix and per-class F1 scores. Then call `predict` on every pixel and save the result as a GeoTIFF. Assess accuracy with **blocked** cross-validation (`GroupKFold` on spatial blocks), not a random split, because neighbouring pixels are nearly identical. This tutorial classifies a 47 × 23 km area of **Sheffield and the Peak District** in northern England, which runs from heather and grass moorland with reservoirs in the west, through the city, to arable farmland and woodland in the east. It uses Sentinel-2 imagery for 2021 and training labels sampled from **ESA WorldCover 2021**. You will learn how to: 1. search a growing season of Sentinel-2 with STAC and build **spring, summer and autumn median composites** with odc-stac and dask; 2. turn them into a **feature stack** of bands plus NDVI, MNDWI, NDBI and NDRE; 3. load **ESA WorldCover** onto the same grid and draw a **stratified random sample** of training pixels; 4. train a **random forest** and test it on held-out **2 km spatial blocks**; 5. compare **random and spatially blocked cross-validation**; 6. read a **confusion matrix** and **per-class F1** scores; 7. rank features with **permutation importance**; 8. **predict a full land cover map** and see how it compares with the reference. All outputs below are real results from [`land_cover_random_forest.py`](https://www.technicalgis.com/code/land-cover-classification-random-forest/), run on 17 October 2026 with Python 3.12, scikit-learn 1.9.1, odc-stac 0.5.3 and xarray 2026.9.0. ## Step 1: search one growing season of Sentinel-2 A single image only shows one moment. Crops are bare in spring and green in summer, while grass is green almost all year, so several dates make classes much easier to separate. WorldCover 2021 describes 2021, so the script uses Sentinel-2 scenes from March to October of that year. The area spans two Sentinel-2 tiles (30UWE and 30UXE), so duplicates are removed per date and tile: `land_cover_random_forest.py` ```python BBOX = (-1.75, 53.25, -1.05, 53.45) # Peak District moors, Sheffield and the farmland to the east CRS, RES = "EPSG:32630", 20 # UTM zone 30N at 20 m: the native grid of the red-edge and SWIR bands YEAR = 2021 # match the year of the ESA WorldCover map used for the labels SEASONS = {"spring": [3, 4, 5], "summer": [6, 7, 8], "autumn": [9, 10]} BANDS = ["B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B11", "B12"] catalog = pystac_client.Client.open(STAC_URL, modifier=planetary_computer.sign_inplace) items = catalog.search(collections=["sentinel-2-l2a"], bbox=BBOX, datetime=f"{YEAR}-03-01/{YEAR}-10-31", query={"eo:cloud_cover": {"lt": 60}}).item_collection() best = {} # keep the newest processing baseline when a date and tile exist twice for it in items: key = (it.datetime.date(), it.properties["s2:mgrs_tile"]) if key not in best or it.properties["s2:processing_baseline"] > best[key].properties["s2:processing_baseline"]: best[key] = it items = sorted(best.values(), key=lambda i: i.datetime) ``` `Output — STAC search` ```text 68 Sentinel-2 L2A items (tiles ['30UWE', '30UXE']) with < 60% cloud, 07 Mar to 11 Oct 2021 ``` The classification works at **20 m**, the native resolution of the red-edge (B05–B07, B8A) and shortwave infrared (B11, B12) bands. The 10 m bands are averaged to 20 m. That quarters the number of pixels compared with 10 m and suits WorldCover, which is resampled to the same grid in Step 4. ## Step 2: seasonal median composites with odc-stac `odc.stac.load` builds a lazy datacube from the 68 items, mosaicking the two tiles on each date (`groupby="solar_day"`). The cloud mask uses the scene classification (SCL) band, as in the [NDVI time series tutorial](https://www.technicalgis.com/tutorials/ndvi-time-series-odc-stac/): keep classes 4–7, and grow cloud and shadow by two pixels (40 m). Then, for each season, the script takes the **median** of the clear observations: `land_cover_random_forest.py` ```python ds = odc.stac.load(items, bands=BANDS + ["SCL"], bbox=BBOX, crs=CRS, resolution=RES, chunks={"time": 1, "x": 2048, "y": 2048}, groupby="solar_day", dtype="uint16", nodata=0, resampling={"*": "average", "SCL": "nearest"}) ds = ds.assign_coords(time=ds.time.dt.floor("D")) gbox = ds.odc.geobox scl = ds["SCL"].compute() cloudy = np.isin(scl.values, [3, 8, 9, 10]) # shadow, cloud (medium, high), cirrus grown = ndimage.binary_dilation(cloudy, iterations=2, axes=(1, 2)) # grow by 40 m clear = scl.isin([4, 5, 6, 7]) & ~xr.DataArray(grown, coords=scl.coords, dims=scl.dims) frac = clear.mean(("y", "x")).to_series() layers = {} for season, months in SEASONS.items(): times = frac.index[(frac >= 0.1) & frac.index.month.isin(months)] sub = ds[BANDS].sel(time=times) refl = ((sub.astype("float32") - offset.sel(time=times)) / 10000).where(clear.sel(time=times)) med = refl.chunk({"time": -1, "y": 512, "x": 512}).median("time").compute() for b in BANDS: layers[f"{season}_{b}"] = med[b] comp = xr.Dataset(layers) ``` `Output — datacube and composites` ```text lazy cube: {'y': 1166, 'x': 2356, 'time': 42}, 2.5 GB if loaded in full dropped 0 dates whose tiles mix processing baselines before and after 04.00 SCL read in 3 s; 36 of 42 dates are at least 10% clear spring 18 dates; pixels with no clear view: 0.31% summer 9 dates; pixels with no clear view: 0.69% autumn 9 dates; pixels with no clear view: 0.78% composites built in 47 s ``` Two details matter: - **Rechunk before the median.** The cube is loaded with one chunk per date, which is right for reading. A median along time needs every date for a pixel in one chunk, so `.chunk({"time": -1, "y": 512, "x": 512})` swaps to small spatial tiles holding the full time series. - **The processing-baseline offset.** Scenes processed with baseline 04.00 or later store reflectance × 10,000 + 1000. The `offset` array holds 0 or 1000 per date, built from each item’s `s2:processing_baseline`. The script drops any date whose two tiles use different conventions, because they are mosaicked together; in 2021 there were none. Even in a wet English summer, the medians leave less than 1% of pixels with no clear view in any season. ![Two maps of a 47 by 23 km area. Top: a natural-colour summer composite with brown and green moorland and dark reservoirs in the west, the grey city of Sheffield in the centre and a patchwork of fields in the east. Bottom: ESA WorldCover for the same area, mostly yellow grassland in the west, red built-up in the centre, pink cropland in the east and green tree cover along valleys and in the south-east.](https://www.technicalgis.com/_astro/01-composite-and-worldcover.DKIccDhl_ZrPTao.webp) _Top: the summer 2021 median composite. Bottom: ESA WorldCover 2021 resampled to the same 20 m grid, which provides the training labels. WorldCover maps the moorland, including heather, as grassland._ ## Step 3: features: bands plus spectral indices Each season contributes its ten bands and four normalised-difference indices. A random forest could in principle learn these ratios from the bands, but giving them directly makes the trees shallower and the model more robust: `land_cover_random_forest.py` ```python def nd(a, b): return (a - b) / (a + b) features = {} for s in SEASONS: for b in BANDS: features[f"{s}_{b}"] = comp[f"{s}_{b}"] features[f"{s}_NDVI"] = nd(comp[f"{s}_B08"], comp[f"{s}_B04"]) # vegetation features[f"{s}_MNDWI"] = nd(comp[f"{s}_B03"], comp[f"{s}_B11"]) # open water features[f"{s}_NDBI"] = nd(comp[f"{s}_B11"], comp[f"{s}_B08"]) # built-up and bare surfaces features[f"{s}_NDRE"] = nd(comp[f"{s}_B8A"], comp[f"{s}_B05"]) # red-edge chlorophyll FEATURES = list(features) X_all = np.stack([features[f].values.astype("float32") for f in FEATURES], axis=-1) # (y, x, feature) H, W, F = X_all.shape ``` `Output — feature stack` ```text 1166 x 2356 pixels of 20 m, 42 features (3 seasons x (10 bands + 4 indices)) pixels with at least one missing feature: 1.56% ``` The result is an ordinary NumPy array of shape (rows, columns, features), which is what scikit-learn needs. At 2.7 million pixels × 42 features × 4 bytes it is about 460 MB, so for a much larger area, process it in tiles. ## Step 4: training labels from ESA WorldCover ESA WorldCover is also on the Planetary Computer. Passing the composite’s `geobox` to `odc.stac.load` puts it on exactly the same 20 m grid. `resampling="mode"` takes the most common 10 m class in each 20 m pixel, which is the right choice for categorical data: `land_cover_random_forest.py` ```python wc_items = catalog.search(collections=["esa-worldcover"], bbox=BBOX, query={"esa_worldcover:product_version": {"eq": "2.0.0"}}).item_collection() wc = odc.stac.load(wc_items, bands=["map"], geobox=gbox, resampling="mode").squeeze("time")["map"].values share = pd.Series(wc.ravel()).value_counts(normalize=True).sort_index() MIN_SHARE = 0.005 # drop classes covering less than 0.5% of the area USE = [c for c in CLASSES if share.get(c, 0) >= MIN_SHARE] ``` `Output — ESA WorldCover 2021 class shares (%)` ```text 10 Tree cover 27.21 20 Shrubland 0.00 30 Grassland 36.99 40 Cropland 20.87 50 Built-up 14.05 60 Bare / sparse 0.14 80 Water 0.72 90 Herbaceous wetland 0.02 classes used: ['Tree cover', 'Grassland', 'Cropland', 'Built-up', 'Water']; dropped: ['Shrubland', 'Bare / sparse', 'Herbaceous wetland'] ``` Version 2.0.0 is the 2021 map. Three classes cover less than 0.5% of the area and are dropped, because a few hundred pixels, possibly misclassified, cannot train a class reliably. Pixels of those classes are still predicted as one of the five classes used. ## Step 5: a stratified random sample of pure pixels With five classes of very different sizes (water is 0.7%, grassland 37%), a simple random sample would contain almost no water. A **stratified** random sample draws the same number from each class. The script also keeps only “pure” pixels whose 3 × 3 neighbourhood is all the same class, which removes mixed pixels on boundaries and many of the reference map’s edge errors: `land_cover_random_forest.py` ```python N_PER_CLASS = 2000 BLOCK = 2000 # 2 km x 2 km spatial blocks for the train/test split rng = np.random.default_rng(SEED) valid = ~np.isnan(X_all).any(-1) rows, cols, ys = [], [], [] for c in USE: # Keep pixels whose 3 x 3 neighbourhood is all the same class: fewer mixed pixels and edge errors. pure = ndimage.binary_erosion(wc == c, structure=np.ones((3, 3))) & valid r, cc = np.nonzero(pure) pick = rng.choice(len(r), size=min(N_PER_CLASS, len(r)), replace=False) rows.append(r[pick]); cols.append(cc[pick]); ys.append(np.full(len(pick), c)) rows, cols, y = np.concatenate(rows), np.concatenate(cols), np.concatenate(ys) X = X_all[rows, cols] block_px = BLOCK // RES block_id = (rows // block_px) * (W // block_px + 1) + cols // block_px ``` `Output — stratified random sample` ```text Tree cover 2000 Grassland 2000 Cropland 2000 Built-up 2000 Water 2000 10000 samples in 288 blocks of 2 km ``` Each sample also records which 2 km × 2 km block it falls in. Those blocks drive the train/test split in the next step. ## Step 6: train a random forest and test on spatial blocks Pixels close together look alike. If you split samples randomly, most test pixels have a training pixel a few tens of metres away, and the model can “recognise” them rather than generalise. That makes accuracy look better than it will be on new areas. The fix is to **split by space**: put every block into one of five folds, and hold out fold 0 (about 20% of blocks) entirely for testing: `land_cover_random_forest.py` ```python blocks = np.unique(block_id) fold_of_block = dict(zip(blocks, rng.permutation(len(blocks)) % 5)) fold = np.array([fold_of_block[b] for b in block_id]) test = fold == 0 # one fifth of the blocks are held out entirely rf = RandomForestClassifier(n_estimators=300, max_features="sqrt", min_samples_leaf=1, n_jobs=-1, random_state=SEED, oob_score=True) rf.fit(X[~test], y[~test]) pred = rf.predict(X[test]) print(classification_report(y[test], pred, labels=USE, target_names=names, digits=3)) ``` `Output — spatial hold-out test (20% of blocks)` ```text trained on 8050 samples in 0.6 s; tested on 1950 samples overall accuracy 0.977 | macro F1 0.977 | out-of-bag accuracy 0.972 precision recall f1-score support Tree cover 0.959 0.961 0.960 361 Grassland 0.986 0.961 0.973 432 Cropland 0.974 0.974 0.974 387 Built-up 0.969 0.991 0.980 447 Water 1.000 1.000 1.000 323 accuracy 0.977 1950 macro avg 0.978 0.977 0.977 1950 weighted avg 0.977 0.977 0.977 1950 area-weighted accuracy (class recall weighted by WorldCover area share): 0.968 ``` Random forests need little tuning. 300 trees and `max_features="sqrt"` (about 6 of the 42 features tried at each split) are standard defaults; more trees rarely help much. The **out-of-bag** score (97.2%) is a free estimate from the samples each tree did not see, but it uses random samples, so it has the same spatial leakage problem as a random split. Because the sample is balanced, overall accuracy treats each class equally. The **area-weighted accuracy** (96.8%) weights each class’s recall by its share of the map, which is closer to what a map user experiences. In remote sensing terms, **precision is user’s accuracy** and **recall is producer’s accuracy**. ![A five by five confusion matrix with WorldCover labels as rows and random forest predictions as columns. The diagonal is dark blue at 96 to 100%. The largest errors are 12 grassland pixels predicted as tree cover and 9 tree cover pixels predicted as built-up.](https://www.technicalgis.com/_astro/03-confusion-matrix.BwxIOxKA_ZOPRzI.webp) _Confusion matrix on the 1950 held-out test samples, with row percentages and counts. The largest confusions are grassland predicted as tree cover (12) and tree cover predicted as built-up (9), which is typical of leafy suburbs and scrubby pasture._ ### Random vs spatially blocked cross-validation One hold-out fold is a small test. Five-fold cross-validation uses every sample once for testing. The script runs it twice, with `KFold` (random) and with `GroupKFold` using the spatial folds as groups: `land_cover_random_forest.py` ```python for name, splitter, groups in [("random", KFold(5, shuffle=True, random_state=SEED), None), ("spatial blocks", GroupKFold(5), fold)]: p = cross_val_predict(RandomForestClassifier(n_estimators=300, n_jobs=-1, random_state=SEED), X, y, cv=splitter, groups=groups) ``` `Output — random vs spatially blocked 5-fold cross-validation` ```text random accuracy 0.973 | macro F1 0.973 | Tree 0.96 Grassland 0.96 Cropland 0.97 Built-up 0.98 Water 1.00 spatial blocks accuracy 0.971 | macro F1 0.971 | Tree 0.96 Grassland 0.96 Cropland 0.97 Built-up 0.97 Water 1.00 random CV overstates accuracy by 0.2 percentage points ``` Here the gap is tiny. That is an honest result, and it has a reason: 10,000 samples spread over about 2.7 million pixels are, on average, hundreds of metres apart, so few test pixels have a near-twin in the training set. The gap grows when samples are dense, for example when every pixel inside hand-drawn training polygons is used, which is the usual way of working in desktop GIS. Then random splits put neighbouring pixels from the same polygon on both sides and can overstate accuracy by 10 points or more. Blocked cross-validation costs nothing extra, so use it by default. ## Step 7: how good is the whole map? The tests above all use pure pixels, which are the easy ones. The real question is how the model does on **every** pixel, including field edges, hedgerows, gardens and mixed pixels. Since WorldCover covers the whole area, the script can compare every pixel: ![A grouped bar chart of F1 score for five classes. For each class, random and spatial cross-validation both give about 0.96 to 1.0, while the full-map comparison gives 0.82 for tree cover, 0.89 for grassland, 0.90 for cropland, 0.85 for built-up and 0.87 for water.](https://www.technicalgis.com/_astro/04-f1-three-estimates.Zr7TEjip_nfE6o.webp) _Per-class F1 against WorldCover from random cross-validation, spatially blocked cross-validation and every pixel of the final map. Overall agreement drops from 97% on pure-pixel samples to 87% over the whole map._ `Output — full map prediction` ```text predicted 2,747,096 pixels in 6 s agreement with WorldCover over the classes used: 87.0%; median top-class probability 0.96 per-class F1 against WorldCover over every pixel: Tree 0.82 Grassland 0.89 Cropland 0.90 Built-up 0.85 Water 0.87 ``` That is a drop of about 10 points, and it is the most important lesson of this tutorial. **A test set drawn the same way as the training set inherits its biases.** Excluding edge pixels made both training and testing cleaner, but the map still has to label those edges. A proper accuracy assessment uses an independent, probability-based sample of the whole map, labelled by people, as described by Olofsson et al. (2014). ## Step 8: feature importance Impurity-based importance (`rf.feature_importances_`) is quick but biased towards features with many distinct values, and it spreads credit across correlated features. **Permutation importance** on the held-out test blocks is more honest: shuffle one feature and see how much macro F1 drops: `land_cover_random_forest.py` ```python pi = permutation_importance(rf, X[test], y[test], n_repeats=5, random_state=SEED, n_jobs=-1, scoring="f1_macro") imp = pd.DataFrame({"permutation": pi.importances_mean, "perm_sd": pi.importances_std, "impurity": rf.feature_importances_}, index=FEATURES).sort_values("permutation", ascending=False) ``` `Output — feature importance (top 6 of 42)` ```text permutation perm_sd impurity spring_MNDWI 0.0121 0.0007 0.0423 autumn_NDRE 0.0039 0.0016 0.0661 spring_B11 0.0034 0.0013 0.0536 spring_B8A 0.0030 0.0016 0.0624 spring_B08 0.0024 0.0013 0.0436 summer_NDRE 0.0020 0.0010 0.0513 impurity importance by season: {'autumn': 0.292, 'spring': 0.39, 'summer': 0.318} ``` ![A horizontal bar chart of the 15 most important features by permutation importance, coloured by season. Spring MNDWI is by far the largest at about 0.012; the rest are below 0.004, led by autumn NDRE, spring B11, spring B8A and spring B08.](https://www.technicalgis.com/_astro/05-feature-importance.Dy7bkZCu_Z1wesJE.webp) _Permutation importance of the top 15 features, coloured by season, with error bars of one standard deviation over five shuffles. Spring features dominate, and the indices (MNDWI, NDRE) rank higher than most raw bands._ The drops are all small, because the 42 features are highly correlated: shuffle one, and the forest uses a near-copy from another season. Treat the ranking as a hint rather than a measurement. Still, it makes sense: in **spring**, arable fields are bare or sparsely covered, so shortwave infrared (B11) and the water index separate them from grass. The **red-edge** index NDRE separates crops, grass and trees by chlorophyll content, and spring features as a group carry the most weight (0.39 of impurity importance). If you need fewer features, drop correlated ones and check that blocked cross-validation does not get worse. ## Predicting the full map The final model is trained on all 10,000 samples and applied to every pixel in chunks of 200,000 to limit memory. Since scikit-learn 1.4, random forests handle `NaN` in features, so the 1.6% of pixels with a missing season still get a prediction. The result is written as a GeoTIFF with the WorldCover colour table, ready for QGIS: `land_cover_random_forest.py` ```python final = RandomForestClassifier(n_estimators=300, max_features="sqrt", n_jobs=-1, random_state=SEED).fit(X, y) flat = X_all.reshape(-1, F) pred_map = np.zeros(H * W, dtype="uint8") for i in range(0, H * W, 200_000): # predict in chunks to keep memory low proba = final.predict_proba(flat[i:i + 200_000]) # scikit-learn >= 1.4 forests accept missing values pred_map[i:i + 200_000] = final.classes_[proba.argmax(1)] pred_map = pred_map.reshape(H, W) with rasterio.open(OUT / "land_cover_rf_2021.tif", "w", driver="GTiff", height=H, width=W, count=1, dtype="uint8", crs=CRS, transform=gbox.affine, compress="deflate", nodata=0) as dst: dst.write(pred_map, 1) ``` `Output — area by class (hectares)` ```text WorldCover_ha random_forest_ha Tree cover 29897 26388 Grassland 40651 40093 Cropland 22928 25171 Built-up 15437 17527 Water 796 705 ``` ![Top: the random forest land cover map, very similar to WorldCover, with yellow grassland in the west, red Sheffield in the centre, pink cropland in the east and green woodland. Bottom: a greyscale image with red pixels where the two maps disagree, scattered thinly along field edges, roads, rivers and through the urban area, with one denser patch on the moors in the west.](https://www.technicalgis.com/_astro/06-predicted-map.BePBefZq_Z2mjDVu.webp) _Top: the random forest land cover map for 2021. Bottom: pixels where it disagrees with WorldCover (red). Disagreement follows boundaries, roads and the fine-grained urban fabric, where 20 m pixels are mixed._ The map looks very like WorldCover, which is what it was trained to do. It has less tree cover (−12%) and more built-up land (+14%) and cropland (+10%). Many of the differences are in the suburbs, where a 20 m pixel of houses and gardens with trees is genuinely mixed, so neither map is clearly right. To learn which is closer to the truth you would need independent reference data. ## Pitfalls and limitations - **The labels are a map, not the truth.** WorldCover’s global overall accuracy is about 77%, and its errors become training data. It also maps most heather moorland in the Peak District as grassland, so this model cannot learn a heath class. Hand-labelled points or polygons are better, and essential for any published map. - **Pure-pixel tests are optimistic.** As Step 7 shows, accuracy on samples drawn away from edges (97.7%) is much higher than agreement over every pixel (87.0%). Report accuracy from an independent sample of the whole map. - **Blocked cross-validation is not a guarantee.** 2 km blocks reduce leakage between nearby pixels, but every block is still in the same landscape. The model may do worse in, say, East Anglia or Scotland. - **Rare classes were dropped.** Shrubland, bare ground and wetland together cover about 0.2% of this area. Their pixels are forced into one of the five classes. - **Balanced training changes the priors.** Training with 2000 samples per class tells the model that water is as common as grassland. That favours small classes in ambiguous pixels. You can rebalance, or use `class_weight`, if area estimates matter. - **Clouds and composites.** Median composites hide most cloud, but SCL misses some haze and shadow, and pixels with few clear dates are noisier. - **One year.** The model is trained for 2021 imagery; applying it to another year needs new composites and, ideally, a check against new reference data. ## Frequently asked questions ### How do I do supervised classification in Python? Build a feature array where each pixel is a row (for example by stacking Sentinel-2 bands and indices with NumPy or xarray), collect labelled pixels, fit a scikit-learn classifier such as `RandomForestClassifier`, check accuracy on spatially separate test data, then call `predict` on every pixel and reshape the result back into an image. Rasterio writes it as a GeoTIFF. This tutorial’s script does every step on real data. ### Why use a random forest for Sentinel-2 classification? Random forests are accurate, fast, need almost no tuning, cope with correlated features and different value ranges without scaling, give class probabilities and handle missing values. Gradient boosting (LightGBM, XGBoost) is often slightly more accurate but needs more tuning. Deep learning (for example U-Nets) can use spatial context but needs far more labelled data. ### How many training samples do I need? Hundreds to a few thousand per class is typical for a random forest. More important than the number is that samples cover the variety within each class (for example, different crops and different kinds of woodland) and the whole area. Thousands of samples from a few polygons are worth less than a few hundred spread widely. ### What is spatial cross-validation and why does it matter? It groups samples by location, for example into 2 km blocks, and holds out whole blocks. A random split puts neighbouring, nearly identical pixels in both training and test sets, which can make accuracy look much better than it really is. In scikit-learn, use `GroupKFold` with block IDs as groups. In R, `CAST` and `blockCV` do the same. ### Can I use ESA WorldCover or Dynamic World as training data? For learning and prototyping, yes, as here. For a real product, it means copying another map’s errors. If you use one, keep only pixels where two maps agree (for example WorldCover and Esri’s IO-LULC), use pure pixels away from edges, and check the result against independent, hand-labelled points. ### How would I do this in Google Earth Engine? The workflow maps across directly: filter `COPERNICUS/S2_SR_HARMONIZED` and mask clouds, make seasonal `median()` composites, add indices, `stratifiedSample` the `ESA/WorldCover/v200` map, train `ee.Classifier.smileRandomForest(300)`, `classify` the image, and use `errorMatrix` for accuracy. Earth Engine runs the computation on Google’s servers, which is convenient for very large areas. We use open STAC data here because Earth Engine needs a Google account and a registered Cloud project, and because the Python version lets you use any scikit-learn model and validation method, including blocked cross-validation, which Earth Engine does not provide out of the box. ## Next steps - Build crop-specific features from a time series instead of seasonal medians, as in [NDVI time series in Python with odc-stac](https://www.technicalgis.com/tutorials/ndvi-time-series-odc-stac/). - Add terrain: elevation and slope help separate moorland, woodland and farmland. See [raster terrain analysis with rasterio](https://www.technicalgis.com/tutorials/raster-terrain-analysis-rasterio/). - Add Sentinel-1 radar features, which see through cloud and are sensitive to structure, as in [Sentinel-1 SAR flood mapping in Python](https://www.technicalgis.com/tutorials/sentinel-1-sar-flood-mapping/). - Turn the map into polygons and summarise land cover by area with [fast vector geoprocessing](https://www.technicalgis.com/tutorials/fast-vector-geoprocessing/), or share it on [an interactive web map with Folium](https://www.technicalgis.com/tutorials/interactive-web-maps-folium/). - Start with the STAC basics in [Measuring Lake Mead’s shrinking shoreline with Sentinel-2 and STAC](https://www.technicalgis.com/tutorials/sentinel-2-stac-lake-mead/). - Segment individual objects such as buildings and fields, rather than classifying pixels, with [Segment Anything on satellite and aerial imagery](https://www.technicalgis.com/tutorials/segment-anything-satellite-imagery/). ## Data and licences - **Sentinel-2 L2A**: contains modified Copernicus Sentinel data 2021, processed by ESA. Accessed as Cloud Optimised GeoTIFFs from the [Microsoft Planetary Computer](https://planetarycomputer.microsoft.com/dataset/sentinel-2-l2a) on 17 October 2026. - **ESA WorldCover 2021 (v200)**: © ESA WorldCover project 2021 / Contains modified Copernicus Sentinel data (2021) processed by the ESA WorldCover consortium. Licensed under [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/). Accessed from the [Planetary Computer](https://planetarycomputer.microsoft.com/dataset/esa-worldcover) on 17 October 2026. Zanaga, D. et al. (2022), ESA WorldCover 10 m 2021 v200, doi:10.5281/zenodo.7254221. - References: Breiman, L. (2001), “Random forests”, _Machine Learning_ 45; Olofsson, P. et al. (2014), “Good practices for estimating area and assessing accuracy of land change”, _Remote Sensing of Environment_ 148; Roberts, D. R. et al. (2017), “Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure”, _Ecography_ 40(8). - Code: [`land_cover_random_forest.py`](https://www.technicalgis.com/code/land-cover-classification-random-forest/), tested with Python 3.12.0, scikit-learn 1.9.1, odc-stac 0.5.3, xarray 2026.9.0, dask 2026.8.0 and pystac-client 0.9.0. ## Key terms - **Supervised classification**: Assigning a class to every pixel with a model trained on examples whose class is known (training labels), as opposed to unsupervised clustering. - **Random forest**: An ensemble of decision trees, each trained on a bootstrap sample of the data and a random subset of features at each split. The forest votes on the class. - **Stratified random sample**: A random sample drawn separately within each class (stratum), so rare classes get enough samples. - **Spatial autocorrelation**: The tendency of nearby pixels to be similar. Random train/test splits leak information when test pixels sit next to training pixels. - **Blocked cross-validation**: Cross-validation in which samples are grouped into spatial blocks and whole blocks are held out, so test data is spatially separate from training data. - **F1 score**: The harmonic mean of precision (user's accuracy) and recall (producer's accuracy) for one class. Macro F1 is the unweighted mean over classes. - **Permutation importance**: The drop in a model's score when the values of one feature are randomly shuffled on test data, which breaks the link between that feature and the label. ## Data for this tutorial - [Microsoft Planetary Computer](https://www.technicalgis.com/data-sources/#src-microsoft-planetary-computer) · Sentinel-2 L2A and ESA WorldCover via STAC - [ESA WorldCover](https://www.technicalgis.com/data-sources/#src-esa-worldcover) · 10 m global land cover for 2020 and 2021 - [Copernicus Data Space Ecosystem](https://www.technicalgis.com/data-sources/#src-copernicus-data-space-ecosystem) · The official source of Sentinel-2 data Licences and formats for each source are listed in the [GIS data sources catalogue](https://www.technicalgis.com/data-sources/). --- Source: https://www.technicalgis.com/tutorials/land-cover-classification-random-forest/