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# 5400.629``

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

``#Two points on equatorxy <- rbind(c(0,0), c(90,0))xy# [,1] [,2]#[1,] 0 0#[2,] 90 0perimeter(xy)/1852# 10819.39perimeter(xy, a=6378137, f=1/298.257223563)/1852# 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# 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 meridianxy <- rbind(c(0,0), c(0,90))xy# [,1] [,2]#[1,] 0 0#[2,] 0 90perimeter(xy)/1852# 10801.26perimeter(xy, a=6378137, f=1/298.257223563)/1852# 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)
## 20037508

## "closing the polygon" with c(0,0) gives same answer
xy <- rbind(c(0,0), c(90,0), c(0,0))
perimeter(xy)
## 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)
## 20037509

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

Thanks for your guidance.

What I ended-up doing was to lag the lat and long fields

``````xy[,3] = lag(xy[,1],k=-1)
xy[,4] = lag(xy[,2],k=-1)
``````

Then compute the segment using distGEO

``````xy\$dDist<-distGeo(as.vector(xy[,1:2]),
as.vector(xy[,3:4]),
a=6378137,
f=1/298.257223563)/1852
``````

Then take the sum

``````sum(xy\$dDist,na.rm = TRUE)
``````

It seems to work!

Top