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.9998, Inf, 179.9998, -89.1,
> Inf, -89.1), .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@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
pEpkey.asc
Description: application/pgp-keys
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
