[R-sig-Geo] PROJ6 rgdal::project/sp::spTransform - issue of points that cannot be inverted

[Daniel Kelley]

>From prior discussions on this thread and elsewhere, I am given to understand that rgdal::project() will not be carried over into PROJ6, and that we are recommendd to use sp::spTransform() instead.

Accordingly, I started work on switching the 'oce' package from rgdal::project() to sp::spTransform().  However, I have a problem with points that cannot be inverted.  With rgdal::project(), we get warning messages if the input contains such points (and the output is Inf for the associated projected points).  However, it seems that with sp::spTransform(), we get an error message and no results, if any of the data contain points that cannot be inverted.  I want to be able to transform a lot of points at the same time (owing to data-set size and speed considerations), so handling the data point-by-point is not an option.

My question is simple: is there a way I can make sp::spTransform() return an equivalent to that of rgdal::project(), with finite data where possible and Inf (or similar) where an inverse cannot be done?

Possibly I am just something.  For anyone who has the patience to look, I have some thoughts at https://github.com/dankelley/oce/issues/1599#issuecomment-562578641 but the gist is as follows: note that 'PROJ' consists of three points, one of which is (of course) Inf, but that 'TRAN' consists only of an error message. Maybe there is an argument to sp::spTransform() that I'm unaware of, that will cause it to return data when it can, and Inf when it cannot?

```
R Under development (unstable) (2019-11-07 r77386) -- "Unsuffered Consequences"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin15.6.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("rgdal")
> library("sp")
> ll <- cbind(rep(180, 3), c(-89, -90, -89))
> lonlat <- sp::SpatialPoints(ll, sp::CRS("+proj=longlat"))
> xy <- sp::spTransform(lonlat, sp::CRS("+proj=moll"))
> # rgdal::project returns a mixture of results and Inf values
> dput(rgdal::project(sp::coordinates(xy), proj="+proj=moll", inv=TRUE))
structure(c(179.999999999998, Inf, 179.999999999998, -89.0000000000001, 
Inf, -89.0000000000001), .Dim = 3:2, .Dimnames = list(NULL, c("coords.x1", 
"coords.x2")))
> ## but sp::spTransform returns no data at all
> dput(try(sp::spTransform(xy, sp::CRS("+proj=longlat")), silent=TRUE))
non finite transformation detected:
    coords.x1     coords.x2                             
 1.104637e-09 -9.020048e+06           Inf           Inf 
structure("Error in sp::spTransform(xy, sp::CRS(\"+proj=longlat\")) : \n  failure in points 2\n", class = "try-error", condition = structure(list(
    message = "failure in points 2", call = sp::spTransform(xy, 
        sp::CRS("+proj=longlat"))), class = c("simpleError", 
"error", "condition")))
> 
```



Dan Kelley
Department of Oceanography
Dalhousie University
Halifax, NS, Canada