I recently built a tile map of land use in Flanders, inspired by Bloomberg’s 2018 US land use graphic. Flanders is rendered as a grid of equal-area 2 km squares, each colored by its dominant land-use class, and scrolling through the page morphs the map into a waffle chart where block sizes match the true area of every class. A second view steps through 2013, 2016, 2022 and 2025 so you can watch the land change over twelve years.
This post is about the data pipeline that makes that possible. The end result is three JSON files that the frontend can load; everything before that is one Python script.
The source data
The Flemish government’s Departement Omgeving publishes the Landgebruik Vlaanderen raster series as open data: one GeoTIFF per snapshot year, with every 10 m cell in the Flemish region classified into one of 19 land-use classes. Each year has its own page in the Datavindplaats catalog — 2013, 2016, 2022 and 2025 — with direct GeoTIFF downloads. The 2013 and 2016 snapshots were reworked in 2023 to match the newer methodology, which is what makes the time series comparable at all.
For reference, Statistiek Vlaanderen publishes the official land-use numbers derived from the same dataset.
Tiling
The rasters are in Belgian Lambert 72 (EPSG:31370), a projected CRS, so a square in coordinates is a square in area. That makes the aggregation simple: cut the raster into 2 km × 2 km squares aligned to the raster’s upper-left corner, and count pixels per class inside each square.
Two refinements:
- Class grouping. The 19 source classes are too fine for a readable map, so they get mapped to 9 groups (cropland, grassland, built-up, forest, transport, nature, industry, recreation, water). The mapping lives in a config file with Dutch and English labels, so no class names are hardcoded in the pipeline.
- Edge tiles. Squares that are mostly outside Flanders would misrepresent the region’s shape, so any tile with less than 50% valid source area is dropped.
Each surviving tile is then assigned its dominant group: the one with the largest pixel count.
The largest-remainder correction
Dominant-class assignment has a known bias: a tile is indivisible, so small classes (water, recreation) get undercounted and large ones (cropland) get inflated. A tile that is 51% cropland and 49% water shows up as one more cropland block.
The fix is a largest-remainder correction:
- Compute each class’s true area in km² straight from the pixel counts.
- Convert that to a target tile count:
area_km2 / tile_area_km2, rounded with the largest-remainder method so the targets sum exactly to the number of retained tiles. - Reassign tiles from surplus classes to deficit classes, lowest confidence first — where confidence is the pixel margin between the dominant and the runner-up class. A tile that is 51/49 gets reassigned before one that is 90/10.
One extra trick: when building multiple years, candidates that already had the target class in the previous year are preferred. Without that anchor, borderline water tiles flicker between water and nature from one snapshot to the next, which reads as change that isn’t there.
In the 2025 build this reassigns roughly 700 of the 3,407 tiles — about one in five — which visibly changes the picture.
Transitions
For the evolution view, the pipeline also computes a transition matrix for every pair of years: for each 10 m cell, which class it had in year A and which in year B. From that the frontend can show “X km² of grassland became built-up between 2016 and 2022”, plus net gains and losses per class.
Output
The script writes one evolution.json with the tile grid, per-year dominant classes, true area shares and transition matrices, plus backward-compatible tiles.json and stats.json for the 2025 view. Every total shown in the UI comes from these files, not from hardcoded numbers.
The result is live here, and all the underlying data is free to download from the Flemish government. Highly recommended material if you want to build a scrollytelling map of your own.