TerraShift

Uncertainty-aware land-cover segmentation and multi-year change analytics for Sentinel-2.

TerraShift is an end-to-end deep learning system that turns raw Sentinel-2 satellite imagery into georeferenced land-cover maps, quantifies how confident the model is at every single pixel, and measures how land cover changes across years. It covers the full lifecycle: satellite data engineering, model research, model fusion, uncertainty quantification, production inference over a real area of interest, and cartographic reporting.

The reference deployment maps the land cover of Vilnius, Lithuania, for 2015, 2020, and 2025, and quantifies urban expansion over that decade.

Final ensemble land-cover prediction and uncertainty maps for Vilnius, 2025.

At a glance

	Summary
Task	7-class semantic segmentation of Sentinel-2 imagery at 10 m ground resolution
Training data	LandCoverNet Europe 2018, reprocessed into cloud-free summer median composites
Models	1 custom Residual U-Net + 3 ImageNet-pretrained architectures + 1 self-supervised Sentinel-2 encoder
Experiments	27 controlled training runs on one fixed 70/15/15 split
Fusion	Weighted soft-probability ensemble of 5 heterogeneous models
Headline result	Pixel Accuracy 0.7938 · mIoU 0.5816 · mDice 0.7064 (+0.0233 mIoU over the best single model)
Uncertainty	Per-pixel predictive entropy, expected entropy, and mutual information, exported as GeoTIFF layers
Applied output	Georeferenced land-cover maps and class-share statistics for Vilnius in 2015, 2020, and 2025
Stack	Python, PyTorch, rasterio, NumPy, Matplotlib, `segmentation_models_pytorch`, `torchgeo`

Skills and capabilities demonstrated

This section is for readers evaluating the engineering behind the project rather than the land-cover result itself. Every claim below maps to code in this repository.

1. Remote-sensing data engineering

Satellite data does not arrive as clean tensors. A large share of the work here is turning a noisy, irregularly sampled, cloud-contaminated multi-spectral archive into a training-ready dataset.

Built a reproducible ingestion pipeline that walks a tiled raw Sentinel-2 archive, resolves inconsistent file naming (files with and without extensions, B8A-style band aliases, alternative label-file patterns), and fails loudly on shape or geometry mismatches instead of silently producing corrupt tensors.
Implemented physically motivated cloud, shadow, cirrus, and snow removal using the Sentinel-2 Scene Classification Layer, followed by a per-pixel temporal median composite across the summer season. Masked pixels become NaN and are excluded from the median rather than biasing it.
Solved a real compatibility problem end to end: LandCoverNet ships no B10 cirrus band, but the Sentinel-2 self-supervised pretrained weights expect a 13-channel input in a fixed band order. The pipeline injects a zero-filled synthetic B10 at the correct spectral position so pretrained weights load without any reindexing.
Preserved geospatial correctness throughout. CRS, affine transform, and nodata metadata are carried from the source rasters into every derived product, so predictions land back on the map exactly where they belong.
Produced three coordinated dataset variants (4-, 10-, and 13-band), a metadata manifest with per-chip provenance, and dedicated inspection and integrity-check tooling.

Code: prepare_landcovernet_median_dataset_multiband.py, check_processed_dataset.py, inspect_processed_dataset.py

2. Deep learning architecture design

Implemented a Residual U-Net from scratch in PyTorch: residual double-convolution blocks with projection shortcuts, a four-level encoder/decoder with skip connections, transposed-convolution upsampling, a per-stage dropout schedule that peaks at the bottleneck, and Kaiming initialization. The blocks, shortcut logic, channel bookkeeping, and weight initialization are all hand-written rather than assembled from a library.
Wrote a custom decoder around a pretrained remote-sensing encoder, extracting the five ResNet-50 feature scales and fusing them through bilinear-resize decoder blocks that tolerate arbitrary input sizes. The wrapper defends against upstream API drift (timm-style act1 versus torchvision-style relu) and inserts a 1×1 input adapter when the channel count does not match the encoder.
Designed a model factory so that five architecturally different networks share one training entry point, one evaluation entry point, and one inference entry point.

Code: model.py, model_torchgeo.py, model_factory.py

3. Training under severe class imbalance

Implemented multi-class Focal Loss and soft Dice Loss from first principles, both with correct ignore_index handling. This is a subtle correctness issue that silently degrades results when done wrong: ignored pixels are masked before reduction, one-hot targets are masked inside the Dice denominator, and a fully ignored batch short-circuits to a zero loss that still carries the right gradient shape.
Combined them into a tunable ComboLoss with self-normalizing weights, then swept the focal/dice balance across nine experiments to find the optimum at 0.2 / 0.8.
Derived inverse-square-root class-frequency alpha weights from measured dataset statistics, and empirically established where they help (rare-class behaviour) and where they do not (overall mIoU).

Code: losses.py

4. Training infrastructure

One configurable CLI trainer serving all five architectures: AdamW with named parameter groups and per-group learning rates, ReduceLROnPlateau driven by validation mIoU, early stopping with a minimum-delta threshold, best/last checkpointing, per-epoch CSV history logging, and a serialized run config for reproducibility.
Implemented discriminative fine-tuning for the self-supervised encoder (encoder 1e-4, decoder 1e-3) and optional staged encoder freezing, including the detail that a frozen encoder must also be held in eval() mode so that BatchNorm running statistics from pretraining are not quietly destroyed during decoder warm-up.
Deterministic seeding across random, NumPy, and PyTorch, with a fixed-seed split so that all 27 experiments are compared on exactly the same train/validation/test partition.
Geometry-preserving joint augmentation of image and mask (flips, 90° rotations, random crops) plus radiometric jitter applied to imagery only, including the contiguity fix that torch.from_numpy requires after np.flip and np.rot90 produce negative strides.

Code: train.py, dataset.py

5. Evaluation rigour

Built a confusion-matrix-based metric tracker that accumulates over the entire dataset with bincount rather than averaging per-batch metrics. Per-batch IoU averaging is a common and quietly optimistic mistake in segmentation reporting. Classes absent from both prediction and ground truth are excluded from the mean instead of contributing a spurious zero.
Selected models on mIoU rather than pixel accuracy, precisely because pixel accuracy is misleading on a dataset where two classes cover more than 60% of the labeled pixels.
Ran a structured 27-experiment program changing one variable at a time, including an explicit reproducibility re-run (exp15 repeated as exp17) after refactoring the model factory.

Code: metrics.py, evaluate.py

6. Model fusion

Implemented and benchmarked four fusion strategies — majority voting, mIoU-weighted hard voting, per-class-IoU-weighted voting, and temperature-scaled weighted soft-probability averaging — across growing subsets of the candidate pool, and established which one wins and why.
Solved the practical problem of ensembling models with different input specifications: two members consume 13-band tensors and three consume 10-band tensors, so the fusion dataset serves paired dual-band inputs for the same chip on the same spatial grid and dispatches each tensor to the correct member.
Made ties in majority voting deterministic and quality-aware through a mIoU-scaled tie-break term instead of depending on argmax index ordering.

Code: ensemble_predict_test.py

7. Uncertainty quantification

Went beyond hard labels to produce an explicit view of where the map should be trusted, using the deep-ensemble decomposition of predictive entropy into aleatoric-like and epistemic-like components.
Implemented weight-consistent uncertainty: the same mIoU weights used for probability averaging are also used for the expected-entropy term, so the decomposition stays internally coherent for a weighted ensemble. Most reference implementations assume uniform weights and skip this.
Added domain-aware class suppression: classes that are physically impossible for the target scene (snow and ice in a Lithuanian summer) are zeroed and the distribution is renormalized before entropy is computed, so uncertainty is measured over the genuinely active class space and normalized by the correct log C.
Exported every uncertainty layer as a georeferenced float32 GeoTIFF with an explicit nodata value and a descriptive band name, so the products are directly consumable in any GIS.

Code: predict_real_ensemble_chips.py

8. Production geospatial inference

