Sf: PROJ 6.3.0 dropping ellps in P4 string

[rsbivand]

More https://github.com/r-spatial/sf/issues/1187 fun.

See https://lists.osgeo.org/pipermail/proj/2019-December/009166.html. Even seems to think that the OGRSpatialReference object may be subject to short-term volatility. rgdal functions P6_SRID_show() and ogrP4S() affected. All works fine in PROJ 6.2.1, for EPSG:27700 +ellps= and +units= omitted in exportToProj4(). I tried to check sf but don't know which branch I'm checking, there is an error in R check. Hope PRØJ gives us a break in 2020.

[rsbivand]

Problem reduced to one of finding out which environment/settings knitr is using in building the PROJ6_GDAL3.Rmd vignette. Failing chunks succeed when run manually in RStudio and in the R console (inside and outside RStudio), but fail with +ellps= and units= when building the vignette for EPSG:27700. Starting a test vignette to see if this can be narrowed down.

[rsbivand]

Test notebook, fails at critical points (renamed *.txt, move back to *.Rmd):

PROJ6_GDAL3_test.txt

[rsbivand]

Snippet failing:

library(sp)
(crs <- CRS("+init=epsg:27700"))
library(rgdal)
(discarded_datum <- showSRID("EPSG:27700", "PROJ"))
(crs <- CRS("+init=epsg:28992"))
(discarded_datum <- showSRID("EPSG:28992", "PROJ"))

> library(sp)
> (crs <- CRS("+init=epsg:27700"))
CRS arguments:
 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
+y_0=-100000 +ellps=airy +units=m +no_defs 
Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
  Discarded datum OSGB_1936 in CRS definition
> library(rgdal)
rgdal: version: 1.5-2, (SVN revision 907)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 3.1.0dev-84d972714a, released 2019/12/30
 Path to GDAL shared files: /usr/local/share/gdal
 GDAL binary built with GEOS: TRUE 
 Loaded PROJ.4 runtime: Rel. 6.3.0, January 1st, 2020, [PJ_VERSION: 630]
 Path to PROJ.4 shared files: /usr/local/share/proj
 Linking to sp version: 1.3-4 
