[GIS] Calculate distance in kilometers between two points

distancershapefile

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

Best Answer

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
Related Question