1

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

1

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.

2

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.