heads-up re map projection

[dankelley]

Below is a msg posted to r-sig-geo. It may affect oce eventually, if/when rgdal changes. (Note that oce does not use the datum= token directly in its calls to rgdal.)

I'll pin this issue, so it doesn't autoclose over time.

Message: 2
Date: Wed, 11 Sep 2019 14:34:59 +0200
From: Roger Bivand Roger.Bivand@nhh.no
To: r-sig-geo@r-project.org
Subject: [R-sig-Geo] PROJ 6 may impact workflows and packages using sf
or rgdal
Message-ID: alpine.LFD.2.21.1909111419590.17735@reclus.nhh.no
Content-Type: text/plain; charset="us-ascii"; Format="flowed"

Changes introduced in PROJ 6, whether using the new proj.h API (sf) or the
old proj_api.h API (rgdal, lwgeom) may impact R workflows and packages
using sf or rgdal for transforming coordinate reference systems.

Because +datum= is effectively defunct, and will be gone in PROJ 7 (March
2020), any data including CRS set correctly for previous PROJ will risk no
longer be transformed correctly. Work is ongoing to establish how users
can adapt to this, but it will be very risky to assume that
transformations that worked in PROJ 4 will necessarily work as well or
better in PROJ >= 6 without manual intervention.

Discussion is welcomed, in particular on:
r-spatial/discuss#28 based on
r-spatial/sf#1146; further analysis described in
the workshop document on https://rsbivand.github.io/geostat19_talk/.

Please examine your workflows and packages, check whether your work is
affected, and contribute to discussions through the issues on github. This
list may be used, but it'll be best to keep things connected.

Roger

--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: Roger.Bivand@nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

[dankelley]

http://rgdal.r-forge.r-project.org/articles/PROJ6_GDAL3.html is a good source for information on the upcoming plans. This makes me wonder whether we ought to be using CRS() inside oce, to check map-projection strings for correctness, so we can handle errors in a gentler way at the oce level.

[dankelley]

Here is a snippet of a posting on r-package-devel-request by R. Bivand:

The development version of rgdal on R-Forge is now at rev 894, and is now
ready for trying out with PROJ6/GDAL3 workflows, and workflows that may
migrate within 6 months to modern CRS representations. The motivating RFC
is also updated to cover coordinate operations, the use of prepared
(pre-searched) coordinate operations, and should be read carefully by
anyone using rgdal::spTransform(). Note further that rgdal::project() will
not be adapted for PROJ6, and is effectively deprecated.

Since oce uses rgdal::project(), this applies to us. I do not fully remember the decision to use rgdal::project(), although I do recall trying some different things. maybe we'll ned to switch to rgdal::spTransform(). I will look into this.

[dankelley]

I did a test, as below, and the procedure does not look arduous. I have sp version 1.3-2 and rgdal version 1.4.7, with R 4.0 i.e. R devel (2019/11/07, r77386) downloaded today.

I'm a bit confused, because I get an error when I try rgdal::spTransform(), but anyway it works with sp::spTransform() so I'm trying that.

requireNamespace("rgdal")
requireNamespace("sp")

data(coastlineWorld, package="oce")
lon <- coastlineWorld[["longitude"]]
lat <- coastlineWorld[["latitude"]]
## cannot have NAs for
#sp::SpatialPoints(cbind(lon, lat)) # fails 'NA values in coordinates'
lonlat <- cbind(lon, lat)
isNA <- is.na(lon)
lonlat[isNA,] <- c(0, 0)
lonlatSpatial <- sp::SpatialPoints(lonlat, sp::CRS("+proj=longlat"))

## NB. rgdal does not export spTransform() but sp does.
lonlatProjected <- sp::spTransform(lonlatSpatial, sp::CRS("+proj=moll"))

xy <- sp::coordinates(lonlatProjected)
xy[isNA] <- c(NA, NA)

if (!interactive()) png("oceproj01.png")
par(mar=c(3,3,1,1))
plot(xy, type="l", asp=1, xlab="easting [m]", ylab="northing [m]")
grid()
polygon(xy, col='lightgray')
if (!interactive()) dev.off()

[dankelley]

NB. git grep -i :project in the R directory yields as follows, so a transition from rgdal::project() to sp::spTransform() may not especially time-consuming.