Took the models from benchmark evaluation to a real, unlabeled area of interest: tiled inference across a whole municipality, chip-level GeoTIFF outputs, mosaicking back into single-raster products, --skip-existing resumability for long runs, and correct per-product nodata conventions (255 for class labels, -9999.0 for float uncertainty layers).
Applied the BOA offset correction required for post-2022 Sentinel-2 products so that the 2015, 2020, and 2025 composites are radiometrically comparable. Without it, the change analysis would partly measure a processing-baseline change rather than a land-cover change.
Handled the domain gap explicitly: training labels come from 2018 Europe-wide chips, while inference runs on a different city, different years, and a different acquisition path.

9. Communication and scientific honesty

Cartographic-quality reporting: hillshaded, boundary-clipped poster maps with class-share bars, scale bars, north arrows, and correct Copernicus attribution.
Results are reported with their limitations attached. The failure on the rare Natural Bare Ground class is stated plainly, model-derived percentages are explicitly not presented as official statistics, and the interpretation of the uncertainty decomposition cites the literature that questions it.
Every design decision in this document is traceable to a numbered experiment.

10. Software engineering practice

Type-annotated Python throughout, dataclass-based configuration, defensive shape and range validation at every module boundary, lazy imports so that optional heavy dependencies (torchgeo, segmentation_models_pytorch) are only required when actually used, an argparse CLI for every pipeline stage, CSV and JSON artifacts for every run, and a clean separation between src/ as an importable library and scripts/ as the executable pipeline.

Objective and tasks
System architecture
Dataset
Processed datasets
Project structure
Model architectures
Loss function
Metrics
Experiments
Main findings from individual models
Ensemble prediction
Prediction visualizations
Final ensemble configuration
Application to real Sentinel-2 data: Vilnius land-cover dynamics
- Vilnius visualizations
Ensemble uncertainty analysis for Vilnius, 2025
Running the pipeline
Limitations
References

Objective and tasks

The main research objective is to estimate how the land cover of Vilnius changes over time by applying deep learning segmentation models to Sentinel-2 summer median composites from different years. In particular, the project investigates whether urban expansion can be detected from satellite imagery and which land-cover classes are most affected by this change.

To support this objective, the work follows two connected directions. First, several segmentation models are trained and evaluated on the processed LandCoverNet Europe 2018 dataset. This includes a fully custom Residual U-Net, ImageNet-pretrained segmentation architectures, and a remote-sensing-specific Sentinel-2 pretrained model. Second, the best-performing models are combined into an ensemble and applied to real Sentinel-2 imagery over Vilnius for 2015, 2020, and 2025.

The main model selection metric is mean Intersection over Union (mIoU), while additional metrics such as pixel accuracy, Dice score, and per-class IoU are used for detailed analysis.

The project includes the following main tasks:

Prepare LandCoverNet Europe 2018 data as summer median Sentinel-2 composites.
Develop and train a fully custom Residual U-Net segmentation model.
Fine-tune general-purpose pretrained segmentation models on the processed dataset.
Fine-tune a remote-sensing-specific model using Sentinel-2 self-supervised pretrained weights.
Compare individual models using validation and test segmentation metrics.
Build and evaluate an ensemble prediction strategy based on soft probability averaging.
Apply the selected ensemble to Sentinel-2 imagery over Vilnius for multiple years.
Estimate land-cover class shares and analyze urban development trends over time.
Analyze ensemble prediction uncertainty using predictive entropy, expected entropy, and mutual information.

System architecture

flowchart TD
    A["Raw Sentinel-2 archive<br/>LandCoverNet Europe 2018"] --> B["SCL cloud / shadow / snow masking"]
    B --> C["Per-pixel summer median composite<br/>June, July, August"]
    C --> D["Processed GeoTIFF chips<br/>4-band / 10-band / 13-band"]
    D --> E["70 / 15 / 15 split, seed 42"]
    E --> F["Training<br/>ComboLoss = Focal + Dice"]
    F --> G1["Residual U-Net<br/>custom"]
    F --> G2["DeepLabV3+ / U-Net++ / FPN<br/>ImageNet encoders"]
    F --> G3["ResNet-50 U-Net<br/>Sentinel-2 DINO encoder"]
    G1 --> H["Test evaluation<br/>mIoU / Dice / per-class IoU"]
    G2 --> H
    G3 --> H
    H --> I["Weighted soft-probability ensemble<br/>5 members, mIoU weights"]
    I --> J["Real Sentinel-2 AOI inference<br/>Vilnius 2015 / 2020 / 2025"]
    J --> K1["Land-cover mosaic<br/>GeoTIFF"]
    J --> K2["Uncertainty mosaics<br/>entropy / MI / confidence"]
    K1 --> L["Class-share statistics<br/>and change analysis"]
    K1 --> M["Cartographic poster maps"]
    K2 --> M

Dataset

The source dataset is LandCoverNet Europe 2018, part of the LandCoverNet global land-cover classification training dataset [1].

Dataset DOI: https://doi.org/10.34911/rdnt.63fxe5
Source page: https://source.coop/radiantearth/landcovernet/landcovernet_eu

LandCoverNet provides land-cover labels for multi-spectral satellite imagery from Sentinel-1, Sentinel-2, and Landsat-8 for 2018. In this project, the experiments use Sentinel-2 imagery and LandCoverNet labels for Europe.

Classes

The original LandCoverNet labels were remapped for training as follows:

Original label	Train ID	Class
0	255	Ignore / no data
1	0	Water
2	1	Artificial Bare Ground
3	2	Natural Bare Ground
4	3	Permanent Snow and Ice
5	4	Woody Vegetation
6	5	Cultivated Vegetation
7	6	Natural Grassland

Pixel class distribution

The following distribution was used for class-frequency weighting experiments. Values are calculated over valid labeled pixels.

Train ID	Class	Pixel fraction	Pixel percent
0	Water	0.046550	4.6%
1	Artificial Bare Ground	0.054580	5.4%
2	Natural Bare Ground	0.010020	1.0%
3	Permanent Snow and Ice	0.014138	1.4%
4	Woody Vegetation	0.265378	26.5%
5	Cultivated Vegetation	0.371394	37.1%
6	Natural Grassland	0.237867	23.8%

The dataset is strongly imbalanced. The rarest valid classes are Natural Bare Ground and Permanent Snow and Ice. The two dominant classes together account for roughly 64% of all labeled pixels, which is exactly why mIoU rather than pixel accuracy drives model selection throughout this project.

Processed datasets

The raw LandCoverNet data was converted into summer median Sentinel-2 composites. The preprocessing pipeline:

uses scenes from June, July, and August;
applies cloud / shadow / snow masking using the Sentinel-2 Scene Classification Layer (SCL);
computes a per-pixel median composite;
saves processed GeoTIFF images and masks.

Masked pixels are set to NaN before compositing, so nanmedian ignores them entirely instead of letting a cloud pull the composite toward bright values. Using the median rather than the mean also makes the composite robust to residual contamination that the SCL mask misses.

Bad SCL classes removed during preprocessing:

SCL ID	Meaning
0	No data
1	Saturated / defective
3	Cloud shadows
8	Cloud medium probability
9	Cloud high probability
10	Thin cirrus
11	Snow / ice

Three processed dataset variants were used:

Dataset	Bands	Channels	Image shape	Purpose
4-band	B02, B03, B04, B08	4	4 × 256 × 256	Early baseline
10-band	B02, B03, B04, B08, B05, B06, B07, B8A, B11, B12	10	10 × 256 × 256	Main dataset for custom and SMP models
13-band Sentinel-2 order	B01, B02, B03, B04, B05, B06, B07, B08, B8A, B09, B10, B11, B12	13	13 × 256 × 256	Sentinel-2 self-supervised pretrained models

For the 13-band dataset, B10 is a zero-filled synthetic channel because the raw LandCoverNet Sentinel-2 data does not include this band. Keeping the channel in its correct spectral position allows 13-channel pretrained weights to load unchanged.

Project structure

