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

[Roger Bivand]

New versions of spdep and spgwr fixing the same bug are now on CRAN.

Roger

On Thu, 8 Dec 2016, Edzer Pebesma wrote:

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

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: Roger.Bivand at nhh.no
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
http://depsy.org/person/434412