# Calculate distance in kilometers between two points

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")

###
###

belgium.points = fortify(belgium, region="POSTCODE")
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
``````

``````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)
{
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

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()