src/
├── dataset.py          # GeoTIFF dataset, normalization modes, joint augmentation, splits
├── evaluate.py         # Test-set evaluation CLI, confusion matrix export
├── losses.py           # Focal / Dice / Combo losses with ignore_index handling
├── metrics.py          # Confusion-matrix metric tracker (accuracy, IoU, Dice)
├── model.py            # Custom Residual U-Net
├── model_factory.py    # Unified constructor for all five architectures
├── model_torchgeo.py   # Sentinel-2 pretrained ResNet-50 encoder + U-Net decoder
└── train.py            # Training CLI, parameter-group LRs, scheduling, early stopping

scripts/
├── prepare_landcovernet_median_dataset_multiband.py  # Raw archive -> median composites
├── inspect_processed_dataset.py                      # Visual sanity check of one chip
├── check_processed_dataset.py                        # Dataset integrity report -> CSV
├── visualize_model_comparison.py                     # Side-by-side model comparison figures
├── ensemble_predict_test.py                          # Benchmark all fusion strategies
└── predict_real_ensemble_chips.py                    # AOI inference + uncertainty + mosaics

data/
├── outputs/
│   └── processed_dataset_check.csv
└── figures/
    ├── comparison_8_samples_with_best_ensemble_seed42.png
    ├── comparison_8_samples_with_best_ensemble_seed1707.png
    ├── comparison_8_samples_with_best_ensemble_seed1736.png
    ├── comparison_8_samples_with_best_ensemble_seed1777.png
    ├── confusion_matrix_best_ensemble_soft_avg_percent.png
    ├── vilnius_2015_prediction_poster.png
    ├── vilnius_2020_prediction_poster.png
    ├── vilnius_2025_prediction_poster.png
    ├── vilnius_2025_rgb_poster.png
    └── vilnius_2025_prediction_and_uncertainty.png

Model architectures

Custom Residual U-Net

The main custom architecture is a Residual U-Net inspired by the original U-Net encoder-decoder structure with skip connections [2] and residual learning [3].

The implemented model uses:

encoder blocks with ResidualDoubleConv, where the shortcut becomes a 1×1 convolution with BatchNorm whenever the channel count changes;
a bottleneck residual convolutional block at double the deepest encoder width;
decoder blocks with transposed convolutions and skip-connection concatenation;
a final 1×1 convolution producing 7-class logits at full input resolution;
Kaiming normal initialization for all convolutions and unit/zero initialization for BatchNorm.

The strongest custom configuration used:

base_features = 64
features = (64, 128, 256, 512)
dropout_encoder = (0.0, 0.05, 0.1, 0.2)
dropout_bottleneck = 0.3
dropout_decoder = (0.2, 0.1, 0.05, 0.0)

The dropout schedule is deliberately asymmetric. Regularization ramps up toward the bottleneck, where the representation is most abstract and most prone to overfitting, and ramps back down toward the output, where the network needs full capacity to recover sharp class boundaries.

ImageNet-pretrained segmentation models

Several pretrained segmentation models were tested through the segmentation_models_pytorch library, which supplies the architectures and the ImageNet encoder weights.

Models used:

Model	Encoder	Pretraining	Architecture / encoder references
DeepLabV3+	ResNet34 / ResNet50	ImageNet	[4], [3]
U-Net++	EfficientNet-B3	ImageNet	[5], [6]
FPN	EfficientNet-B3	ImageNet	[7], [6]

These models were selected for architectural diversity rather than for individual accuracy: an atrous-convolution model, a densely nested skip-connection model, and a feature-pyramid model make different errors, which is exactly what a probability-averaging ensemble needs.

Sentinel-2 self-supervised pretrained model

A Sentinel-2 pretrained ResNet-50 encoder from the TorchGeo model zoo was also tested [8], using a ResNet-50 backbone [3].

The best setup used:

ResNet50 encoder, Sentinel-2 pretrained
weights = SENTINEL2_ALL_DINO
input channels = 13
normalization = reflectance
decoder/head LR = 1e-3
encoder LR = 1e-4

The encoder is combined with a lightweight custom U-Net-like decoder. Encoder features are taken at five scales (stem, layer1–layer4), and each decoder block resizes bilinearly to the skip-connection resolution before concatenating, so the model accepts non-power-of-two input sizes such as the 192-pixel training crops.

The key result here is that a domain-specific self-supervised encoder only pays off when the fine-tuning schedule respects it: the pretrained encoder needs an order-of-magnitude lower learning rate than the randomly initialized decoder, otherwise the first epochs destroy the pretrained representation.

Loss function

The main loss function is a weighted combination of Focal Loss and Dice Loss:

ComboLoss = focal_weight × FocalLoss + dice_weight × DiceLoss

The best general setting was:

focal_weight = 0.2
dice_weight  = 0.8
gamma = 2.0
ignore_index = 255

Some experiments additionally used class-frequency alpha weights for Focal Loss, computed as the normalized inverse square root of the measured class frequencies.

The two terms play complementary roles. Focal Loss operates per pixel and down-weights easy, already-correct pixels so that the gradient signal keeps coming from hard boundaries. Dice Loss operates per class over the whole batch and is inherently scale-invariant, so a rare class contributes as much to the loss as a dominant one. The Dice-heavy 0.8 weighting is what the ablation converged on, and it is consistent with mIoU being the selection metric.

Metrics

The following metrics were used:

Metric	Meaning
Pixel Accuracy	Fraction of correctly classified valid pixels
IoU per class	Intersection over Union for each class
Mean IoU	Average IoU over valid classes
Dice per class	Dice coefficient for each class
Mean Dice	Average Dice over valid classes

All metrics are computed from a confusion matrix accumulated over the full dataset, never averaged batch by batch. The main model-selection metric was Mean IoU (mIoU) because it is far more robust than pixel accuracy under class imbalance.

Experiments

All experiments were evaluated on the same 70/15/15 train/validation/test split with seed=42.

Full experiment table

Exp	Dataset	Model	Main setup	Test Acc	Test mIoU	Test mDice	Notes
1	4-band	ResUNet	no crop, no alpha, Focal/Dice 0.5/0.5, 30 epochs	0.7028	0.4650	0.5887	baseline
2	4-band	ResUNet	crop224, no alpha, 0.5/0.5	0.7247	0.4783	0.5955	crop improved result
3	4-band	ResUNet	crop224, alpha, 0.5/0.5	0.6949	0.4871	0.6144	alpha improved mIoU/mDice
4	4-band	ResUNet	crop192, alpha, 0.5/0.5	0.7002	0.4893	0.6152	crop192 better than crop224
5	4-band	ResUNet	crop128, alpha, 0.5/0.5	0.6859	0.4891	0.6154	close to exp4
6	4-band	ResUNet	crop192, alpha, 0.4/0.6	0.6955	0.4833	0.6105	worse
7	4-band	ResUNet	crop192, alpha, 0.3/0.7	0.7123	0.4900	0.6145	better than exp6
8	4-band	ResUNet	crop192, alpha, 0.2/0.8	0.7180	0.4947	0.6167	best 4-band
9	4-band	ResUNet	crop192, alpha, 0.1/0.9	0.7077	0.4886	0.6156	too much Dice weight
10	10-band	ResUNet f32	alpha, crop192, 0.2/0.8	0.7618	0.5445	0.6723	strong gain from 10 bands
11	10-band	ResUNet f32	alpha, crop192, 0.5/0.5	0.7439	0.5280	0.6591	worse than Dice-heavy
12	10-band	ResUNet f32	no alpha, crop192, 0.2/0.8	0.7524	0.5482	0.6750	no alpha improved mIoU
13	10-band	ResUNet f48	no alpha, crop192, 0.2/0.8	0.7523	0.5384	0.6655	wider model did not help
14	10-band	ResUNet f64	no alpha, crop192, 0.2/0.8	0.7612	0.5479	0.6778	wider, but no clear gain
15	10-band	ResUNet f64	no alpha, crop192, stronger dropout, 0.2/0.8	0.7573	0.5536	0.6827	best early custom model
16	10-band	ResUNet f64	no alpha, crop224, stronger dropout, 0.2/0.8	0.7583	0.5498	0.6826	crop224 slightly worse
17	10-band	ResUNet f64	repeat of exp15 after model factory changes	0.7599	0.5486	0.6766	reproducibility check
18	10-band	DeepLabV3+	ResNet34, ImageNet, crop192, no alpha, 0.2/0.8	0.7412	0.5117	0.6324	weaker than ResUNet
19	10-band	U-Net++	EfficientNet-B3, ImageNet, crop192, no alpha, 0.2/0.8	0.7605	0.5365	0.6571	best ImageNet-pretrained baseline
20	10-band	FPN	EfficientNet-B3, ImageNet, crop192, no alpha, 0.2/0.8	0.7588	0.5104	0.6336	useful diversity for ensemble
21	10-band	DeepLabV3+	ResNet50, ImageNet, crop192, no alpha, 0.2/0.8	0.7413	0.5237	0.6470	better than ResNet34 variant
22	13-band	ResUNet f64	13-band input, reflectance norm	0.7516	0.5429	0.6706	13 bands alone did not improve custom model
23	13-band	ResNet50 U-Net (S2 SSL)	Sentinel-2 ALL DINO, same LR	0.7457	0.5198	0.6395	weak baseline
24	13-band	ResNet50 U-Net (S2 SSL)	Sentinel-2 ALL DINO, encoder LR 1e-4, decoder LR 1e-3	0.7846	0.5602	0.6815	best single model overall
25	13-band	ResNet50 U-Net (S2 SSL)	separate LR + Sentinel-2 stats norm	0.7697	0.5404	0.6561	S2 stats norm hurt
26	13-band	ResNet50 U-Net (S2 SSL)	separate LR + S2 stats norm + freeze encoder 5 epochs	0.7751	0.5391	0.6495	freezing did not help
27	13-band	ResNet50 U-Net (S2 SSL)	separate LR + reflectance norm + class alpha	0.7742	0.5583	0.6815	improved rare class vs exp24

