How the Klimaateffectatlas lookup works

Coordinate systems

PDOK returns both forms of the same point: centroide_ll is WGS84 and centroide_rd is RD New. The ArcGIS tab only needs centroide_ll, because ArcGIS converts it. The self-hosted tab uses centroide_rd, because its tiles are in RD New and the browser does no conversion either.

APIs we use

PDOK Locatieserver

Base URL: https://api.pdok.nl/bzk/locatieserver/search/v3_1

Step 1: find the address

Send the typed text to /suggest. fq=type:adres limits results to full addresses, and fl picks the fields to return.

const PDOK = 'https://api.pdok.nl/bzk/locatieserver/search/v3_1';

const q = 'Rokin 111 Amsterdam';
const resp = await fetch(`${PDOK}/suggest?q=${encodeURIComponent(q)}&fq=type:adres&rows=7&fl=id,weergavenaam,type,score`);
const docs = (await resp.json()).response.docs;

Response (trimmed):

{ "response": { "docs": [ {
  "type": "adres",
  "weergavenaam": "Rokin 111, 1012KN Amsterdam",
  "id": "adr-636952f00c72313cf8d7b88dcb4ca35f",
  "score": 10.274468
} ] } }

Step 2: get the postcode and coordinates

Pass the id to /lookup and ask for the fields you need.

const id = 'adr-636952f00c72313cf8d7b88dcb4ca35f';
const fields = 'weergavenaam,postcode,buurtcode,buurtnaam,wijkcode,wijknaam,gemeentenaam,provincienaam,centroide_ll,centroide_rd';
const resp = await fetch(`${PDOK}/lookup?id=${encodeURIComponent(id)}&fl=${fields}`);
const doc = (await resp.json()).response.docs[0];

Response (trimmed):

{
  "weergavenaam": "Rokin 111, 1012KN Amsterdam",
  "postcode": "1012KN",
  "buurtcode": "BU0363AE04",
  "buurtnaam": "Nes e.o.",
  "wijknaam": "Burgwallen-Oude Zijde",
  "gemeentenaam": "Amsterdam",
  "provincienaam": "Noord-Holland",
  "centroide_ll": "POINT(4.89278973 52.3695768)",
  "centroide_rd": "POINT(121328 486970)"
}

The points arrive as text, so pull the numbers out. Longitude comes first; in RD New, x comes first.

const match = doc.centroide_ll.match(/POINT\(([\d.]+)\s+([\d.]+)\)/);
const lon = parseFloat(match[1]);   // 4.89278973
const lat = parseFloat(match[2]);   // 52.3695768

const rd = doc.centroide_rd.match(/POINT\(([\d.]+)\s+([\d.]+)\)/);
const x = parseFloat(rd[1]);        // 121328 (only the self-hosted tab needs these)
const y = parseFloat(rd[2]);        // 486970

Converting the coordinates yourself

The lookup page never needs this, because ArcGIS converts for us. If you do need RD New metres from a GPS point (or the other way round) without calling a server, these two functions do it with plain arithmetic. They use the well-known approximation formulas by Schreutelkamp and Strang van Hees.

// WGS84 (GPS) -> RD New, approximation by Schreutelkamp & Strang van Hees
function wgs84ToRd(lon, lat) {
  const dLat = 0.36 * (lat - 52.15517440);
  const dLon = 0.36 * (lon - 5.38720621);
  // [power of dLat, power of dLon, coefficient]
  const X = [[0,1,190094.945],[1,1,-11832.228],[2,1,-114.221],[0,3,-32.391],[1,0,-0.705],
             [3,1,-2.340],[1,3,-0.608],[0,2,-0.008],[2,3,0.148]];
  const Y = [[1,0,309056.544],[0,2,3638.893],[2,0,73.077],[1,2,-157.984],[3,0,59.788],
             [0,1,0.433],[2,2,-6.439],[1,1,-0.032],[0,4,0.092],[1,4,-0.054]];
  const sum = (terms) => terms.reduce((t, [p, q, c]) => t + c * dLat ** p * dLon ** q, 0);
  return { x: 155000 + sum(X), y: 463000 + sum(Y) };
}

// RD New -> WGS84 (GPS)
function rdToWgs84(x, y) {
  const dX = (x - 155000) * 1e-5;
  const dY = (y - 463000) * 1e-5;
  // [power of dX, power of dY, coefficient]
  const LAT = [[0,1,3235.65389],[2,0,-32.58297],[0,2,-0.24750],[2,1,-0.84978],[0,3,-0.06550],
               [2,2,-0.01709],[1,0,-0.00738],[4,0,0.00530],[2,3,-0.00039],[4,1,0.00033],[1,1,-0.00012]];
  const LON = [[1,0,5260.52916],[1,1,105.94684],[1,2,2.45656],[3,0,-0.81885],[1,3,0.05594],
               [3,1,-0.05607],[0,1,0.01199],[3,2,-0.00256],[1,4,0.00128],[0,2,0.00022],
               [2,0,-0.00022],[5,0,0.00026]];
  const sum = (terms) => terms.reduce((t, [p, q, c]) => t + c * dX ** p * dY ** q, 0);
  return { lon: 5.38720621 + sum(LON) / 3600, lat: 52.15517440 + sum(LAT) / 3600 };
}

wgs84ToRd(4.89278973, 52.3695768);   // { x: 121327.82, y: 486969.77 }   PDOK says 121328, 486970
rdToWgs84(121328, 486970);           // { lon: 4.892792, lat: 52.369579 }

Checked against PDOK's own RD values for ten addresses from Groningen to Maastricht: the result was 0.30 m off every time, in both directions. That is fine for finding the right neighbourhood or postcode, but not for surveying work.

How the two versions compare

Checked for nine addresses from Groningen to Maastricht. Both tabs always found the same neighbourhood.

Some things are only in one tab:

• Only in ArcGIS: rainfall, rainfall deficit and tropical days per postcode (the Klimaatdata service), flood depth at the postcode, the number of severely lonely people aged 75 and over, and the detailed green breakdown in m². Klimaatdata is not a Klimaateffectatlas layer at all: Esri Nederland built it for klimaat.esri.nl by joining Klimaateffectatlas maps to postcodes. The Klimaateffectatlas itself only offers those maps as rasters (GeoTIFF on request via Data opvragen), so they are not in the self-hosted version yet.
• Only self-hosted: the share of the neighbourhood per flood-depth class, buildings with a dry floor, the feels-like temperature spread over the outdoor area, tree canopy, shade on all land, rental homes, WOZ value, and RIVM figures for loneliness, money worries and obesity among adults.