SteveG - 1 month ago 4x
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

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
``````
Source (Stackoverflow)