Laying an airport on real ground with the IGN terrain model

The first version of Chronotopia was flat. Everything sat on a plane at height zero, which works for the runways at El Prat, on the Llobregat delta, and badly for everything around them. At Ibiza it was worse: the 1972 Iberia crash on the sa Talaiassa ridge had a comment in the code saying “the model has no terrain, so the ridge is not drawn”.

Commit b637751 added real ground to El Prat and Ibiza. The heights come from the Instituto Geográfico Nacional’s MDT25 terrain model, a build script turns them into a 40 m grid 34 km across with the airport levelled in the middle, and the viewer lays every map layer, building and approach-light mast on it. This article follows the heights from the IGN’s server to the screen, and ends with what the result gets wrong.

Where the heights come from

The MDT25 is a bare-earth model of Spain: ground height without buildings or trees, on a grid of about 25 m. The IGN serves it through a WCS, a web service that returns a piece of the raster for a bounding box. scripts/fetch-terrain.mjs builds the request by hand: coverage Elevacion4258_25 (ETRS89 latitude and longitude) as a GeoTIFF, for the grid plus a 0.01° margin, cached until the bounding box changes. The downloads for El Prat and Ibiza are about 5.3 to 5.5 MB each.

The TIFF reader is written by hand too, about 40 lines. It accepts exactly one kind of file, uncompressed single-band 16-bit data in strips, and throws on anything else. That is what the WCS returns, and reading a handful of tags directly was simpler than adding a GeoTIFF library for one format. Because the grid is in degrees, a pixel is 25 m north–south but only about 19 m east–west at these latitudes.

The service returns whole metres, so on gentle slopes the raw data has one-metre terraces. The script smooths it with a small Gaussian blur (sigma 0.9 pixels) and then samples it bilinearly onto the output grid: 851 × 851 points, 40 m apart, 34 km on each side. The grid is centred per site so that more of it falls on land. At El Prat it reaches Collserola, Montjuïc and the Garraf; at Ibiza it covers the south of the island, from sa Talaiassa to Dalt Vila.

Positions use the same flat projection as the rest of the app: metres east and south of the airport’s origin, with a degree of longitude worth 111,320 m × cos(latitude). Over 34 km that is accurate enough, and much simpler than the projections we use for whole continents in Globdrop.

Levelling the airport

A real airfield isn’t flat. The cells under Ibiza’s runway, aprons and terminal range from 1.1 to 16.9 m above sea level. Drawn as they are, runways would bend and aircraft would taxi uphill, so the script flattens the airport and blends it back into the ground around it.

First it marks the platform: every apron, terminal, hangar and tower polygon, the runways plus 90 m on each side of the centreline, and the taxiways plus 20 m. Then it takes the median height of those cells as the platform level and ramps back to natural ground over 220 m. From scripts/fetch-terrain.mjs:

const levels = [];
for (let k = 0; k < N; k++) if (platform[k] && !sea[k]) levels.push(hgt[k]);
levels.sort((a, b) => a - b);
const base = opts.base ?? Math.round((levels[Math.floor(levels.length / 2)] ?? 0) * 10) / 10;
const toPlatform = distanceTo(platform);
const smoothstep = (t) => (t <= 0 ? 0 : t >= 1 ? 1 : t * t * (3 - 2 * t));
for (let k = 0; k < N; k++) {
  if (sea[k]) continue;
  const w = smoothstep(toPlatform[k] / opts.ramp);
  hgt[k] = base + (hgt[k] - base) * w;
}

distanceTo is a two-pass distance transform over the grid, so w is 0 on the platform and reaches 1 at 220 m from its nearest cell. The median rather than the mean keeps a few odd cells, a cutting or an embankment, from pulling the level. The platform becomes the model’s height zero, and everything else is stored relative to it.

At El Prat the levelling barely changes anything: the platform is at 2.5 m, over 5,402 cells that range from 0.1 to 7.0 m. At Ibiza it matters. The platform ends up at 6.3 m over 1,244 cells, so the higher north-west edge of the airfield is cut down by up to 10 m and the lower south-east side is raised by up to 5 m.

Where the land meets the sea

The sea is a single opaque plane at sea level, 200 km across. Wherever the ground mesh is below it, you see water, so the coastline is wherever the mesh crosses zero. Left alone, that line would follow the terrain model’s coast, which doesn’t match the OpenStreetMap coastline the rest of the map is drawn from, and beaches and breakwaters would end up half under water or floating in front of it.

The script fixes this in three steps. Land is never lower than 0.3 m. Sea cells get an invented depth that grows with the distance from land, 1.5 m plus 2 cm per metre, down to 12 m; nobody sees it, because the plane is opaque. And sea cells within 60 m of the OSM coastline get a depth chosen so that the line between them and their land neighbour crosses zero exactly on the mapped coast:

