PROJ6 rgdal::project/sp::spTransform - issue of points that cannot be inverted

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: PROJ6 rgdal::project/sp::spTransform - issue of points that cannot be inverted

edzer
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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

pEpkey.asc (3K) Download Attachment