In acos(cos(dp)/cos(xtr)) : NaNs produced

[amael-ls]

I ran into a strange bug that can be reproduced with R/3.6.0, but not with R/3.5.1, although the packages are at the same version. Here is the reproducible example (zip file at the end of the message):

library(geosphere)
library(sp)

pt = readRDS("./bug_coords.rds")
line = readRDS("./line.rds")

R.Version()$version.string
packageVersion("sp")
packageVersion("geosphere")

dist2Line(pt, line, distfun = geosphere::distGeo)

plot(line)
points(pt["x"], pt["y"], pch = 15, col = "red")

I got the following error on R/3.6.0:

> R.Version()$version.string
[1] "R version 3.6.0 (2019-04-26)"
> packageVersion("sp")
[1] ‘1.3.1’
> packageVersion("geosphere")
[1] ‘1.5.10’
>
> dist2Line(pt, line, distfun = geosphere::distGeo)
Error in if (trackdist1[j] < trackdist2[j]) { :
  missing value where TRUE/FALSE needed
In addition: Warning message:
In acos(cos(dp)/cos(xtr)) : NaNs produced

while I got the following on R/3.5.1

> R.Version()$version.string
[1] "R version 3.5.1 (2018-07-02)"
> packageVersion("sp")
[1] ‘1.3.1’
> packageVersion("geosphere")
[1] ‘1.5.10’
>
> dist2Line(pt, line, distfun = geosphere::distGeo)
     distance       lon      lat ID
[1,] 239791.1 -73.52718 55.60157  1

bug_sp.zip

[rhijmans]

It do not get that bug on R 3.6.0 nor on current R (3.6.3).

library(geosphere)
library(sp)

pt = readRDS("./bug_coords.rds")
line = readRDS("./line.rds")

R.Version()$version.string
#[1] "R version 3.6.3 (2020-02-29)"
packageVersion("geosphere")
#[1] ‘1.5.10’

dist2Line(pt, line, distfun = geosphere::distGeo)
#     distance       lon      lat ID
#[1,] 239791.1 -73.52718 55.60157  1

plot(line)
points(pt["x"], pt["y"], pch = 15, col = "red")

Maybe there are packages loaded? Did you completely clean your workspaces and compare your sessionInfo()'s?

[amael-ls]

The bug occurs on a supercomputer, but the sessionInfo() of my own laptop is quite similar. Do you think it could be related to the version of gdal it is linked?

Session info on the super computer:

> sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /cvmfs/soft.computecanada.ca/easybuild/software/2017/Core/imkl/2018.3.222/compilers_and_libraries_2018.3.222/linux/mkl
/lib/intel64_lin/libmkl_gf_lp64.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8

 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] sp_1.3-1         geosphere_1.5-10

loaded via a namespace (and not attached):
[1] compiler_3.6.0  grid_3.6.0      lattice_0.20-38

[amael-ls]

The only notable difference with the session of my laptop is that my laptop also loads tools_3.5.1

[amael-ls]

I decided to contact the supercomputer help centre, I will keep you posted. I could not reproduce the problem on R/3.6.0 at home either...

[amael-ls]

The error comes from alongTrackDistance <- function(p1, p2, p3, r=6378137), line 28:

dist <- bearing * acos(cos(dp) / cos(xtr)) * r

For one value (in my example, value 39), cos(dp) / cos(xtr) is bigger than 1 (but it is obviously 1, its difference with 1 is 2.220446e-16). Because cos is between -1 and 1, it produces a NaN.

A (quick and dirty) solution I coded in the function alongTrackDistance:

vec <- cos(dp) / cos(xtr)
if (any(abs(vec) > 1))
{
    l_sup <- length(vec[vec > 1]) != 0
    l_inf <- length(vec[vec < -1]) != 0

    if (l_sup & all.equal(vec[vec > 1], rep(1, length(vec[vec > 1]))))
        vec[vec > 1] <- 1
    if (l_inf & all.equal(vec[vec < -1], rep(-1, length(vec[vec < -1]))))
        vec[vec < -1] <- -1
}

dist <- bearing * acos(vec) * r

Therefore, depending how the computer coded numbers, this error might occur