if (sea[k]) {
  // Slope of the neighbouring shore, so the crossing lands at the right fraction of the cell.
  const i = k % n;
  let shore = 0.6;
  for (const o of [-1, 1, -n, n]) {
    const q = k + o;
    if (q < 0 || q >= N || (o === -1 && i === 0) || (o === 1 && i === n - 1) || sea[q]) continue;
    shore = Math.max(shore, hgt[q] / Math.max(S * 0.25, near[q] === Infinity ? S : near[q]));
  }
  hgt[k] = Math.min(hgt[k], -Math.max(0.3, near[k] * Math.min(shore, 0.5)));
}

near[k] is the cell’s distance to the coastline, so a sea point 20 m offshore goes down by 20 m times the slope of the land beside it. The rendered surface drops through zero right where OSM says the water starts.

Standing things on the ground

The file stores heights in decimetres above the platform, as the difference from the previous value, compressed with zlib: a JSON header, then 724,201 numbers. El Prat’s is 585,828 bytes and Ibiza’s 388,513, smaller because more of its grid is flat sea. The header carries the attribution, “MDT25 © Instituto Geográfico Nacional (CC BY 4.0)”. A site without the file is drawn flat, as before.

In the browser, the grid becomes one mesh, with each 40 m cell split into two triangles along the same diagonal. Everything else has to sit on those triangles, not on a smoother guess of the ground, or objects would float or sink by the difference. So height() in src/engine/terrain.ts doesn’t interpolate bilinearly; it finds which of the cell’s two triangles the point is in and interpolates on that one:

height(x: number, z: number): number {
  const n = this.n;
  const gu = Math.min(n - 1, Math.max(0, (x - this.x0) / this.step));
  const gv = Math.min(n - 1, Math.max(0, (z - this.z0) / this.step));
  const i = Math.min(n - 2, Math.floor(gu));
  const j = Math.min(n - 2, Math.floor(gv));
  const fu = gu - i;
  const fv = gv - j;
  const h = this.h;
  const a = h[j * n + i];
  const d = h[(j + 1) * n + i + 1];
  if (fu >= fv) {
    const b = h[j * n + i + 1];
    return a + fu * (b - a) + fv * (d - b);
  }
  const c = h[(j + 1) * n + i];
  return a + fv * (c - a) + fu * (d - c);
}

Flat map layers, such as fields, roads, water and rail, are harder, because a large polygon spans many cells. drape() cuts each of its triangles along the grid columns, then the rows, then each cell’s diagonal, so that every piece lies inside a single ground triangle, and lifts the pieces with height(). Where the ground under a triangle is level to within 2 cm, which covers the whole airport platform, it skips the cutting and just lifts it. The mesh itself is drawn 12 cm below the true height, so draped layers never flicker against it.

Solid things each have a rule:

• Buildings stand on the lowest ground point under their footprint’s corners, so on a slope the uphill side is buried rather than the downhill side floating. On the platform that’s zero.
• Trees are sunk 40 cm, so trunks on slopes don’t hang in the air on their downhill side.
• Viaducts keep their grade across a valley, and piers go down to the ground every 28 m or so wherever the deck is more than 3 m above it.
• Approach-light masts grow taller as the ground drops away beneath the approach, and are never less than 2.5 m above the ground under their bar.

On phones the viewer uses every second grid point, 80 m apart, which cuts the ground from about 1.45 million triangles to 360,000.

What the relief gets wrong

• It’s bare earth. The MDT25 has no buildings or trees; those come from OpenStreetMap and stand on it.
• It smooths the peaks. Between the source resolution, the blur and the 40 m grid, ridges lose a metre or more. Sa Talaiassa, 475 m in the raw data, is 473.9 m in the model. Anything narrower than about 80 m, such as an embankment or a river bank, is smeared.
• The airport is one plane. Half of Ibiza’s airfield is raised and half lowered to the median, and the real 16 m of slope across it is gone.
• The salt pans are land. Ground at or below sea level is lifted to 0.3 m, so ses Salines at Ibiza renders as flat land just above the water, with the mapped ponds draped on top.
• The ground doesn’t change with time. A 1950 runway is drawn on today’s relief. On the flat delta that hardly shows, but ground that was reshaped later, such as the Llobregat’s new channel, is drawn as it is now in every year.
• The seabed is invented. Depths depend only on the distance to land. That’s fine under an opaque sea, and wrong for anything that would need to see through it.

The levelled platform also keeps the runways simple: the viewer draws them at height zero without asking the terrain, because by construction that is where the ground is. How that strip grows and shrinks over the years is the subject of a runway is one line and two fractions.