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

[Edzer Pebesma]

I don't know why Roger wants to deprecate rgdal::project, but for
various reasons I implemented a similar function in sf:

library(sf)
# Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
ll <- cbind(rep(180, 3), c(-89, -90, -89))
xy = sf_project("+proj=longlat", "+proj=moll", ll)
sf_project("+proj=moll", "+proj=longlat", xy)
#      [,1] [,2]
# [1,]  180  -89
# [2,]  Inf  Inf
# [3,]  180  -89
# Warning message:
# In CPL_proj_direct(as.character(c(from[1], to[1])), as.matrix(pts)) :
#   one or more projected point(s) not finite

As it stands now, I want to keep this function in sf, and would be happy
to add an argument to suppress the warning msg.



On 12/6/19 3:25 PM, Daniel Kelley wrote:
> 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
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

-------------- next part --------------
A non-text attachment was scrubbed...
Name: pEpkey.asc
Type: application/pgp-keys
Size: 3110 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20191206/dcf0e6bd/attachment.bin>