Main findings from individual models

Adding more Sentinel-2 bands was the largest single improvement. The best 4-band model reached mIoU = 0.4947, while the best 10-band custom model reached mIoU = 0.5536. Red-edge and SWIR bands carry vegetation and moisture information that simply is not recoverable from RGB and NIR.
The custom Residual U-Net remained competitive. ImageNet-pretrained segmentation models did not outperform the best custom Residual U-Net. Natural-image priors transfer poorly to 10-band top-of-canopy reflectance.
Sentinel-2 self-supervised pretraining helped, but only with careful fine-tuning. Using a lower learning rate for the pretrained encoder and a higher learning rate for the decoder/head was decisive: the same weights scored 0.5198 with a shared LR (exp23) and 0.5602 with split LRs (exp24).
Sentinel-2 mean/std normalization did not help in this setup. Simple reflectance normalization performed better for the processed summer median composites, most likely because the published statistics were computed on single-date scenes rather than seasonal medians.
Class alpha improved some rare-class behavior but did not improve the best overall mIoU. It is a recall/precision trade, not a free win.
Model capacity was not the bottleneck. Widening the custom model from 32 to 48 to 64 base features produced no consistent gain, whereas input bands, loss balance, and fine-tuning schedule all did.

Ensemble prediction

After evaluating individual models, several ensemble strategies were tested.

Candidate models:

Exp	Model	Input
22	ResUNet 13-band	13-band
27	ResNet50 U-Net, S2 DINO + alpha	13-band
19	U-Net++ EfficientNet-B3	10-band
21	DeepLabV3+ ResNet50	10-band
20	FPN EfficientNet-B3	10-band

For the same chip, both 10-band and 13-band inputs are loaded from paired dataset roots. Each model predicts logits or masks on the same spatial grid, so the fusion happens pixel to pixel with no resampling.

Tested ensemble methods

Method	Description
Majority voting	Each model votes for one class per pixel; ties are broken by model mIoU
Weighted hard voting	Each model vote is weighted by its test mIoU
Class-aware hard voting	Each model vote for class `c` is weighted by the model's class-specific IoU for class `c`
Soft probability averaging	Model logits are converted to probabilities with softmax; probabilities are averaged with model-level mIoU weights

The best method was weighted soft probability averaging.

Ensemble results

Models	Strategy	Pixel Acc	mIoU	mDice
exp22 + exp27	majority	0.7652	0.5538	0.6826
exp22 + exp27	weighted	0.7742	0.5583	0.6815
exp22 + exp27	class-aware	0.7652	0.5389	0.6628
exp22 + exp27	soft_avg	0.7777	0.5606	0.6851
exp22 + exp27 + exp19	majority	0.7779	0.5651	0.6906
exp22 + exp27 + exp19	weighted	0.7784	0.5640	0.6879
exp22 + exp27 + exp19	class-aware	0.7776	0.5472	0.6621
exp22 + exp27 + exp19	soft_avg	0.7850	0.5707	0.6950
exp22 + exp27 + exp19 + exp21	majority	0.7791	0.5669	0.6924
exp22 + exp27 + exp19 + exp21	weighted	0.7790	0.5676	0.6928
exp22 + exp27 + exp19 + exp21	class-aware	0.7770	0.5471	0.6620
exp22 + exp27 + exp19 + exp21	soft_avg	0.7867	0.5711	0.6945
exp22 + exp27 + exp19 + exp21 + exp20	majority	0.7877	0.5746	0.6994
exp22 + exp27 + exp19 + exp21 + exp20	weighted	0.7880	0.5752	0.6999
exp22 + exp27 + exp19 + exp21 + exp20	class-aware	0.7821	0.5399	0.6568
exp22 + exp27 + exp19 + exp21 + exp20	soft_avg	0.7938	0.5816	0.7064

Two patterns are worth reading out of this table. First, soft averaging wins at every ensemble size, and its lead does not shrink as members are added. Second, adding FPN (exp20) improves the ensemble even though exp20 is the weakest individual model in the pool at mIoU = 0.5104. That is the ensemble working as intended: what matters is decorrelated errors, not individual rank.

Best ensemble

The best final ensemble is:

Models:
exp22 + exp27 + exp19 + exp21 + exp20

Strategy:
weighted soft probability averaging

Result:
Pixel Acc = 0.7938
mIoU      = 0.5816
mDice     = 0.7064

Compared with the best single model from the final set, exp27:

Model	Pixel Acc	mIoU	mDice
exp27	0.7742	0.5583	0.6815
Best ensemble soft_avg	0.7938	0.5816	0.7064
Improvement	+0.0196	+0.0233	+0.0249

Best ensemble per-class metrics

The final selected ensemble was evaluated per class as follows:

Train ID	Class	IoU	Dice
0	Water	0.8517	0.9199
1	Artificial Bare Ground	0.6946	0.8198
2	Natural Bare Ground	0.1589	0.2742
3	Permanent Snow and Ice	0.3947	0.5660
4	Woody Vegetation	0.6864	0.8140
5	Cultivated Vegetation	0.7662	0.8676
6	Natural Grassland	0.5187	0.6830

The remaining weakness is the rare Natural Bare Ground class. However, the ensemble improves several important classes compared with the best single model, including Cultivated Vegetation, Natural Grassland, and Permanent Snow and Ice.

Confusion matrix

For interpretability, the final ensemble confusion matrix is visualized as a row-normalized percentage matrix. Each row corresponds to a ground-truth class and sums to 100%. The diagonal therefore shows the percentage of pixels from each true class that were correctly classified, while off-diagonal cells show exactly where each class leaks.

How to read this figure. Rows are the ground-truth class, columns are the predicted class, cell values are the percentage of that row's pixels, and the viridis colour scale runs from 0% (dark purple) to 100% (yellow). Because rows are normalized independently, a rare class is shown on the same footing as a dominant one, which is what makes the failure modes visible at all.

What the matrix shows, class by class:

