SteveG - 9 months ago 43

R Question

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

Source (Stackoverflow)