While building support for geosphere distance functions in sf and comparing results from geosphere::distMeeuws with those of sp::spDists, I found a typo in the C code underlying the great circle distance computation functions in sp (sp::spDists, sp::spDistsN1) [*]:
https://github.com/edzer/sp/commit/d8374ff7efc6735cba9a054748c602bed0672f23 Before: > library(sf) Linking to GEOS 3.5.0, GDAL 2.1.0, proj.4 4.9.2 > library(sp) > library(units) > > x = st_sfc( + st_point(c(0,0)), + st_point(c(1,0)), + st_point(c(2,0)), + st_point(c(3,0)), + crs = 4326 + ) > > y = st_sfc( + st_point(c(0,10)), + st_point(c(1,0)), + st_point(c(2,0)), + st_point(c(3,0)), + st_point(c(4,0)), + crs = 4326 + ) > > st_distance(x, y) Units: m [,1] [,2] [,3] [,4] [,5] [1,] 1105855 111319.5 222639.0 333958.5 445278.0 [2,] 1111387 0.0 111319.5 222639.0 333958.5 [3,] 1127822 111319.5 0.0 111319.5 222639.0 [4,] 1154693 222639.0 111319.5 0.0 111319.5 > > d.sf = st_distance(x, y, geosphere::distMeeus) > d.sp = spDists(as(x, "Spatial"), as(y, "Spatial")) > units(d.sp) = make_unit("km") > d.sf - d.sp Units: m [,1] [,2] [,3] [,4] [,5] [1,] 1851.990 0.000000e+00 0.000000e+00 0.000000e+00 0 [2,] 1842.938 0.000000e+00 0.000000e+00 2.910383e-11 0 [3,] 1816.564 0.000000e+00 0.000000e+00 1.455192e-11 0 [4,] 1775.032 2.910383e-11 1.455192e-11 0.000000e+00 0 > After: > library(sf) Linking to GEOS 3.5.0, GDAL 2.1.0, proj.4 4.9.2 > library(sp) > library(units) > > x = st_sfc( + st_point(c(0,0)), + st_point(c(1,0)), + st_point(c(2,0)), + st_point(c(3,0)), + crs = 4326 + ) > > y = st_sfc( + st_point(c(0,10)), + st_point(c(1,0)), + st_point(c(2,0)), + st_point(c(3,0)), + st_point(c(4,0)), + crs = 4326 + ) > > st_distance(x, y) Units: m [,1] [,2] [,3] [,4] [,5] [1,] 1105855 111319.5 222639.0 333958.5 445278.0 [2,] 1111387 0.0 111319.5 222639.0 333958.5 [3,] 1127822 111319.5 0.0 111319.5 222639.0 [4,] 1154693 222639.0 111319.5 0.0 111319.5 > > d.sf = st_distance(x, y, geosphere::distMeeus) > d.sp = spDists(as(x, "Spatial"), as(y, "Spatial")) > units(d.sp) = make_unit("km") > d.sf - d.sp Units: m [,1] [,2] [,3] [,4] [,5] [1,] -2.328306e-10 0.000000e+00 0.000000e+00 0.000000e+00 0 [2,] 2.328306e-10 0.000000e+00 0.000000e+00 2.910383e-11 0 [3,] 0.000000e+00 0.000000e+00 0.000000e+00 1.455192e-11 0 [4,] 0.000000e+00 2.910383e-11 1.455192e-11 0.000000e+00 0 > [*] assuming the source code of http://www.abecedarical.com/javascript/script_greatcircle.html is correct; I don't have access to the Meeuws book mentioned there under (1). The same C function is used by gstat; I corrected that as well. -- Edzer Pebesma Institute for Geoinformatics (ifgi), University of Münster Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081 Journal of Statistical Software: http://www.jstatsoft.org/ Computers & Geosciences: http://elsevier.com/locate/cageo/
library(sf)
library(sp)
library(units)
x = st_sfc(
st_point(c(0,0)),
st_point(c(1,0)),
st_point(c(2,0)),
st_point(c(3,0)),
crs = 4326
)
y = st_sfc(
st_point(c(0,10)),
st_point(c(1,0)),
st_point(c(2,0)),
st_point(c(3,0)),
st_point(c(4,0)),
crs = 4326
)
st_distance(x, y)
d.sf = st_distance(x, y, geosphere::distMeeus)
d.sp = spDists(as(x, "Spatial"), as(y, "Spatial"))
units(d.sp) = make_unit("km")
d.sf - d.sp
#summary(unclass(d.sf) - d.sp)
st_crs(x) = st_crs(y) = NA
d.sf = st_distance(x, y)
d.sp = spDists(as(x, "Spatial"), as(y, "Spatial"))
d.sf - d.sp
signature.asc
Description: OpenPGP digital signature
_______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
