Calculate distance in kilometers between two points

Question

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

Answer 1

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

Answer 2

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