Compute a distance matrix in R with geosphere

Question

I am trying to compute the a distance matrix between two sets of points in two shapefiles (S1 and S2) in R. I am new to geographical analysis in R, having previously only used QGIS and ArcMap.

I am using to do so the packages geosphere and rgdal. My two shapefiles are in the projection EPSG:102023. When I try to compute the distances using distm, I get the error:

 Error in .pointsToMatric(x): longitude <-360

Does this mean that you have to use Geographic coordinate systems to use distm?

Here is my code:

proj4string(S1)

[1] "+proj=eqdc +lat_0=0 +lon_0=25 +lat_1=20 +lat_2=-23 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

head(coordinates(S1))

  coords.x1 coords.x2

[1,] 330091.20 -3277120

[2,] -672066.99 -3753558

[3,] -671040.65 -3755774

[4,] -676182.43 -3741369

[5,] 61188.02 -3759233

[6,] 156078.69 -3687183

proj4string(S2)

[1] "+proj=eqdc +lat_0=0 +lon_0=25 +lat_1=20 +lat_2=-23 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

head(coordinates(S2))

  coords.x1  coords.x2

[1,] -1246820.3 -586142.4

[2,] 1831068.8 1024375.2

[3,] -1093944.1 -617261.2

[4,] -1075037.3 -614924.2

[5,] -652579.9 -3745182.5

[6,] -498183.2 261947.4

test<-distm(coordinates(S1),coordinates(S2), fun=distHaversine)

Error in .pointsToMatrix(x) : longitude < -360

Answer 1

From documentation:

Arguments:

x longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object

y Same as x. If missing, y is the same as x

You need to use lat/long CRS. So:

library(sp)
S1b <- spTransform(S1, CRS("+init=epsg:4326"))
S2b <- spTransform(S2, CRS("+init=epsg:4326"))
test <- distm(coordinates(S1b),coordinates(S2b), fun=distHaversine)

Answer 2

By passing coordinates(S1) you are dropping the coordinate reference system and passing a simple matrix. distm then assumes you have lat-long and tries to do the calculation, but discovers your lat-long are huge values, outside the range of -360 to 360 (and -90 to +90) so you get that error.

If you pass the spatial objects themselves you get a different error:

> test<-distm(s1,s2, fun=distHaversine)
Error in .pointsToMatrix(x) : 
  Points are projected. They should be in degrees (longitude/latitude)

which should make it clearer you need to spTransform to degrees lat-long before using distm.