Water — 95.4% correct. The brightest cell on the diagonal. Water has a distinctive low-reflectance spectral signature across the NIR and SWIR bands, and it is spatially compact, so it is essentially solved. Its 2.0% leak into Snow/Ice is the classic bright-water and thin-cloud-over-water confusion.
Woody Vegetation — 86.8% and Cultivated Vegetation — 86.4%. The two large vegetation classes are both handled well. Their residual error goes almost entirely into each other and into grassland (3.1% and 9.5% for woody; 3.1% and 7.9% for cultivated), which is the expected behaviour at 10 m resolution where a field edge or a tree line falls inside a single pixel.
Artificial Bare Ground — 84.3%. Strong, which matters because this is the class that carries the urban-expansion analysis later in this document. Its errors go mainly to Cultivated Vegetation (7.8%) and Natural Grassland (5.9%), that is, to the green space interleaved with buildings in low-density residential areas.
Natural Grassland — 63.7%. The weakest of the common classes. Nearly a fifth of its pixels (19.1%) are predicted as Woody Vegetation and a further 11.3% as Cultivated Vegetation. In a summer median composite, unmanaged grassland, meadow, and a mown field can all sit in a very similar part of the spectral space, and the label boundary between them is partly a management distinction rather than a physical one.
Permanent Snow and Ice — 78.1%. Surprisingly high for a class covering only 1.4% of pixels, because where snow survives a European summer it is spectrally unambiguous.
Natural Bare Ground — 20.4%. The clear failure, and the figure explains why rather than just reporting it. Its pixels are scattered almost evenly across three other classes: 26.7% to Snow/Ice, 26.4% to Cultivated Vegetation, and 15.0% to Water. That signature is characteristic of a class with too few training examples (1.0% of pixels) and no coherent spectral identity: bright dry sand resembles snow, sparsely vegetated bare soil resembles cropland, and wet or shadowed bare ground resembles water. This is a data problem, not a capacity problem, and no amount of architecture tuning fixed it.

The class-aware hard voting strategy did not work well, and this matrix is the reason. Weighting each vote by the per-class IoU concentrates the decision on whichever model happens to score best for a class, which is far too sensitive when that per-class score is itself unreliable — as it is at IoU = 0.16 for Natural Bare Ground. Soft probability averaging preserves the full confidence distribution instead, and produces the best overall segmentation quality.

Prediction visualizations

The following figures place the individual models and the final ensemble side by side on identical test chips, so that differences in behaviour are visible directly rather than only through aggregate metrics.

Each figure has the same layout: 8 rows are 8 randomly drawn test chips (labelled on the left with their LandCoverNet chip ID, for example 32TPT_03 or 35UMU_07), and 8 columns are, from left to right:

RGB — a true-colour rendering of the summer median composite, built from B04/B03/B02 with a 2–98 percentile stretch per band;
Ground truth — the LandCoverNet reference mask;
ResUNet 13b (exp22) — the custom Residual U-Net;
DeepLabV3+ R50 (exp21);
Unet++ EffB3 (exp19);
FPN EffB3 (exp20);
TorchGeo DINO α (exp27) — the Sentinel-2 self-supervised model;
Best ensemble (soft avg) — the final weighted soft-probability fusion.

The class colour key runs along the bottom of each figure: blue = Water, brown = Artificial Bare Ground, tan = Natural Bare Ground, white = Permanent Snow and Ice, dark green = Woody Vegetation, yellow-green = Cultivated Vegetation, mid green = Natural Grassland.

Three seeds are shown so that the comparison is not cherry-picked from a single lucky draw of chips. Each seed selects a different random sample of 8 chips from the same held-out test split. Figures are stored in data/figures/.

Seed 1707

Reading across the rows in this figure makes several individual-model behaviours obvious. In the coastal/estuary rows, FPN produces visibly noisier, speckled output and over-predicts Snow/Ice on bright water, while DeepLabV3+ tends to smooth small structures away entirely. The Sentinel-2 self-supervised model resolves narrow linear features such as river channels and field tracks more crisply than the ImageNet-pretrained models, which is the expected payoff of pretraining on the same sensor. In the agricultural rows, the models disagree most about the grassland/cropland split — exactly the confusion the matrix above quantifies at 11.3% — with some models painting a whole field as cultivated and others breaking it into grassland patches. The rightmost ensemble column consistently sits closest to the ground-truth column: it keeps the sharp features that the self-supervised model finds, without inheriting the speckle that FPN introduces.

Seed 1736

Seed 1777

Across all three seeds the qualitative pattern is stable. The ensemble produces smoother and more spatially coherent predictions than any individual model, while still preserving useful minority-class regions that a simple majority vote would have erased. Where all five members agree, the ensemble output is effectively identical to theirs; the visible benefit appears precisely in the mixed and transitional areas where they disagree — and those are the same regions that light up in the uncertainty maps later in this document.

Final ensemble configuration

The best result is obtained not by a single model, but by a heterogeneous ensemble that combines:

a custom Residual U-Net trained on 13-band data;
a Sentinel-2 DINO self-supervised pretrained model;
three ImageNet-pretrained segmentation architectures trained on 10-band Sentinel-2 composites.

The five members differ along three axes at once — input band count, encoder pretraining domain, and decoder topology — which is what makes their errors decorrelated enough for fusion to pay off.

The final weighted soft probability ensemble achieved:

Pixel Accuracy = 0.7938
Mean IoU       = 0.5816
Mean Dice      = 0.7064

Application to real Sentinel-2 data: Vilnius land-cover dynamics

After model selection, the best ensemble was applied to real Sentinel-2 L2A imagery for the Vilnius municipality. For each target year, a summer median Sentinel-2 composite was generated, tiled into 256 × 256 chips, processed by the ensemble, and mosaicked back into a georeferenced land-cover map.

The real-data inference pipeline uses:

Sentinel-2 L2A imagery accessed through Microsoft Planetary Computer;
13-band Sentinel-2 composites at 10 m resolution;
BOA offset correction for post-2022 Sentinel-2 products, so that all three years remain radiometrically comparable;
the best five-model ensemble with soft probability averaging;
suppression of the Permanent Snow and Ice class, which is physically impossible in a Lithuanian summer composite, with probabilities renormalized over the remaining six classes;
city-boundary clipping for final visualization and area statistics.

The analyzed years were 2015, 2020, and 2025. The table below reports the predicted share of the main land-cover classes inside the Vilnius boundary. Values are percentages of classified pixels within the city boundary.

Year	Natural Grassland	Cultivated Vegetation	Woody Vegetation / Forest	Artificial Bare Ground	Water
2015	6.4%	18.8%	43.8%	29.6%	1.4%
2020	7.1%	15.5%	43.8%	32.2%	1.3%
2025	8.1%	12.2%	43.8%	34.5%	1.4%

Change summary

Change	2015 -> 2020	2020 -> 2025	2015 -> 2025
Natural Grassland	+0.7 pp	+1.0 pp	+1.7 pp
Cultivated Vegetation	-3.3 pp	-3.3 pp	-6.6 pp
Woody Vegetation / Forest	0.0 pp	0.0 pp	0.0 pp
Artificial Bare Ground	+2.6 pp	+2.3 pp	+4.9 pp
Water	-0.1 pp	+0.1 pp	0.0 pp

The results suggest a steady expansion of urbanized surfaces in Vilnius. The predicted artificial class increased from 29.6% in 2015 to 34.5% in 2025, which corresponds to approximately +4.9 percentage points over ten years, or roughly 0.5 percentage points per year. The main decrease is observed in cultivated vegetation, which dropped from 18.8% to 12.2%. In contrast, predicted forest cover remained stable at 43.8% across the analyzed years.

This indicates that urban expansion in the analyzed period was predicted mostly at the expense of cultivated or open agricultural areas rather than forested areas. In qualitative map inspection, major green and park-like zones remain largely preserved, while artificial surfaces expand around already urbanized or peri-urban parts of the municipality.

These values should be interpreted as model-derived estimates, not as official land-cover statistics. The analysis is still useful for showing the practical transfer of the trained ensemble from benchmark evaluation to multi-year remote-sensing inference over a real urban area of interest.

Vilnius visualizations

The final maps are stored in data/figures/. All four poster figures share one cartographic template, so they can be compared directly:

a dark background with the raster clipped to the Vilnius municipal boundary, drawn as a dark red outline;
DEM-based hillshade blended under the classification, which is why terrain relief, valley slopes, and the river corridor remain visible through the flat class colours;
a stacked 100% class-share bar on the left, giving the exact percentage of every class in that year without needing the table above;
a method annotation block stating the input (Sentinel-2 summer median composite, 10 m, 13-band) and the model (ensemble of five segmentation models, fused by soft probability averaging);
a 5 km scale bar and Copernicus attribution.

