SteveG SteveG - 3 months ago 14
R Question

Difference in Great Circle calculation

I have some problems using the perimeter function from the geosphere package.

From first principles 1 deg of arc equals 1 nm
so 90deg of arc = 90 * 60 = 5400nm

When I use distGeo and two points on the equator I get an answer that corresponds to this figure.

library(geosphere)
distGeo(c(0,0),c(0,90), a=6378137, f=1/298.257223563)/1852
#[1] 5400.629


If I then try to do the same calculation with perimeter I get twice this distance!

#Two points on equator
xy <- rbind(c(0,0), c(90,0))
xy
# [,1] [,2]
#[1,] 0 0
#[2,] 90 0
perimeter(xy)/1852
#[1] 10819.39
perimeter(xy, a=6378137, f=1/298.257223563)/1852
#[1] 10819.39


If we take another point this time between the equator and the North Pole along the prime meridian the results are...

distGeo(c(0,0),c(90,0), a=6378137, f=1/298.257223563)/1852
#[1] 5409.694


This shows the effect of the ellipsoid.

When I use perimeter with these new points I still get the same problem as before!

#From the equator to the north Pole along the Prime meridian
xy <- rbind(c(0,0), c(0,90))
xy
# [,1] [,2]
#[1,] 0 0
#[2,] 0 90
perimeter(xy)/1852
#[1] 10801.26
perimeter(xy, a=6378137, f=1/298.257223563)/1852
#[1] 10801.26


So what am I doing wrong?

Steve

Answer

You are doing nothing wrong. The function perimeter computes the perimeter distance around a closed polygon, and it will close that polygon if the end point is not the first, so it is interpreting your two points as being a polygon with no "height".

To illustrate (results in meters to see small additions):

## your results
xy <- rbind(c(0,0), c(90,0))
perimeter(xy)
##[1] 20037508

## "closing the polygon" with c(0,0) gives same answer
xy <- rbind(c(0,0), c(90,0), c(0,0)) 
perimeter(xy)
##[1] 20037508

## adding one meter (0.00001 deg latitude) north as another point, 
## perimeter will close the resulting triangle, result increases by one meter!
xy <- rbind(c(0,0), c(90,0), c(90, 0.00001))
perimeter(xy)
##[1] 20037509

## closing the triangle yourself, same answer
xy <- rbind(c(0,0), c(90,0), c(90,0.00001), c(0,0)) 
perimeter(xy)
##[1] 20037509