> (discarded_datum <- showSRID("EPSG:27700", "PROJ"))
[1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +no_defs"
Warning message:
In showSRID("EPSG:27700", "PROJ") :
  Discarded datum OSGB_1936 in CRS definition
> (crs <- CRS("+init=epsg:28992"))
CRS arguments:
 +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889
+k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs 
Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
  Discarded datum Amersfoort in CRS definition
> (discarded_datum <- showSRID("EPSG:28992", "PROJ"))
[1] "+proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +no_defs"
Warning message:
In showSRID("EPSG:28992", "PROJ") :
  Discarded datum Amersfoort in CRS definition

[rsbivand]

Perhaps resolved with rev 910, declaring OGRSpatialReference object without setting to (OGRSpatialReference) NULL. However, not sure - the objects seem to be sticky.

[rsbivand]

Not resolved, CMD check passing but +ellps= and +units= being dropped on exportToProj4() when OGRSpatialReference created from EPSG. No ideas. Rev 911 succeeds in building vignettes but example above fails on command line.

[rsbivand]

Now failing vignette builds for scot_BNG <- readOGR(system.file("vectors", package="rgdal"), "scot_BNG") and in examples, but succeeds in console, at rev. 912 and PROJ 6.3.0/GDAL 3.1 master. Same revision passes check for PROJ 6.2.1/GDAL 3.0.2.

[rsbivand]

Rev 914 checks for PROJ < 6.3.0 and branches in vignette line 847 to avoid trying to rebuild the CRS from the broken PROJ string; in spTransform the grid example for scot_BNG does the same on line 144.

[rsbivand]

Rev 915 reverts the branching in 914 on >= 6.3.0, by intervening earlier. If the PROJ string is missing +ellps= where it is present in the input object read from vector data file, the PROJ string is rebuilt from the WKT2 string using showSRID(wkt2, "PROJ") in rgdal::OGRSpatialRef(). Passing on 630/31master, 621/302, https://r-forge.r-project.org/R/?group_id=884 gives status for 520/240.

[edzer]

I'll build a docker image with these GDAL/PROJ versions, and will check with sf.

[rsbivand]

Ok, thanks. I think it is about not handling or freeing OGRSpatialReference
right, and some change in proj makes the datum node vanish.

Roger Bivand
Falsensvei 32
5063 Bergen

ons. 1. jan. 2020, 20.14 skrev Edzer Pebesma [e-mail]:

I'll build a docker image with these GDAL/PROJ versions, and will check
with sf.

—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub
https://github.com/r-spatial/sf/issues/1231?email_source=notifications&email_token=ACNZ3BGEZOF4Q3GYZELB463Q3TTQLA5CNFSM4KBWKGWKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEH5K2GA#issuecomment-570076440,
or unsubscribe
https://github.com/notifications/unsubscribe-auth/ACNZ3BC73GGPFGOTOXHX2R3Q3TTQLANCNFSM4KBWKGWA
.

[edzer]

First results: with

library(sf)
# Linking to GEOS 3.8.0, GDAL 3.1.0dev-bb67ec3, PROJ 6.3.0

the following breaks now:

library(sf)
p1 = st_point(c(7,52))
p2 = st_point(c(-30,20))
sfc = st_sfc(p1, p2, crs = 4326)
st_transform(sfc, "+init=epsg:3857")

with

Error in CPL_transform(x, crs$proj4string) : 
  OGRCreateCoordinateTransformation() returned NULL: PROJ available?
In addition: Warning messages:
1: In CPL_crs_from_proj4string(x) :
  GDAL Message 1: +init=epsg:XXXX syntax is deprecated. It might return a CRS with a non-EPSG compliant axis order.
2: In CPL_transform(x, crs$proj4string) :
  GDAL Error 1: PROJ: proj_create_operations: At least one of the operation lacks a source and/or target CRS
3: In CPL_transform(x, crs$proj4string) :
  GDAL Error 6: Cannot find coordinate operations from `GEOGCRS["unknown",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]],CS[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude",north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]]' to `PROJCRS["WGS 84 / Pseudo-Mercator",BASEGEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["unnamed",METHOD["Popular Visualisation Pseudo Mercator",ID["EPSG",1024]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False eastin [... truncated]

The following does NOT break:

library(sf)
p1 = st_point(c(7,52))
p2 = st_point(c(-30,20))
sfc = st_sfc(p1, p2, crs = 4326)
try(st_transform(sfc, "+init=epsg:3857"))
st_transform(sfc, 3857)

but in the following, the LAST command now breaks too:

library(sf)
p1 = st_point(c(7,52))
p2 = st_point(c(-30,20))
sfc = st_sfc(p1, p2, crs = 4326)
try(st_transform(sfc, "+init=epsg:3857"))
st_transform(sfc, 3857)

valgrind gives strange memory errors, the whole thing smells fishy.

[rsbivand]

Thanks, now whack a fish after whack a mole? Is it replicable in C++? To
get Even to look, that would help, but my attempt in the proj list thread
wasn't successful.

Roger Bivand
Falsensvei 32
5063 Bergen

ons. 1. jan. 2020, 22.14 skrev Edzer Pebesma [e-mail]:

First results: with

library(sf)# Linking to GEOS 3.8.0, GDAL 3.1.0dev-bb67ec3, PROJ 6.3.0

the following breaks now:

library(sf)p1 = st_point(c(7,52))p2 = st_point(c(-30,20))sfc = st_sfc(p1, p2, crs = 4326)
st_transform(sfc, "+init=epsg:3857")

with

Error in CPL_transform(x, crs$proj4string) :
OGRCreateCoordinateTransformation() returned NULL: PROJ available?
In addition: Warning messages:
1: In CPL_crs_from_proj4string(x) :
GDAL Message 1: +init=epsg:XXXX syntax is deprecated. It might return a CRS with a non-EPSG compliant axis order.
2: In CPL_transform(x, crs$proj4string) :
GDAL Error 1: PROJ: proj_create_operations: At least one of the operation lacks a source and/or target CRS
3: In CPL_transform(x, crs$proj4string) :
GDAL Error 6: Cannot find coordinate operations from GEOGCRS["unknown",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]],CS[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude",north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]]' toPROJCRS["WGS 84 / Pseudo-Mercator",BASEGEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["unnamed",METHOD["Popular Visualisation Pseudo Mercator",ID["EPSG",1024]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False eastin [... truncated]

The following does NOT break:

library(sf)p1 = st_point(c(7,52))p2 = st_point(c(-30,20))sfc = st_sfc(p1, p2, crs = 4326)
try(st_transform(sfc, "+init=epsg:3857"))
st_transform(sfc, 3857)

but in the following, the LAST command now breaks too:

library(sf)p1 = st_point(c(7,52))p2 = st_point(c(-30,20))sfc = st_sfc(p1, p2, crs = 4326)
try(st_transform(sfc, "+init=epsg:3857"))
st_transform(sfc, 3857)

valgrind gives strange memory errors, the whole thing smells fishy.

—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub
https://github.com/r-spatial/sf/issues/1231?email_source=notifications&email_token=ACNZ3BFJSFLRC4GIYXKJA4TQ3UBULA5CNFSM4KBWKGWKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEH5MV7A#issuecomment-570084092,
or unsubscribe
https://github.com/notifications/unsubscribe-auth/ACNZ3BFNBLR73BG2J3BWIU3Q3UBULANCNFSM4KBWKGWA
.

[edzer]

No: I see

root@f2b1ca036605:/# echo "151.2093 -33" | cs2cs "+proj=longlat +datum=WGS84 +no_defs" "+init=epsg:3857"
16832542.28 -3895303.96 0.00
root@f2b1ca036605:/# echo "151.2093 -33" | cs2cs "+proj=longlat +datum=WGS84 +no_defs" "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"
16832542.28 -3895303.96 0.00

[rsbivand]

Right, staying just on the PROJ side seems consistent, I think. But when
GDAL gets involved, I think odder things start seeming to appear that
weren't there with 6.2.1.

Roger Bivand
Falsensvei 32
5063 Bergen

ons. 1. jan. 2020, 22.39 skrev Edzer Pebesma [e-mail]:

No: I see

root@f2b1ca036605:/# echo "151.2093 -33" | cs2cs "+proj=longlat +datum=WGS84 +no_defs" "+init=epsg:3857"
16832542.28 -3895303.96 0.00
root@f2b1ca036605:/# echo "151.2093 -33" | cs2cs "+proj=longlat +datum=WGS84 +no_defs" "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"
16832542.28 -3895303.96 0.00

—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub
https://github.com/r-spatial/sf/issues/1231?email_source=notifications&email_token=ACNZ3BELGCZVSO6CA2HYMOTQ3UEP3A5CNFSM4KBWKGWKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEH5NBVI#issuecomment-570085589,
or unsubscribe
https://github.com/notifications/unsubscribe-auth/ACNZ3BEEP5CFSNVJ5NU6W4TQ3UEP3ANCNFSM4KBWKGWA
.

[edzer]

On your script (posted above) I see:

root@85645a0ba6f0:/# R

R version 3.6.2 (2019-12-12) -- "Dark and Stormy Night"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (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.

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(sp)
> (crs <- CRS("+init=epsg:27700"))
CRS arguments:
 +init=epsg:27700 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717
+x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs 
> library(rgdal)
rgdal: version: 1.5-3, (SVN revision (unknown))
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 3.1.0dev-bb67ec3, released 2019/12/31
 Path to GDAL shared files: /usr/local/share/gdal
 GDAL binary built with GEOS: FALSE 
 Loaded PROJ.4 runtime: Rel. 6.3.0, January 1st, 2020, [PJ_VERSION: 630]
 Path to PROJ.4 shared files: /usr/local/share/proj
 Linking to sp version: 1.3-2 
> (discarded_datum <- showSRID("EPSG:27700", "PROJ"))
[1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs"
Warning message:
In showSRID("EPSG:27700", "PROJ") :
  Discarded datum OSGB_1936 in CRS definition
> (x <- list_coordOps(paste0(discarded_datum, " +type=crs"), "EPSG:4326"))
Candidate coordinate operations found:  1 
Strict containment:  FALSE 
Visualization order:  TRUE 
Source: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
        +y_0=-100000 +ellps=airy +units=m +no_defs +type=crs
Target: EPSG:4326 
Best instantiable operation has only ballpark accuracy 
Description: Inverse of unknown + Ballpark geographic offset from unknown to
             WGS 84 + axis order change (2D)
Definition:  +proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2
             +k=0.9996012717 +x_0=400000 +y_0=-100000
             +ellps=airy +step +proj=unitconvert +xy_in=rad
             +xy_out=deg
> scot_BNG <- readOGR(system.file("vectors", package="rgdal"), "scot_BNG")
OGR data source with driver: ESRI Shapefile 
Source: "/usr/local/lib/R/site-library/rgdal/vectors", layer: "scot_BNG"
with 56 features
It has 13 fields
Warning messages:
1: In OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS = dumpSRS,  :
  Discarded ellps Airy 1830 in CRS definition
2: In OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS = dumpSRS,  :
  Discarded datum OSGB_1936 in CRS definition: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +no_defs
> (load_status <- get_transform_wkt_comment())
[1] TRUE

(now building a container with GDAL from git & PROJ 6.2.0)

[rsbivand]

Done changing all rgdal OGRSpatialReference declarations to pointers, committed at rev 918. Don't see any changes so far but need to check properly.

exportToProj4() starting line 10007 gdal/ogr/ogrspatialreference.cpp calls proj_as_proj_string() starting line 1414 proj-6.3.0/src/iso19111/c_api.cpp twice (added see https://lists.osgeo.org/pipermail/gdal-dev/2019-December/051344.html) but not present in 3.0.2, so likely not the problem. The commits to the file containing proj_as_proj_string() are https://github.com/OSGeo/PROJ/commits/master/src/iso19111/c_api.cpp, but it may very well be another change somewhere else in a function called by it, or further out in the search tree. The commentary introducing proj_as_proj_string() is perhaps suggestive:

Get a PROJ string representation of an object. The returned string is valid while the input obj parameter is valid, and until a next call to proj_as_proj_string() with the same input object.

[rsbivand]

Rev 919 to catch a missed ./-> change in GDAL 2 code found on R-Forge.

[rsbivand]

With sf master I can replicate with PROJ 6.3.0 and GDAL 3.1 master:

  Running ‘crs.R’
  Comparing ‘crs.Rout’ to ‘crs.Rout.save’ ...8c8
<   GDAL Error 1: PROJ: proj_create_from_database: crs not found
---
>   GDAL Error 6: EPSG PCS/GCS code -1 not found in EPSG support files.  Is this a valid EPSG coordinate system?
13c13
<   GDAL Error 1: PROJ: proj_create_from_database: crs not found
---
>   GDAL Error 6: EPSG PCS/GCS code 999999 not found in EPSG support files.  Is this a valid EPSG coordinate system?
15,16c15
< proj_create: unrecognized format / unknown name
< Error : invalid crs: error, reason: generic error of unknown origin
---
> Error : invalid crs: error, reason: no arguments in initialization list
21c20

[edzer]

With sf master and GDAL 3.0.3 and PROJ 6.3.0 released, the problem I saw above no longer happens. I hope that means it's fixed!

[rsbivand]

Now +towgs84= ? https://github.com/r-spatial/sf/issues/1328

[nevilamos]

Hi I am encountering the same error.
Error in CPL_transform(x, crs, aoi, pipeline, reverse) :
OGRCreateCoordinateTransformation() returned NULL: PROJ available?

following change of GDAL version 2.4.2 from ( I think) 2.2.x on ubuntu 18.04, current cran verions of rgdal(1.4-8), and sf(0.9-2)

library(rgdal)
rgdal: version: 1.4-8, (SVN revision 845)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
Path to GDAL shared files: /usr/share/gdal
GDAL binary built with GEOS: TRUE
Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
Path to PROJ.4 shared files: (autodetected)
Linking to sp version: 1.4-1

from Edzer's commnet above if I update GDAL and proj 4 ( not sure how to acheive the latter) then the issue is fixed? is this correct?

thanks

[jcvdav]

Not sure it belongs here, but similar issues were raised in the history of the issue. Happy to post in the correct one.

Having trouble too, but I am on GEOS 3.8.1, GDAL 3.0.4, PROJ 7.0.0. (I am, however, using renv and I've noticed it doesn't play nicely with the rspatial packages that link to projlib and gdal.)

In short, I am trying to st_transform.

Modifying the example in ?st_transform:

library(sf)

p1 = st_point(c(7,52))
p2 = st_point(c(-30,20))
sfc = st_sfc(p1, p2, crs = 4326)
sfc
st_transform(sfc, 54009) # <<------- I'm trying to reproject to Mollweide (https://epsg.io/54009)

And get this error message:

Error in CPL_transform(x, crs, aoi, pipeline, reverse) : 
  crs not found: is it missing?
In addition: Warning message:
In CPL_crs_from_input(x) :
  GDAL Error 1: PROJ: proj_create_from_database: crs not found

FWIW, when I explicitly include the entire proj4string, it does work.

---

Relevant info here:

R version 3.6.3 (2020-02-29)
Platform: x86_64-apple-darwin19.3.0 (64-bit)
Running under: macOS Catalina 10.15.4

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /usr/local/Cellar/openblas/0.3.9/lib/libopenblasp-r0.3.9.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] sf_0.9-3

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.4.6       rstudioapi_0.11    magrittr_1.5       knitr_1.28.3      
 [5] units_0.6-6        R6_2.4.1           rlang_0.4.5        fansi_0.4.1       
 [9] tools_3.6.3        grid_3.6.3         pkgbuild_1.0.6     xfun_0.13         
[13] KernSmooth_2.23-16 e1071_1.7-3        DBI_1.1.0          cli_2.0.2         
[17] withr_2.1.2        class_7.3-15       htmltools_0.4.0    remotes_2.1.1     
[21] yaml_2.2.1         assertthat_0.2.1   rprojroot_1.3-2    digest_0.6.25     
[25] crayon_1.3.4       processx_3.4.2     callr_3.4.3        ps_1.3.2          
[29] curl_4.3           glue_1.4.0         evaluate_0.14      rmarkdown_2.1     
[33] compiler_3.6.3     backports_1.1.6    prettyunits_1.1.1  classInt_0.4-3    
[37] renv_0.9.3  

[Nowosad]

@jcvdav your problem is different. To specify CRS using the ESRI code, you should use:

st_transform(sfc, 'ESRI:54009')

[jcvdav]

Ah, that makes sense. Sorry about that.

This used to work before, but this time I was running it on a new machine, so I assumed it had to do with how I had set up and linked GDAL or PROJ. Was this a recent change in sf?

[edzer]

No, in PROJ.