4

I have a few questions regarding the calculation of the distance between 2 postal codes. Take a look at the following code.

###
# load the different packages we need
###

require("rgdal")
require("maptools")
require("ggplot2")
require("plyr")
require("rgeos")

###
# load the shape file
###

belgium = readOGR(dsn="XXX", layer="NAME")
belgium.points = fortify(belgium, region="POSTCODE") 
head(belgium.points)
tail(belgium.points)

This gives ;

      long      lat     id
1 151584.6 175068.7    1000
2 150829.0 174739.9    1000
3 150614.4 174580.1    1000
4 150178.8 174208.4    1000
5 149715.7 173823.0    1000
6 149126.5 173269.7    1000

      long      lat    id
51839 83871.46 215487.5 9992
51840 84519.46 215217.5 9992
51841 83925.51 214823.3 9992
51842 83578.20 214704.9 9992
51843 82796.53 215022.7 9992
51844 81811.86 214716.5 9992

What I now want to do is calculating the distance between the centroid of each postal code. The postal code is given by the variable id in the dataset belgium.points. Is the code below valid? Can I convert it to kilometers? (Please if you need additional info just ask me.

lat = sapply(split(belgium.points, belgium.points$id), function(x) mean(x$lat))/1000
long = sapply(split(belgium.points, belgium.points$id), function(x) mean(x$long))/1000
lat.dist = outer(lat,lat,function(lat1,lat2) (lat1-lat2)^2)
long.dist = outer(long,long,function(long1,long2) (long1-long2)^2)
spatial.dist = lat.dist + long.dist
1
lat = sapply(split(belgium.points, belgium.points$id), function(x) mean(x$lat))/1000

... You sure this gives you a correct latitude value? If you know the initial projection, then why not convert it into a spatialpointdataframe and use spproject with a lat-long projection to get the actual values.

To calculate distances:

  1. Ether use the PointDistance tool in the raster package. It even works with non-latlon point locations
  2. Otherwise use this custom R function to approximate the distance

    # Calculate distance in kilometers between two points
    earth.dist <- function (long1, lat1, long2, lat2)
    {
    rad <- pi/180
    a1 <- lat1 * rad
    a2 <- long1 * rad
    b1 <- lat2 * rad
    b2 <- long2 * rad
    dlon <- b2 - a2
    dlat <- b1 - a1
    a <- (sin(dlat/2))^2 + cos(a1) * cos(b1) * (sin(dlon/2))^2
    c <- 2 * atan2(sqrt(a), sqrt(1 - a))
    R <- 6378.145
    d <- R * c
    return(d)
    }
    
  • I used that particular function to calculate the centroid of my postcode areas. For my application it don't has to be a very exact result. – Roger Dec 10 '13 at 17:30
  • I found out that the latitude/longitude values were in fact not expressed in degrees but in meters. I noticed that the difference in 1 degree longitude was more or less 111200 meters. The difference in latitude was 70000 à 74000 meters (dependent on the longitude). I'll mark you answer as correct though, because it answers the original question in my opinion. – Roger Dec 10 '13 at 17:35
8

Use gCentroid from rgeos to calculate centroids. Use spDists from sp (as mentioned by @Jot eN above) to calculate distance matrices.

I made up a simple example using admin boundaries. Note that the parameters of spDists change, depending on your map units. The example is for degrees. If you have meters use something like: spDists(centroids)/1000

library(sp)
library(rgdal)
library(rgeos)

# getting sample data of the three main Belgian regions
tmpdir <- tempdir()
download.file("http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_0_map_subunits.zip", "world.zip")
unzip("world.zip", exdir = tmpdir )
shapefile <- paste0(tmpdir,"/ne_10m_admin_0_map_subunits.shp")
world <- readOGR(shapefile, "ne_10m_admin_0_map_subunits")
unlink(tmpdir)
belgium <- world[world$ADM0_A3 == "BEL", ]

centroids <- gCentroid(belgium, byid=TRUE)
spDists(centroids, longlat=TRUE)  #distance in km

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.