Note that the poster palette is deliberately different from the benchmark-figure palette used in the model comparison section above. On the posters: red = Artificial, dark green = Forest, yellow = Cultivated, light green = Grassland, blue = Water. Snow/Ice and Natural Bare Ground do not appear because Snow/Ice is suppressed for summer inference and Natural Bare Ground is effectively absent from the area of interest.

Sentinel-2 RGB reference, 2025

This is the model input, not a model output, and it is included so the predictions can be checked against reality by eye. It is the true-colour rendering of the same cloud-free summer median composite the ensemble consumed for 2025, at 10 m per pixel, clipped to the same municipal boundary.

Even without any classification, the structure the model has to recover is visible here: the dense grey-white core of the old town and the central business district; the regular geometry of the Soviet-era apartment districts spreading north-west and south; the dark, continuous forest blocks in the east and north-east; the pale rectangular field patterns in the south-west; and the Neris river threading through the middle of the city. Comparing this panel against the 2025 prediction below is the fastest sanity check available — every red area on the prediction map should correspond to a grey or bright built-up surface here.

Predicted land cover, 2015

The baseline year. The share bar reads Artificial 29.6%, Forest 43.8%, Cultivated 18.8%, Grassland 6.4%, Water 1.4%. Note how much yellow (cultivated land) is present in this map, particularly in the southern and south-western sectors and around the outer edges of the built-up area. The red urban mass is concentrated in the centre and along the main radial corridors, and there are visible gaps of yellow and green between the outer districts. This map is the reference against which the following two years should be read.

Predicted land cover, 2020

The mid-point. Artificial has risen to 32.2% and cultivated has fallen to 15.5%. The change is not a uniform expansion of the city edge; it appears as an infill pattern, where the yellow gaps between existing red districts close first. This is consistent with development following existing infrastructure rather than leapfrogging into open countryside.

Predicted land cover, 2025

The end state. Artificial reaches 34.5% while cultivated drops to 12.2%. Placing this map next to the 2015 map is the clearest single statement of the project's finding: the red area has grown noticeably and the yellow area has shrunk by roughly a third of its 2015 extent, while the dark green forest blocks in the east and north-east are visually almost unchanged — matching the 43.8% forest share that holds constant across all three years.

One caveat worth stating alongside these maps. Part of the cultivated-to-artificial transition is genuine construction, but part of it is a classification boundary effect: as a peri-urban field becomes fragmented by roads, plots, and scattered houses, its pixels become mixed, and mixed pixels near built-up areas resolve toward the artificial class. The uncertainty analysis in the next section is what allows those two situations to be told apart, because genuine dense construction is predicted with low entropy while the fragmented transition zones are exactly where entropy is highest.

Ensemble uncertainty analysis for Vilnius, 2025

In addition to the final class prediction, the ensemble can also be used to estimate spatial prediction uncertainty. This is possible because the final method is based on weighted soft probability averaging, so the model does not only produce hard class labels, but also class probability distributions. This follows the general deep ensemble idea of combining probabilistic predictions from multiple models for predictive uncertainty estimation [9].

For each pixel ($x$), every ensemble member ($m$) produces a probability distribution over land-cover classes:

$$p_m(y=c \mid x)$$

where ($c$) is a land-cover class and ($m = 1, \dots, M$) is an ensemble model. In this project, the final ensemble contains five models:

$$M = 5$$

The ensemble-averaged probability for class ($c$) is calculated as a weighted average:

$$\bar{p}(y=c \mid x) = \frac{\sum_{m=1}^{M} w_m p_m(y=c \mid x)} {\sum_{m=1}^{M} w_m}$$

where ($w_m$) is the model-level weight based on test-set mIoU.

The final predicted class is then:

$$\hat{y}(x) = \arg\max_c \bar{p}(y=c \mid x)$$

To analyze uncertainty, three information-theoretic metrics were computed. Similar entropy-based uncertainty measures are commonly used in Bayesian deep learning and ensemble uncertainty estimation [10], [11], [12]. Here, the five ensemble models are used as a practical way to approximate uncertainty over model predictions.

Predictive entropy: total uncertainty

Predictive entropy is calculated from the final ensemble-averaged probability distribution:

$$H[\bar{p}(y \mid x)] = -\sum_{c=1}^{C} \bar{p}(y=c \mid x) \log \bar{p}(y=c \mid x)$$

This is the standard entropy of a categorical predictive distribution and is commonly used as a direct measure of predictive uncertainty in classification [10], [11].

The value is normalized by the maximum possible entropy:

$$H_{\text{pred,norm}} = \frac{H[\bar{p}(y \mid x)]}{\log C}$$

This gives values approximately in the range from 0 to 1. Note that ($C$) here is the number of active classes after domain-based class suppression, not the nominal seven, so the normalization stays correct when a class is excluded.

Predictive entropy can be interpreted as total predictive uncertainty. It is high when the final ensemble probability is distributed across several competing classes, and low when one class clearly dominates.

In the Vilnius map, low predictive entropy is visible over large homogeneous areas such as dense urban districts, water bodies, and major forested zones. These objects have more stable spectral and spatial patterns, so the ensemble predicts them more confidently.

Higher predictive entropy appears mostly in transitional and fragmented areas: suburban zones, edges of built-up areas, cultivated vegetation, natural grassland, river banks, and areas along major roads. This is expected because these locations often contain mixed pixels or gradual transitions between land-cover types.

Expected entropy: data / aleatoric-like uncertainty

Expected entropy is calculated by computing entropy for each model separately and then averaging these entropy values using the ensemble weights:

$$\mathbb{E}[H[p_m(y \mid x)]] = \frac{\sum_{m=1}^{M} w_m H[p_m(y \mid x)]} {\sum_{m=1}^{M} w_m}$$

where

$$H[p_m(y \mid x)] = -\sum_{c=1}^{C} p_m(y=c \mid x) \log p_m(y=c \mid x)$$

The normalized version is:

$$H_{\text{exp,norm}} = \frac{\mathbb{E}[H[p_m(y \mid x)]]}{\log C}$$

Using the same weights ($w_m$) here as in the probability average is what keeps the decomposition self-consistent; mixing weighted averaging with unweighted expected entropy would make the mutual-information term below meaningless.

Expected entropy shows how uncertain the individual models are on average. In this project, it is used as a practical proxy for aleatoric-like uncertainty, meaning uncertainty that comes from the data itself.

This interpretation follows the common uncertainty-decomposition view where expected conditional entropy captures the part of uncertainty associated with data ambiguity or aleatoric uncertainty [11], [12], [13].

This kind of uncertainty is especially important in Sentinel-2 land-cover segmentation. Many areas are not clean semantic objects with sharp boundaries. A single 10 m pixel can include trees, grass, roofs, roads, gardens, shadows, or river-bank vegetation. This is especially visible in private housing areas with dense vegetation, where small objects are mixed inside the same spatial unit.

In the Vilnius uncertainty maps, expected entropy is high in many of the same places where predictive entropy is high. This suggests that a large part of the uncertainty is caused by the physical and semantic ambiguity of the scene itself: fragmented land cover, mixed pixels, and transitions between visually or spectrally similar classes.

Mutual information: model / epistemic-like uncertainty

Mutual information is computed as the difference between predictive entropy and expected entropy:

$$MI(y, m \mid x) = H[\bar{p}(y \mid x)] - \frac{\sum_{m=1}^{M} w_m H[p_m(y \mid x)]}{\sum_{m=1}^{M} w_m}$$

Using normalized entropy values, this becomes:

$$MI_{\text{norm}} = H_{\text{pred,norm}} - H_{\text{exp,norm}}$$

This follows the same information-theoretic structure used in Bayesian active learning and Bayesian deep learning, where mutual information is expressed as the difference between predictive entropy and expected entropy and is used to measure model disagreement [10], [11], [12], [14].

