search

r-spatial / sf

Simple Features for R
https://r-spatial.github.io/sf/
Other
1.35k stars 300 forks source link

`st_transform` docs are outdated in relation to projections 'not supported by GDAL' #1956

Closed Robinlovelace closed 2 years ago

Robinlovelace commented 2 years ago

The documentation of st_transform states that

Projecting simple feature geometries to projections not supported by GDAL may be done by st_transform_proj, part of package lwgeom.

As discussed with @Nowosad it is not clear that this is still true. We have tried many of the projections listed here: https://proj.org/operations/projections/index.html

All of the ones we've tried work with st_transform. A simple fix to the docs, if this is indeed the case (I'm not 100% sure), would be to remove that line in the docs. If there are examples where only st_transform_proj it may be worth mentioning them.

Context is updating Geocomputation with R following reviewer comments:

https://github.com/Robinlovelace/geocompr/commit/32e89551c34e92058124a463dd6ee1e273021cf4

This seems to be an improvement, reducing the number of functions people need to remember for reprojecting sf objects, if I understand correctly.

edzer commented 2 years ago

if this is indeed the case

Evidence that this is indeed the case would be helpful!

Robinlovelace commented 2 years ago

The only documented evidence I have of this so far is:

world_wintri = lwgeom::st_transform_proj(world, crs = "+proj=wintri")
world_wintri2 = st_transform(world, crs = "+proj=wintri")
waldo::compare(world_wintri$geom[1], world_wintri2$geom[1])

From the commit above. I tried a few other CRS projections that used to fail. All that is needed is a single example of when st_transform() fails but st_transform_proj works for the current wording in the docs to be correct. If anyone knows of such an example let us know, if not I would assume that it is indeed the case!

rsbivand commented 2 years ago

Please see ?rgdal::project for code to generate a list of known projections:

projs <- as.character(projInfo()$name)
res <- logical(length(projs))
names(res) <- projs
msgs <- character(length(projs))
names(msgs) <- projs
owarn <- options("warn")$warn
options(warn=2L)
for (i in seq(along=res)) {
  iprs <- paste("+proj=", projs[i], sep="")
  xy <- try(project(cbind(0, 0), iprs, legacy=TRUE, use_aoi=FALSE), silent=TRUE)
  if (inherits(xy, "try-error")) {
    res[i] <- NA
    msgs[i] <- paste("fwd:", strsplit(xy, "\n")[[1]][2])
  } else if(any(abs(xy) > 1e+08)) {
    res[i] <- NA
    msgs[i] <- paste("fwd: huge value")
  } else {
    out <- try(project(xy, iprs, inv=TRUE, legacy=TRUE, use_aoi=FALSE), silent=TRUE)
    if (inherits(out, "try-error")) {
      res[i] <- NA
      msgs[i] <- paste("inv:", strsplit(out, "\n")[[1]][2])
    } else {
      res[i] <- isTRUE(all.equal(cbind(0,0), out))
    }
  }
}
options(warn=owarn)
df <- data.frame(res=unname(res), msgs=unname(msgs), row.names=names(res))
# projection and inverse projection failures
# fwd: missing parameters
# inv: mostly inverse not defined
df[is.na(df$res),]
# inverse not equal to input
# (see http://lists.maptools.org/pipermail/proj/2011-November/006015.html)
df[!is.na(df$res) & !df$res,]
# inverse equal to input
row.names(df[!is.na(df$res) & df$res,])

which started to see which declared projections provide inverses. Output from current PROJ (many forward projections fail in rgdal (would fail in lwgeom because also using PROJ directly like rgdal) because extra information is required):

> df <- data.frame(res=unname(res), msgs=unname(msgs), row.names=names(res))
> # projection and inverse projection failures
> # fwd: missing parameters
> # inv: mostly inverse not defined
> df[is.na(df$res),]
             res
adams_hemi    NA
adams_ws1     NA
aea           NA
airy          NA
alsk          NA
apian         NA
august        NA
axisswap      NA
bacon         NA
bertin1953    NA
boggs         NA
bonne         NA
ccon          NA
chamb         NA
defmodel      NA
deformation   NA
denoy         NA
eqdc          NA
euler         NA
geoc          NA
geos          NA
gins8         NA
gn_sinu       NA
gs50          NA
guyou         NA
helmert       NA
hgridshift    NA
horner        NA
imw_p         NA
isea          NA
labrd         NA
larr          NA
lask          NA
lcc           NA
lcca          NA
lsat          NA
misrsom       NA
molobadekas   NA
molodensky    NA
murd1         NA
murd2         NA
murd3         NA
nicol         NA
nsper         NA
nzmg          NA
ob_tran       NA
oea           NA
omerc         NA
ortel         NA
pconic        NA
peirce_q      NA
pipeline      NA
rpoly         NA
sch           NA
tcc           NA
tinshift      NA
tissot        NA
topocentric   NA
tpeqd         NA
tpers         NA
unitconvert   NA
urm5          NA
urmfps        NA
utm           NA
vandg2        NA
vandg3        NA
vandg4        NA
vitk1         NA
vgridshift    NA
wag7          NA
xyzgridshift  NA
                                                                         msgs
