[R-sig-Geo] sp::spDists*

[Edzer Pebesma]

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/
-------------- next part --------------
A non-text attachment was scrubbed...
Name: dist.R
Type: text/x-r-source
Size: 562 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20161208/a3918d17/attachment.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 473 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20161208/a3918d17/attachment-0001.bin>