I want to show data on a map of the Netherlands. To ensure my Dutch participants' privacy, I only asked them to specify the first two numbers of their postcode.

I found latitude and longitude for all Dutch postcodes at http://postcodedata.nl, but these are for the full postcode, and I don't know how to map these more 'spatially specific' latitudes and longitudes to my less fine-grained data.


Well, coincidentally I have been working on this issue myself recently :-)

Of course, Stack Exchange has the answer. A post at Stack Overflow, specifically, provides the algorithm to 'average' latitude and longitude data. For your convenience, and that of anybody who may search for this same issue, I have written an R script to achieve your goals:

### Compute 'central' longitude/latitude based on
###    https://stackoverflow.com/questions/6671183 and
###    http://www.geomidpoint.com/calculation.html

### Load ddply
if (!require(plyr)) {
  stop("The split-apply-combine system of plyr is required, but you don't ",
       "seem to have that package installed.");
} else {

  ### Read data (first download the zip file from http://postcodedata.nl and
  ### extract it; of course, don't forget to adjust the path here below.
  pcDat <- read.csv2("B:/Data/statistics/postcodes/postcode_NL_head.csv",
                     stringsAsFactors = FALSE);

  ### Extract the three relevant variables
  pcDat.2digits <- pcDat[, c('pnum', 'lat', 'lon')];

  ### Remove original (large) dataframe

  ### Remove last two numbers
  pcDat.2digits$pnum <- trunc(pcDat.2digits$pnum / 100);

  ### Convert lat and lon to radians
  pcDat.2digits$lat.radian <- as.numeric(pcDat.2digits$lat) * pi/180;
  pcDat.2digits$lon.radian <- as.numeric(pcDat.2digits$lon) * pi/180;

  ### Add X, Y and Z
  pcDat.2digits$x <- cos(pcDat.2digits$lat.radian) * cos(pcDat.2digits$lon.radian);
  pcDat.2digits$y <- cos(pcDat.2digits$lat.radian) * sin(pcDat.2digits$lon.radian);
  pcDat.2digits$z <- sin(pcDat.2digits$lat.radian);

  ### 'Average' coordinates
  pcDat <- ddply(.data = pcDat.2digits, .variables = 'pnum',
                 .fun = function(dat) {
    ### Average x, y and z
    meanX <- mean(dat$x);
    meanY <- mean(dat$y);
    meanZ <- mean(dat$z);
    ### Compute means
    meanLon.radian <- atan2(meanY, meanX);
    hyp <- sqrt(meanX^2 + meanY^2);
    meanLat.radian <- atan2(meanZ, hyp);
    ### Convert back to decimal system and store in dataframe
    res <- data.frame(pnum = mean(dat$pnum),
                      lat = meanLat.radian * 180/pi,
                      lon = meanLon.radian * 180/pi);
    ### Return result

  ### Load ggmap to show the results
  if (!require(ggmap)) {
    stop("To show the results, ggmap is required, but you don't ",
         "seem to have that package installed.");
  } else {

    ### Get map of the Netherlands from Google Maps
    netherlands <- get_map(location = c(lon = 5.2507557, lat = 52.1680345),
                           color = "color",
                           source = "google",
                           maptype = "roadmap",
                           zoom = 7)

    ### Map the resulting locations to the map as big red dots.
    ggmap(netherlands, extent='device') +
                 aes(x=lon, y=lat),

  ### You can now use merge() to combine these data with your dataframe.


Note: admins: I wasn't sure where to post this; Academia considered it off-topic, so this is my next attempt. If also considered off-topic here, I'll just relocate to Stack Overflow.

| improve this answer | |
  • What is the purpose of asking the question if you already know the answer? – Tim Jan 28 '16 at 15:37
  • 1
    @Tim There is no objection in principle to answering your own question. But I see this as (a) more geographical information science than statistics (b) too focused on code in specific software. It's not really a programming problem either. – Nick Cox Jan 28 '16 at 15:42
  • Sorry, didn't know this community existed! Posted just to share this code in case anybody else runs into this same problem (I've been looking for a solution for this for a few years now). Thank you very much for relocating @NickCox! – Matherion Jan 28 '16 at 16:43

For a country the size of the Netherlands you don't need to worry about the spherical nature of the Earth or the edge effects in that link. So a simple average will be fine.

You could actually use the rd_x and rd_y which I suspect are a projected coordinate system (I don't speak enough Dutch to read the metadata doc's to find out which one). This would allow you to no need to worry about the curvature of the Earth and you should be able to use other answers on this site to find out how to either map them directly or reproject them back to lat/lon if that is what your stats package requires.

If you are feeling particularly adventurous you could weight the postcode centroids by the count of properties in them (max-min?) to tip the centroid towards where people live.

| improve this answer | |
  • Good to know, thanks! The rd_x and rd_y data apparently relate to the "RD format", the "Rijksdriehoek systeem" (roughly means 'state triangle system'). I'm afraid I have no idea what that means . . . I'm afraid I'll pass on your suggestion for adventurous moods for now, this journey into what is apparently Geographic Information Systems was enough for today :-) – Matherion Jan 28 '16 at 16:46
  • Looks like RD is EPSG:28992 based on nl.wikipedia.org/wiki/Rijksdriehoeksco%C3%B6rdinaten#EPSG – Ian Turton Jan 28 '16 at 16:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.