argo.R:1531:#' a projection in the notation used by the [rgdal::project] function
coastline.R:1421:#' call to the [rgdal::project()] function in the \CRANpkg{rgdal} package.
map.R:13:#' Wrapper to rgdal::project()
map.R:16:#' changes to the [rgdal::project()] function in the \CRANpkg{rgdal}
map.R:23:#'    to [rgdal::project()], both functions provided by the
map.R:25:#' 2. 2016 Apr: rgdal::project started returning named quantities
map.R:33:#' @param xy,proj,inv,use_ob_tran,legacy  As for the [rgdal::project()] function in the
map.R:57:                           XY <- unname(rgdal::project(xy, proj=proj, inv=inv))
map.R:67:                               XY <- unname(rgdal::project(xy, proj=proj, inv=inv, legacy=legacy, allowNAs_if_not_legacy=TRUE))
map.R:73:                               XY <- unname(rgdal::project(xy=xy, proj=proj, inv=inv, legacy=legacy))
map.R:2628:#' `oce` uses [rgdal::project()] in the \CRANpkg{rgdal}

[dankelley]

Hm, transferring from rgdal::project() to sp::spTransform() is not going to be easy, because the former returns good results where possible, whereas the latter returns just an error and no results, e.g.

lonlat <- sp::SpatialPoints(cbind(rep(180,3), c(-89,-89.2,-89)), sp::CRS("+proj=longlat"))
xy <- sp::spTransform(lonlat, sp::CRS("+proj=moll"))
rgdal::project(sp::coordinates(xy), proj="+proj=moll", inv=TRUE)

gives

     coords.x1 coords.x2
[1,]       180       -89
[2,]       Inf       Inf
[3,]       180       -89

whereas

try(sp::spTransform(xy, sp::CRS("+proj=longlat")), silent=TRUE)

gives

[1] "Error in sp::spTransform(xy, sp::CRS(\"+proj=longlat\")) : \n  failure in points 2\n"
attr(,"class")
[1] "try-error"
attr(,"condition")
<simpleError in sp::spTransform(xy, sp::CRS("+proj=longlat")): failure in points 2>

i.e. we get no results if any point is problematic.

I've done some experimentation to find out how close we can get to the "edge of the earth", but oce will become a wild tangle of data-editing code if we have to try to avoid such spots for general projections, where the edge of the earth may be at any longitude and latitude (e.g. consider the effect of +lon_0 and so forth).

[dankelley]

In looking at the rgal-1.4-8 code (e.g. R/project.R line 179) I see that there are no args to spTransform.SpatialPoints() that are used to avoid this error; the test at R/project.R line 181 will call stop() even for a single non-finite result.

Although I think the best plan is to use a package's exported functions, we could perhaps solve th problem with direct calls, as tested in the following:

library("rgdal")
library("sp")
ll <- cbind(rep(180, 3), c(-89, -90, -89))
CRSll <- sp::CRS("+proj=longlat")
lonlat <- sp::SpatialPoints(ll, CRSll)
CRSxy <- sp::CRS("+proj=moll")

## test direct call to transform
xyDEK <- .Call("transform", proj4string(lonlat), slot(CRSxy, "projargs"),
               n=dim(ll)[1], as.double(ll[,1]), as.double(ll[,2]), NULL,
               PACKAGE="rgdal")
xy <- sp::spTransform(lonlat, sp::CRS("+proj=moll"))
stopifnot(all.equal(xyDEK[[1]], coordinates(xy)[,1]))
stopifnot(all.equal(xyDEK[[2]], coordinates(xy)[,2]))

[dankelley]

Below shows (for this example) how to do both forward and inverse. The is ugly, with bad variable names, but I see little reason to clean it up before making a decision on whether to access this rgdal function with a .Call (and I need to check CRAN policies on even doing that).

> library("rgdal")
> library("sp")
> ll <- cbind(rep(180, 3), c(-89, -90, -89))
> CRSll <- sp::CRS("+proj=longlat")
> lonlat <- sp::SpatialPoints(ll, CRSll)
> CRSxy <- sp::CRS("+proj=moll")
> 
> ## Check 'forward' ddirect call
> xyDEK <- .Call("transform", proj4string(lonlat), slot(CRSxy, "projargs"),
+                n=dim(ll)[1], as.double(ll[,1]), as.double(ll[,2]), NULL,
+                PACKAGE="rgdal")
> xDEK <- xyDEK[[1]]
> yDEK <- xyDEK[[2]]
> xy <- coordinates(sp::spTransform(lonlat, CRSxy))
> stopifnot(all.equal(xyDEK[[1]], xy[,1]))
> stopifnot(all.equal(xyDEK[[2]], xy[,2]))
> 
> ## Check 'inverse' direct call
> owarn <- options("warn")$warn
> options(warn=-1)
> llDEK <- .Call("transform", slot(CRSxy, "projargs"), slot(CRSll, "projargs"),
+                n=length(xDEK), as.double(xDEK), as.double(yDEK), NULL,
+                PACKAGE="rgdal")
> LL <- rgdal::project(sp::coordinates(xy), proj="+proj=moll", inv=TRUE)
> options(warn=owarn)
> stopifnot(all.equal(llDEK[[1]], LL[,1]))
> stopifnot(all.equal(llDEK[[2]], LL[,2]))
> 
> data.frame(lon=ll[,1], lat=ll[,2],
+            xDEK=xyDEK[[1]], yDEK=xyDEK[[2]],
+            xRGDAL=xy[,1], yRGDAL=xy[,2],
+            lonDEK=llDEK[[1]], latDEK=llDEK[[2]])
  lon lat         xDEK     yDEK       xRGDAL   yRGDAL lonDEK latDEK