Mutual information is used as a proxy for epistemic-like uncertainty, or model disagreement. It becomes high when individual models are confident but disagree with each other. This interpretation is common in uncertainty quantification, although recent work also discusses limitations of treating conditional entropy and mutual information as a perfect aleatoric/epistemic decomposition [15].

For example, if one model confidently predicts woody vegetation, another confidently predicts grassland, and a third confidently predicts cultivated vegetation, then the average ensemble distribution becomes uncertain even though each individual model is confident. This indicates model disagreement rather than only ambiguous input data.

In the Vilnius 2025 map, mutual information is relatively low over most of the area. This is an interesting result because it suggests that the ensemble members are generally consistent with each other. The models do not strongly disagree over most of the city. Instead, most uncertainty seems to come from mixed or transitional land-cover areas rather than from severe model disagreement.

The highest mutual information values appear mainly in complex zones where different architectures may interpret the same pixel slightly differently: suburban areas with small buildings and vegetation, edges between forest and grassland, cultivated-to-grassland transitions, and some fragmented peri-urban regions.

Visual uncertainty map

The figure below shows the final 2025 prediction together with the three uncertainty maps.

How to read this figure. It is a 2 × 2 panel grid over the full inference tile, so it extends beyond the municipal boundary; the bright green outline is the Vilnius boundary drawn over the raster, and each panel carries the same alphanumeric reference grid (columns A–M, rows 1–13) plus a north arrow, so any location can be named and compared across all four panels.

The four panels are:

Predicted land cover (top left) — the final ensemble class map, with the full class legend beneath it.
Predictive entropy (top right) — total uncertainty, on a 0–1 colour scale where dark purple is confident and bright yellow is maximally uncertain.
Expected entropy (bottom left) — data / aleatoric-like uncertainty, on the same 0–1 scale.
Mutual information (bottom right) — model / epistemic-like uncertainty, on the same 0–1 scale.

All three uncertainty panels deliberately share one colour scale and one normalization. That is what makes the single most important observation in this figure legible at a glance: the bottom-right panel is almost entirely dark while the two entropy panels are not. Since predictive entropy = expected entropy + mutual information, a near-black mutual-information panel means the total uncertainty is almost entirely explained by the aleatoric term. The five ensemble members essentially agree with each other; what they are collectively unsure about is the scene, not each other.

Reading the panels against one another also reveals structure that no single map shows:

The black patches in the entropy panels correspond exactly to water bodies and to the densest urban blocks in the class map — the two extremes of spectral distinctiveness.
The bright filamentary networks in the predictive-entropy panel trace the road network, the Neris river corridor, and forest–field boundaries. These are linear mixed-pixel features, one to two pixels wide, and they light up because a 10 m pixel straddling two land-cover types genuinely has no single correct label.
The diffuse bright regions in the north-west and around the outer suburbs coincide with the low-density residential and peri-urban zones in the class map, confirming that the fragmented housing-plus-garden landscape is the hardest part of this task.
Where mutual information does rise above the floor, it does so faintly at forest–grassland and cultivated–grassland boundaries. These are the same class pairs that dominate the off-diagonal cells of the confusion matrix, which is a useful consistency check: the models disagree with each other in exactly the places where they collectively disagree with the ground truth.

The uncertainty analysis

The uncertainty analysis gives a more detailed view of the final predictions than the class map alone.

The urbanized center of Vilnius and large built-up districts in the southern and south-eastern parts of the city are predicted quite confidently. Dense urban areas with apartment blocks, industrial buildings, and large artificial surfaces appear relatively stable in the uncertainty maps.

Water bodies are also predicted with high confidence. The Neris river channel, lakes, and ponds are clearly visible and have low uncertainty. Large forest and forest-park areas are also generally stable and confidently classified by the ensemble.

The most uncertain areas are not randomly distributed. They mostly appear where the land cover is naturally mixed or fragmented. Low-density residential areas are less homogeneous because private houses are often mixed with gardens, trees, grass, small roads, and other small objects. At Sentinel-2 resolution, such areas are difficult to separate into a single clean class.

Cultivated vegetation and natural grassland also show higher uncertainty. This is expected because these classes can be spectrally similar in summer composites, especially when fields, meadows, and unmanaged open land appear in similar phenological stages.

Areas along the river and major roads are among the brighter regions in the entropy maps. This likely happens for two reasons. First, these are linear transition zones where different land-cover types meet. Second, the surrounding areas may include mixed vegetation, wet or shadowed surfaces, shrubs, bare soil, and small artificial objects. As a result, the model sees several plausible classes rather than one obvious answer.

The comparison between expected entropy and mutual information is also useful. Expected entropy is much more visible than mutual information, while the mutual information map is mostly dark. This suggests that most uncertainty comes from the ambiguity of the land-cover signal itself rather than from strong disagreement between ensemble models.

In other words, the ensemble appears to be relatively consistent. The remaining uncertainty is mainly concentrated in areas where the segmentation task is physically difficult: boundaries between classes, mixed pixels, small objects, and fragmented suburban landscapes.

This uncertainty layer is useful because it changes how the final maps should be interpreted. Instead of treating every predicted pixel equally, it becomes possible to distinguish between confident predictions and areas where the model result should be read more carefully. In an operational setting, the entropy raster is the layer that tells an analyst which parts of the map can be used directly and which need manual review.

Running the pipeline

Requirements

pip install torch torchvision
pip install -r requirements.txt

torchgeo and segmentation_models_pytorch are imported lazily, so the custom Residual U-Net path runs without them installed.

1. Build the processed datasets

python scripts/prepare_landcovernet_median_dataset_multiband.py \
  --raw-dir /path/to/landcovernet_eu_2018 \
  --output-dir data/processed/landcovernet_10band \
  --band-preset 10bands \
  --months 6 7 8

python scripts/prepare_landcovernet_median_dataset_multiband.py \
  --raw-dir /path/to/landcovernet_eu_2018 \
  --output-dir data/processed/landcovernet_13band \
  --band-preset torchgeo13 \
  --months 6 7 8

Verify the result before training anything:

python scripts/check_processed_dataset.py \
  --data-dir data/processed/landcovernet_13band \
  --output-csv data/outputs/processed_dataset_check.csv

python scripts/inspect_processed_dataset.py \
  --data-dir data/processed/landcovernet_13band

2. Train

Custom Residual U-Net, best 10-band configuration (exp15):

python -m src.train \
  --data-root data/processed/landcovernet_10band \
  --output-dir outputs/exp15_resunet_f64 \
  --model resunet \
  --in-channels 10 \
  --base-features 64 \
  --random-crop-size 192 \
  --focal-weight 0.2 --dice-weight 0.8 \
  --normalization-mode reflectance \
  --epochs 30 --batch-size 8 --learning-rate 1e-3 \
  --early-stopping-patience 8

Sentinel-2 self-supervised model with discriminative fine-tuning (exp27):

python -m src.train \
  --data-root data/processed/landcovernet_13band \
  --output-dir outputs/exp27_torchgeo_dino_alpha \
  --model torchgeo_resnet50_unet \
  --in-channels 13 \
  --torchgeo-weights sentinel2_all_dino \
  --normalization-mode reflectance \
  --random-crop-size 192 \
  --focal-weight 0.2 --dice-weight 0.8 \
  --use-class-alpha \
  --learning-rate 1e-3 \
  --encoder-learning-rate 1e-4 \
  --epochs 30 --batch-size 8

3. Evaluate a single model

python -m src.evaluate \
  --data-root data/processed/landcovernet_13band \
  --checkpoint outputs/exp27_torchgeo_dino_alpha/best_model.pth \
  --output-dir outputs/exp27_torchgeo_dino_alpha/eval \
  --model torchgeo_resnet50_unet \
  --in-channels 13 \
  --normalization-mode reflectance \
  --use-class-alpha

The --seed, --model, --in-channels, and --normalization-mode arguments must match the training run, otherwise the evaluation split or the input scaling will differ from training.

4. Benchmark the ensemble on the test split