adams_hemi                                    inv:   malformed pipeline found
adams_ws1                                     inv:   malformed pipeline found
aea          fwd:   target crs creation failed: Invalid value for an argument
airy                                          inv:   malformed pipeline found
alsk                                                          fwd: huge value
apian                                         inv:   malformed pipeline found
august                                        inv:   malformed pipeline found
axisswap      fwd:   target crs creation failed: Mutually exclusive arguments
bacon                                         inv:   malformed pipeline found
bertin1953                                    inv:   malformed pipeline found
boggs                                         inv:   malformed pipeline found
bonne        fwd:   target crs creation failed: Invalid value for an argument
ccon         fwd:   target crs creation failed: Invalid value for an argument
chamb        fwd:   target crs creation failed: Invalid value for an argument
defmodel                  fwd:   target crs creation failed: Missing argument
deformation               fwd:   target crs creation failed: Missing argument
denoy                                         inv:   malformed pipeline found
eqdc         fwd:   target crs creation failed: Invalid value for an argument
euler        fwd:   target crs creation failed: Invalid value for an argument
geoc             fwd:   source crs creation failed: Unknown error (code 4096)
geos         fwd:   target crs creation failed: Invalid value for an argument
gins8                                         inv:   malformed pipeline found
gn_sinu                   fwd:   target crs creation failed: Missing argument
gs50                                                          fwd: huge value
guyou                                         inv:   malformed pipeline found
helmert          fwd:   source crs creation failed: Unknown error (code 4096)
hgridshift                fwd:   target crs creation failed: Missing argument
horner                    fwd:   target crs creation failed: Missing argument
imw_p        fwd:   target crs creation failed: Invalid value for an argument
isea                                          inv:   malformed pipeline found
labrd        fwd:   target crs creation failed: Invalid value for an argument
larr                                          inv:   malformed pipeline found
lask                                          inv:   malformed pipeline found
lcc          fwd:   target crs creation failed: Invalid value for an argument
lcca         fwd:   target crs creation failed: Invalid value for an argument
lsat         fwd:   target crs creation failed: Invalid value for an argument
misrsom      fwd:   target crs creation failed: Invalid value for an argument
molobadekas               fwd:   target crs creation failed: Missing argument
molodensky                fwd:   target crs creation failed: Missing argument
murd1        fwd:   target crs creation failed: Invalid value for an argument
murd2        fwd:   target crs creation failed: Invalid value for an argument
murd3        fwd:   target crs creation failed: Invalid value for an argument
nicol                                         inv:   malformed pipeline found
nsper        fwd:   target crs creation failed: Invalid value for an argument
nzmg                                                          fwd: huge value
ob_tran                   fwd:   target crs creation failed: Missing argument
oea          fwd:   target crs creation failed: Invalid value for an argument
omerc        fwd:   target crs creation failed: Invalid value for an argument
ortel                                         inv:   malformed pipeline found
pconic       fwd:   target crs creation failed: Invalid value for an argument
peirce_q      inv:   (converted from warning) 1 projected point(s) not finite
pipeline        fwd:   target crs creation failed: Invalid PROJ string syntax
rpoly                                         inv:   malformed pipeline found
sch                       fwd:   target crs creation failed: Missing argument
tcc                                           inv:   malformed pipeline found
tinshift                  fwd:   target crs creation failed: Missing argument
tissot       fwd:   target crs creation failed: Invalid value for an argument
topocentric               fwd:   target crs creation failed: Missing argument
tpeqd        fwd:   target crs creation failed: Invalid value for an argument
tpers        fwd:   target crs creation failed: Invalid value for an argument
unitconvert      fwd:   source crs creation failed: Unknown error (code 4096)
urm5                      fwd:   target crs creation failed: Missing argument
urmfps                    fwd:   target crs creation failed: Missing argument
utm              fwd:   target crs creation failed: Unknown error (code 4096)
vandg2                                        inv:   malformed pipeline found
vandg3                                        inv:   malformed pipeline found
vandg4                                        inv:   malformed pipeline found
vitk1        fwd:   target crs creation failed: Invalid value for an argument
vgridshift                fwd:   target crs creation failed: Missing argument
wag7                                          inv:   malformed pipeline found
xyzgridshift              fwd:   target crs creation failed: Missing argument
> # inverse not equal to input
> # (see http://lists.maptools.org/pipermail/proj/2011-November/006015.html)
> df[!is.na(df$res) & !df$res,]
         res msgs