1 180 -89 1.281330e+06 -8997267 1.281330e+06 -8997267    180    -89
2 180 -90 1.104637e-09 -9020048 1.104637e-09 -9020048    Inf    Inf
3 180 -89 1.281330e+06 -8997267 1.281330e+06 -8997267    180    -89

[dankelley]

I see that using .Call() to a different package is a violation of CRAN rules (https://cran.r-project.org/doc/manuals/R-exts.html#Writing-portable-packages)

It is not portable to call compiled code in R or other packages via .Internal, .C, .Fortran, .Call or .External, since such interfaces are subject to change without notice and will probably result in your code terminating the R process.

[dankelley]

A note that Edzer Pebesma posted on r-sig-geo@r-project.org (Fri, 6 Dec 2019 15:53:37 +0100) suggests that sf::sf_project() may be a better approach, since it is not slated for withdrawal, as rgdal::project() is, and since it does not have this sp:spTransform() difficulty of losing all output if even a single point cannot be projected.

I've just read the following, and it makes me think that sf is a good package for oce users to be familiar with, as well.

@article{RJ-2018-009,
  title = {Simple {{Features}} for {{R}}: {{Standardized Support}} for {{Spatial Vector Data}}},
  volume = {10},
  url = {https://doi.org/10.32614/RJ-2018-009},
  doi = {10.32614/RJ-2018-009},
  number = {1},
  journaltitle = {The R Journal},
  date = {2018},
  pages = {439-446},
  author = {Pebesma, Edzer}
}

[dankelley]

Having read a few more articles about sf, I am becoming inclined to propose that oce switch to it, for things now done with rgdal etc. I suspect (but have not checked) that we can use sf for all the following (a - means not changed yet; a + means changed)

git grep -En '^[ a-zA-Z].(raster::)|(sp::)|(raster::)' R
+R/argo.R:534:                      keep <- 1==sp::point.in.polygon(lon, lat, lonp, latp)
-R/coastline.R:222:              box <- as(raster::extent(W, E, S, N), "SpatialPolygons")
-R/coastline.R:239:                      A <- sp::Polygon(cbind(lon, lat))
-R/coastline.R:240:                      B <- sp::Polygons(list(A), "A")
-R/coastline.R:241:                      C <- sp::SpatialPolygons(list(B))
-R/coastline.R:242:                      i <- raster::intersect(box, C)
-R/map.R:786:                            erase <- 1==sp::point.in.polygon(xc, yc,
-R/section.R:678:                      keep <- 1==sp::point.in.polygon(lon, lat, lonp, latp)
-R/section.R:2936:                                                    input_data=sp::SpatialPointsDataFrame(xy, data.frame(F)),
-R/section.R:2937:                                                    new_data=sp::SpatialPoints(expand.grid(xg/xr, yg/yr)))

[dankelley]

I have started a new branch called sf in which I'll test the things listed in the previous comment. So far, I've just done one test, which is for subset,argo-method. This is in commit ebfab31 of the sf branch. To test it, do as follows:

options("oce:test_sf"=1)

and then run the 4th example for ?'subset,argo-method':

par(mfrow=c(1, 2))
plot(argo, which="map")
## Can get a boundary with e.g. locator(4)
bdy <- list(x=c(-65, -40, -40, -65), y=c(65, 65, 45, 45))
argoSubset <- subset(argo, within=bdy)
plot(argoSubset, which="map")

which prints a few lines before and after the test. (The test works, which makes me think that this will work in any case ... either I figured out how to do an intersection in sf or I did not, and this suggests that I did.) I will likely do the other point-in-polygon cases sometime during the upcoming week.

[dankelley]

I am closing this since it is superseded by #1629.