问题描述:

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
网友答案:

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