python scripts/ensemble_predict_test.py \
  --data-root-10 data/processed/landcovernet_10band \
  --data-root-13 data/processed/landcovernet_13band \
  --output-dir outputs/ensemble_test \
  --models exp22 exp27 exp19 exp21 exp20 \
  --exp22-checkpoint outputs/exp22/best_model.pth \
  --exp27-checkpoint outputs/exp27/best_model.pth \
  --exp19-checkpoint outputs/exp19/best_model.pth \
  --exp21-checkpoint outputs/exp21/best_model.pth \
  --exp20-checkpoint outputs/exp20/best_model.pth

5. Run inference on a real area of interest

python scripts/predict_real_ensemble_chips.py \
  --chips-root-10 data/aoi/vilnius_2025_10band \
  --chips-root-13 data/aoi/vilnius_2025_13band \
  --output-dir outputs/vilnius_2025 \
  --models exp22 exp27 exp19 exp21 exp20 \
  --exp22-checkpoint outputs/exp22/best_model.pth \
  --exp27-checkpoint outputs/exp27/best_model.pth \
  --exp19-checkpoint outputs/exp19/best_model.pth \
  --exp21-checkpoint outputs/exp21/best_model.pth \
  --exp20-checkpoint outputs/exp20/best_model.pth \
  --exclude-classes 3 \
  --save-uncertainty \
  --make-mosaic \
  --make-uncertainty-mosaics \
  --skip-existing

This produces the class-label mosaic plus predictive-entropy, expected-entropy, mutual-information, and confidence mosaics, all georeferenced and ready for GIS.

6. Regenerate the comparison figures

python scripts/visualize_model_comparison.py \
  --data-root-10 data/processed/landcovernet_10band \
  --data-root-13 data/processed/landcovernet_13band \
  --output-path data/figures/comparison_8_samples_with_best_ensemble_seed1707.png \
  --seed 1707 --num-samples 8 \
  --exp22-checkpoint outputs/exp22/best_model.pth \
  --exp27-checkpoint outputs/exp27/best_model.pth \
  --exp19-checkpoint outputs/exp19/best_model.pth \
  --exp21-checkpoint outputs/exp21/best_model.pth \
  --exp20-checkpoint outputs/exp20/best_model.pth

Limitations

Stated plainly, because a model card without them is marketing rather than engineering.

Natural Bare Ground is not usable. At IoU = 0.16 this class should be treated as unresolved. It covers 1.0% of training pixels and has no coherent spectral identity at 10 m.
Natural Grassland is only moderately reliable (IoU = 0.52, 63.7% recall), and it is systematically confused with woody and cultivated vegetation.
Temporal domain gap. The models are trained on 2018 labels and applied to 2015, 2020, and 2025 imagery. Radiometric differences between processing baselines are mitigated by the BOA offset correction, but sensor and atmospheric-correction changes across a decade are not fully removed.
Geographic domain gap. Training chips are distributed across Europe; inference is over a single Lithuanian city. Nothing in the training set is specific to Vilnius.
Absolute percentages are model estimates. The year-over-year trend is considerably more trustworthy than any individual year's absolute class share, because systematic model bias largely cancels when differencing two years produced by the same pipeline.
The uncertainty decomposition is a proxy. A five-member ensemble is not a posterior. Treating expected entropy as aleatoric and mutual information as epistemic is a widely used approximation whose limitations are discussed in the literature [15].
Scope of this repository. The tracked code covers dataset preparation, training, evaluation, ensembling, area-of-interest inference, and the model-comparison figures. The area-of-interest acquisition from Microsoft Planetary Computer and the cartographic poster rendering were run outside these scripts; predict_real_ensemble_chips.py consumes chips that have already been tiled and written to disk.

References

[1] Radiant Earth Foundation.
LandCoverNet. Source Cooperative / Radiant Earth dataset.
DOI: 10.34911/rdnt.63fxe5 | Source: LandCoverNet Europe

[2] Ronneberger, O., Fischer, P., & Brox, T. (2015).
U-Net: Convolutional Networks for Biomedical Image Segmentation. MICCAI 2015.
DOI: 10.1007/978-3-319-24574-4_28 | Preprint: arXiv:1505.04597

[3] He, K., Zhang, X., Ren, S., & Sun, J. (2016).
Deep Residual Learning for Image Recognition. CVPR 2016.
DOI: 10.1109/CVPR.2016.90 | Preprint: arXiv:1512.03385

[4] Chen, L.-C., Zhu, Y., Papandreou, G., Schroff, F., & Adam, H. (2018).
Encoder-Decoder with Atrous Separable Convolution for Semantic Image Segmentation. ECCV 2018.
DOI: 10.1007/978-3-030-01234-2_49 | Preprint: arXiv:1802.02611

[5] Zhou, Z., Rahman Siddiquee, M. M., Tajbakhsh, N., & Liang, J. (2018).
UNet++: A Nested U-Net Architecture for Medical Image Segmentation. DLMIA 2018.
DOI: 10.1007/978-3-030-00889-5_1 | Preprint: arXiv:1807.10165

[6] Tan, M., & Le, Q. V. (2019).
EfficientNet: Rethinking Model Scaling for Convolutional Neural Networks. ICML 2019.
Proceedings: PMLR v97 | Preprint: arXiv:1905.11946

[7] Lin, T.-Y., Dollár, P., Girshick, R., He, K., Hariharan, B., & Belongie, S. (2017).
Feature Pyramid Networks for Object Detection. CVPR 2017.
DOI: 10.1109/CVPR.2017.106 | Preprint: arXiv:1612.03144

[8] Stewart, A. J., Robinson, C., Corley, I. A., Ortiz, A., Lavista Ferres, J. M., & Banerjee, A. (2022).
TorchGeo: Deep Learning With Geospatial Data. ACM SIGSPATIAL 2022.
DOI: 10.1145/3557915.3560953 | Preprint: arXiv:2111.08872

[9] Lakshminarayanan, B., Pritzel, A., & Blundell, C. (2017).
Simple and Scalable Predictive Uncertainty Estimation using Deep Ensembles. NeurIPS 2017.
Proceedings: NeurIPS 2017 | Preprint: arXiv:1612.01474

[10] Smith, L., & Gal, Y. (2018).
Understanding Measures of Uncertainty for Adversarial Example Detection. UAI 2018.
Proceedings: UAI 2018 | Preprint: arXiv:1803.08533

[11] Malinin, A., & Gales, M. (2018).
Predictive Uncertainty Estimation via Prior Networks. NeurIPS 2018.
Proceedings: NeurIPS 2018 | Preprint: arXiv:1802.10501

[12] Depeweg, S., Hernández-Lobato, J. M., Doshi-Velez, F., & Udluft, S. (2018).
Decomposition of Uncertainty in Bayesian Deep Learning for Efficient and Risk-sensitive Learning. ICML 2018.
Proceedings: PMLR v80 | Preprint: arXiv:1710.11263

[13] Kendall, A., & Gal, Y. (2017).
What Uncertainties Do We Need in Bayesian Deep Learning for Computer Vision? NeurIPS 2017.
Proceedings: NeurIPS 2017 | Preprint: arXiv:1703.04977

[14] Houlsby, N., Huszár, F., Ghahramani, Z., & Lengyel, M. (2011).
Bayesian Active Learning for Classification and Preference Learning. arXiv.
DOI: 10.48550/arXiv.1112.5745 | Preprint: arXiv:1112.5745

[15] Wimmer, L., Sale, Y., Hofman, P., Bischl, B., & Hüllermeier, E. (2023).
Quantifying Aleatoric and Epistemic Uncertainty in Machine Learning: Are Conditional Entropy and Mutual Information Appropriate Measures? UAI 2023.
Proceedings: PMLR v216 | Preprint: arXiv:2307.03055

Released under the MIT License.

Contains modified Copernicus Sentinel-2 data. Sentinel-2 L2A imagery accessed via Microsoft Planetary Computer.

Name		Name	Last commit message	Last commit date
data		data
scripts		scripts
src		src
.gitignore		.gitignore
LICENSE		LICENSE
README.md		README.md
requirements.txt		requirements.txt

Latest commit

History

Folders and files

Repository files navigation