# 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
``````

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()
unzip("world.zip", exdir = tmpdir )
belgium <- world[world\$ADM0_A3 == "BEL", ]

centroids <- gCentroid(belgium, byid=TRUE)
spDists(centroids, longlat=TRUE)  #distance in km
``````
``````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)
{