Confusion regarding distance calculation in R (euclidean distance, "great circle distance")

Question

OK so considering these two cases:

ln1<-SpatialLinesDataFrame(SpatialLines(list(Lines(Line(matrix(c(53.3604464,53.36062,-6.2424442, -6.242413),ncol=2)),ID="a"))),data=data.frame(dummy="a"),match.ID=F)
proj4string(pt1) <- CRS("+init=epsg:4326")

SpatialLinesLengths(ln1,longlat=T)*1000
SpatialLinesLengths(spTransform(ln1, CRS("+init=epsg:3857")),longlat=F)


ln2<-SpatialLinesDataFrame(SpatialLines(list(Lines(Line(matrix(c(15.43911,15.43914,47.00849, 47.00837),ncol=2)),ID="a"))),data=data.frame(dummy="a"),match.ID=F)
proj4string(ln2) <- CRS("+init=epsg:4326")

SpatialLinesLengths(ln2,longlat=T)*1000
SpatialLinesLengths(spTransform(ln2, CRS("+init=epsg:3857")),longlat=F)

I calculate the lengths of the lines (ln1 and ln2) in meters.

The first calculation being the "great circle distance" the second Euclidean distance. Well I read that those distances should lie pretty close to each other when calculated for small distances. That is true for the first case:

Great Circle:

SpatialLinesLengths(ln1,longlat=T)*1000
[1] 19.51758

Euclidean

SpatialLinesLengths(spTransform(ln1, CRS("+init=epsg:3857")),longlat=F)
[1] 19.63836  

But in the second case the length difference is pretty great. I mean its over 40%...

Great Circle:

SpatialLinesLengths(ln2,longlat=T)*1000
[1] 13.52404

Euclidean

SpatialLinesLengths(spTransform(ln2, CRS("+init=epsg:3857")),longlat=F)
[1] 19.87276 

Well I understand the difference between both methods (straight line vs. "as the crow flies" etc.) but reading (and understanding so) that the difference on small scale should not be to big. I worry seeing something like that...

Is it just because of the distance to the Equator? (What I can't imagine) Is it a rounding issue? Is my code wrong? (Well the same effect takes place using gLength(rgeos) or spDists/spDistsN1(sp) or any other distance calculation out there for R)

So what's going on here?

Answer 1

EPSG 3857 is a Mercator projection, and this is not a suitable choice for distance calculations, generally.

This can be a subtle and tricky topic, and your example brings out an problem very clearly. To be as accurate as ellipsoidal calculations generally, i.e. for any particular distance between two arbitrary points you must choose an equidistant projection that has one of the points as its centre. A reasonable choice is Azimuthal Equidistant, and in R you can achieve it like this:

SpatialLinesLengths(spTransform(pt2, CRS(sprintf("+proj=aeqd +lon_0=%f +lat_0=%f +ellps=WGS84 +no_defs", m2[1, 1], m2[1, 2]))))
## [1] 13.53418

Note that this is a bit impractical for a large set of points since you have to create a local projection for every single one. The complication is that you need a local projection for particular sets of point pairs, and so grouping your points into optimal sets gets tricky. You can choose a local projection for distance measures where the properties will be sufficient, but it's specific to the locality obviously.

This is not true for some equal area projections, they really are equal across the entire range available. But you can still get caught by topological limitations, if your lines/edges don't traverse a curved space properly, also true for distance calcs.

Generally, for any arbitrary points in the world you need to use ellipsoid methods, and these are effectively vectorized in R in a way that local projections for pairs of points cannot.

Answer 2

library(rgdal) library(sp)

COORDINATES #DATAFRAME THAT CONTENT THE LATITUDE (LAT) AND LONGITUDE (LON) IN THE COORDINATE REFERENT SYSTEM (CRS) WGS84.

coordinates(COORDINATES) <- ~ LON + LAT

proj4string(COORDINATES) <- CRS("+proj=longlat +datum=WGS84") #ASSIGN THE CRS

Zone <- input$Zone #UTM ZONE FOR YOUR COUNTRY

COORDINATES <- spTransform(COORDINATES, CRS(paste("+proj=utm", " +zone=", Zone, " +ellps=WGS84", " +datum=WGS84", " +units=m", sep=""))) #REPROJECT THE CRS

COORDINATES <- as.data.frame(COORDINATES)
X <- COORDINATES$LON  #EXTRACT THE LOGITUDE VECTOR
Y <- COORDINATES$LAT  #EXTRACT THE LATITUDE VECTOR
MX1 <- X %*% t(X) #CREATE A MATRIX FOR LONGITUDE VECTOR
MX2 <- matrix(rep(t(X),nrow(COORDINATES)), ncol = nrow(COORDINATES), nrow = nrow(COORDINATES)) #CREATE A MATRIX FOR REPEAT LONGITUDE VECTOR
MX <- MX1/MX2 #DEFENITIVE MATRIX FOR LONGITUDE VECTORS
MX <- abs((MX-MX2)**2) #SQUARE SUM OF LONGITUDE VECTORS
colnames(MX)<- paste(COORDINATES$STATION) #ASSIGN COLNAMES
rownames(MX)<- paste(COORDINATES$STATION) #ASSIGN ROWNAMES
MY1 <- Y %*% t(Y) #CREATE A MATRIX FOR LATITUDE VECTOR
MY2 <- matrix(rep(t(Y), nrow(COORDINATES)), ncol = nrow(COORDINATES), nrow = nrow(COORDINATES)) #CREATE A MATRIX FOR REPEAT LATITUDE VECTOR
MY <- MY1/MY2 #DEFENITIVE MATRIX FOR LATITUDE VECTORS
MY <- abs((MY-MY2)*2) #SQUARE SUM OF LONGITUDE VECTORS
colnames(MY)<- paste(COORDINATES$STATION) #ASSIGN COLNAMES
rownames(MY)<- paste(COORDINATES$STATION) #ASSIGN ROWNAMES
EUCLIDEAND <- round((sqrt(MX+MY)/1000), digits = 0) #EUCLIDEAN DISTANCE FOR THESE COORDINATES
EUCLIDEAND <- as.data.frame(EUCLIDEAND)