Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

PROJ 6.3.0 dropping ellps in P4 string #1231

Closed
rsbivand opened this issue Dec 31, 2019 · 25 comments
Closed

PROJ 6.3.0 dropping ellps in P4 string #1231

rsbivand opened this issue Dec 31, 2019 · 25 comments

Comments

@rsbivand
Copy link
Member

rsbivand commented Dec 31, 2019

More #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
Copy link
Member Author

rsbivand commented Jan 1, 2020

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
Copy link
Member Author

rsbivand commented Jan 1, 2020

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

PROJ6_GDAL3_test.txt

@rsbivand
Copy link
Member Author

rsbivand commented Jan 1, 2020

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
Copy link
Member Author

rsbivand commented Jan 1, 2020

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

@rsbivand
Copy link
Member Author

rsbivand commented Jan 1, 2020

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
Copy link
Member Author

rsbivand commented Jan 1, 2020

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
Copy link
Member Author

rsbivand commented Jan 1, 2020

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
Copy link
Member Author

rsbivand commented Jan 1, 2020

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
Copy link
Member

edzer commented Jan 1, 2020

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

@rsbivand
Copy link
Member Author

rsbivand commented Jan 1, 2020 via email

@edzer
Copy link
Member

edzer commented Jan 1, 2020

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
Copy link
Member Author

rsbivand commented Jan 1, 2020 via email

@edzer
Copy link
Member

edzer commented Jan 1, 2020

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
Copy link
Member Author

rsbivand commented Jan 1, 2020 via email

@edzer
Copy link
Member

edzer commented Jan 1, 2020

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
Copy link
Member Author

rsbivand commented Jan 2, 2020

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
Copy link
Member Author

rsbivand commented Jan 2, 2020

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

@rsbivand
Copy link
Member Author

rsbivand commented Jan 2, 2020

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
Copy link
Member

edzer commented Jan 13, 2020

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
Copy link
Member Author

Now +towgs84= ? #1328

@nevilamos
Copy link

nevilamos commented Apr 16, 2020

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
Copy link

jcvdav commented Apr 23, 2020

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
Copy link
Contributor

Nowosad commented Apr 24, 2020

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

st_transform(sfc, 'ESRI:54009')

@jcvdav
Copy link

jcvdav commented Apr 24, 2020

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
Copy link
Member

edzer commented Apr 24, 2020

No, in PROJ.

@edzer edzer closed this as completed Apr 24, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

5 participants