gs48   FALSE     
lee_os FALSE     
> # inverse equal to input
> row.names(df[!is.na(df$res) & df$res,])
  [1] "adams_ws2"  "aeqd"       "affine"     "aitoff"     "bipc"      
  [6] "calcofi"    "cart"       "cass"       "cc"         "cea"       
 [11] "collg"      "col_urban"  "comill"     "crast"      "eck1"      
 [16] "eck2"       "eck3"       "eck4"       "eck5"       "eck6"      
 [21] "eqearth"    "eqc"        "etmerc"     "fahey"      "fouc"      
 [26] "fouc_s"     "gall"       "geogoffset" "gnom"       "goode"     
 [31] "hammer"     "hatano"     "healpix"    "rhealpix"   "igh"       
 [36] "igh_o"      "kav5"       "kav7"       "krovak"     "laea"      
 [41] "lagrng"     "lonlat"     "latlon"     "leac"       "loxim"     
 [46] "mbt_s"      "mbt_fps"    "mbtfpp"     "mbtfpq"     "mbtfps"    
 [51] "merc"       "mil_os"     "mill"       "moll"       "natearth"  
 [56] "natearth2"  "nell"       "nell_h"     "noop"       "ocea"      
 [61] "ortho"      "patterson"  "poly"       "pop"        "push"      
 [66] "putp1"      "putp2"      "putp3"      "putp3p"     "putp4p"    
 [71] "putp5"      "putp5p"     "putp6"      "putp6p"     "qua_aut"   
 [76] "qsc"        "robin"      "rouss"      "s2"         "set"       
 [81] "sinu"       "somerc"     "stere"      "sterea"     "gstmerc"   
 [86] "tcea"       "times"      "tmerc"      "tobmerc"    "ups"       
 [91] "vandg"      "wag1"       "wag2"       "wag3"       "wag4"      
 [96] "wag5"       "wag6"       "webmerc"    "weren"      "wink1"     
[101] "wink2"      "wintri"
rsbivand commented 2 years ago

Roughly:

library(sf)
projs <- sf_proj_info(type = "proj")$name
res <- logical(length(projs))
names(res) <- projs
msgs <- character(length(projs))
names(msgs) <- projs
owarn <- options("warn")$warn
options(warn=2L)
source <- st_sfc(list(st_point(matrix(c(0,0), nrow=1))), crs="OGC:CRS84")
for (i in seq(along=res)) {
  iprs <- paste("+proj=", projs[i], sep="")
  xy <- try(st_transform(source, iprs))
  if (inherits(xy, "try-error")) {
    res[i] <- NA
    msgs[i] <- paste("fwd:", strsplit(xy, "\n")[[1]][2])
}
}
library(lwgeom)
res1 <- logical(length(projs))
names(res1) <- projs
msgs1 <- character(length(projs))
names(msgs1) <- projs
for (i in seq(along=res1)) {
   iprs <- paste("+proj=", projs[i], sep="")
  xy <- try(st_transform_proj(source, iprs))
  if (inherits(xy, "try-error")) {
    res1[i] <- NA
    msgs1[i] <- paste("fwd:", strsplit(xy, "\n")[[1]][2])
}}
options(warn=owarn)

Output objects in RDA file in zipfile: proj_comparison.zip

rsbivand commented 2 years ago

It looks as though PROJ-direct and PROJ-through-GDAL agree:

> all.equal(nchar(msgs) == 0L, nchar(msgs1) == 0L)
[1] TRUE
> msgs["vitk1"]
                                                                                                                                                             vitk1 
"fwd:   (converted from warning) GDAL Error 1: PROJ: proj_create: Error 1027 (Invalid value for an argument): vitk1: Missing parameter: lat_1 should be specified" 
> msgs1["vitk1"]
                                                           vitk1 
"fwd:   st_lwgeom_transform: one of the proj strings is invalid"

with messages (here all errors) only varying in that GDAL messages are more explanatory. So removing the reference to lwgeom::st_transform_proj() seems justified.

Robinlovelace commented 2 years ago

So removing the reference to lwgeom::st_transform_proj() seems justified.

Done in https://github.com/r-spatial/sf/pull/1957