should oce use rgdal::spTransform()?

[dankelley]

I sort of think this is the recommended method, as opposed to rgdal::project(). However, I am having a hard time decoding the docs to find out if there are parameters that will make rgdal::spTransform() return NA or Inf for points that cannot be projected. I wonder, does @paleolimbot happen to know? (I can also ask on r-sig-geo@r-project.org but thought Dewey might know off the top of his head.)

# Demo that spTransform() returns nothing, if a point cannot be projected.

library(rgdal)
#> Loading required package: sp
#> rgdal: version: 1.5-23, (SVN revision 1121)
#> Geospatial Data Abstraction Library extensions to R successfully loaded
#> Loaded GDAL runtime: GDAL 3.1.4, released 2020/10/20
#> Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.1/Resources/library/rgdal/gdal
#> GDAL binary built with GEOS: TRUE 
#> Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
#> Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.1/Resources/library/rgdal/proj
#> Linking to sp version:1.4-5
#> To mute warnings of possible GDAL/OSR exportToProj4() degradation,
#> use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.

# 1. Works
n <- 5
lon <- rep(30, n)
lat <- seq(-60, 60, length.out=n)
ll <- data.frame(x=lon, y=lat)
coordinates(ll) <- c("x", "y")
proj4string(ll) <- CRS("+proj=longlat")
spTransform(ll, CRS("+proj=ortho +lat_0=-20"))
#> SpatialPoints:
#>            x        y
#> [1,] 1594534 -4245917
#> [2,] 2761814 -1360656
#> [3,] 3189068  1889192
#> [4,] 2761814  4632833
#> [5,] 1594534  6135109
#> Coordinate Reference System (CRS) arguments: +proj=ortho +lat_0=-20
#> +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

# 2. Fails on 30E, 70N. Q: is there an arg that makes spTransform() return
# Inf/NA as needed?
n <- 5
lon <- rep(30, n)
lat <- seq(-70, 70, length.out=n)
ll <- data.frame(x=lon, y=lat)
coordinates(ll) <- c("x", "y")
proj4string(ll) <- CRS("+proj=longlat")
spTransform(ll, CRS("+proj=ortho +lat_0=-20"))
#> Warning in spTransform(ll, CRS("+proj=ortho +lat_0=-20")): 1 projected point(s)
#> not finite
#> non finite transformation detected:
#>   x   y         
#>  30  70 Inf Inf
#> Error in spTransform(ll, CRS("+proj=ortho +lat_0=-20")): failure in points 5

^{Created on 2021-07-16 by the reprex package (v2.0.0)}

[paleolimbot]

Is sf::sf_project() not going to work here? I think you probably want to remove all references to rgdal since it's getting retired next year. Pretty sure it can handle NA (but if not you can do a subset assign).

sf::sf_project("WGS84", "EPSG:3857", matrix(c(0, 0, -64, NA), nrow = 2, byrow = TRUE))
#>          [,1] [,2]
#> [1,]        0    0
#> [2,] -7124447  NaN

pts <- matrix(c(0, 0, -64, NA), nrow = 2, byrow = TRUE)
pts_finite <- is.finite(pts[, 1]) & is.finite(pts[, 2])
pts[pts_finite, ] <- sf::sf_project("WGS84", "EPSG:3857", pts[pts_finite, , drop = FALSE])
pts[!pts_finite, ] <- NA_real_
pts
#>      [,1] [,2]
#> [1,]    0    0
#> [2,]   NA   NA

^{Created on 2021-07-16 by the reprex package (v2.0.0)}

[dankelley]

Oh, I didn't realize rgdal was being retired. Wow. Anyway, the reason I'm investigating this is that it seems to me (and I'm likely just wrong) that sf::sf_project() might have some problems; see #1848 (comment) and r-spatial/sf#1729.

[dankelley]

I am going to close this issue now, because I am getting good results in sf, and I surely don't want to switch a lot of work to rgdal if it's going to be retired soon. The answer to the question in the issue title is, therefore, "no".