PROJ6 and GDAL 2.5 are coming

[edzer]

See also #545, the changes this time are much more far reaching, and also involve GDAL interfaces, potentially even axis swapping. This affects both rgdal and sf. @rsbivand and I are pondering about submitting an R Consortium proposal that includes addressing this, and communicating the changes to R users.

Links:

• https://trac.osgeo.org/gdal/wiki/rfc73_proj6_wkt2_srsbarn
• https://gdalbarn.com/#why
• http://blog.cleverelephant.ca/2019/02/proj4-postgis.html

[edzer]

affecting both sf and rgdal:

Step 28/35 : RUN R CMD INSTALL rgdal_*.tar.gz
 ---> Running in d016b7b7357f
* installing to library '/usr/local/lib/R/site-library'
* installing *source* package 'rgdal' ...
configure: R_HOME: /usr/lib/R
configure: CC: gcc -std=gnu99
configure: CXX: g++
configure: C++11 support available
configure: rgdal: 1.3-9
checking for /usr/bin/svnversion... yes
cat: inst/SVN_VERSION: No such file or directory
configure: svn revision: 
checking for gdal-config... /usr/local/bin/gdal-config
checking gdal-config usability... yes
configure: GDAL: 2.4.0
checking C++11 support for GDAL >= 2.3.0... yes
checking GDAL version >= 1.11.4... yes
checking GDAL version <= 2.5.0... yes
checking gdal: linking with --libs only... yes
checking GDAL: gdal-config data directory readable... yes
checking GDAL: /usr/local/share/gdal/pcs.csv readable... yes
configure: pkg-config proj exists, will use it
configure: PROJ version: 6.0.0
checking proj_api.h presence and usability... no
configure: error: proj_api.h not found in standard or given locations.
ERROR: configuration failed for package 'rgdal'
* removing '/usr/local/lib/R/site-library/rgdal'
The command '/bin/sh -c R CMD INSTALL rgdal_*.tar.gz' returned a non-zero code: 1

[rsbivand]

How do I trigger a test run, from R-Forge? If I commit to R-Forge, can this CI see it? I'll be going with -DACCEPT_USE_OF_DEPRECATED_PROJ_API_H for now if the first character in the Proj version string is "6". Is GDAL built against 6.0.0? (no, I think, you'd need the master branch for 2.5.0?)

[edzer]

That seems like the best short term fix.

I used this Dockerfile. I find developing with travis or appveyor as the test platform a pain.

[rsbivand]

rgdal committed to R-Forge passing check on PROJ 6.0.0 and GDAL 2.4.0 built against PROJ 6.0.0. sf is stuck for me at:

checking GDAL: checking whether PROJ is available fur running:... yes
configure: pkg-config proj exists, will use it
checking proj_api.h usability... no
checking proj_api.h presence... no
checking for proj_api.h... no
configure: error: proj_api.h not found in standard or given locations.

which may mean that you need to look at how rgdal detects `proj_api.h the clunky way (own code, not autoconf's code) after setting the define. The file is in the declared location, but it barfs if the define is not set:

#ifndef ACCEPT_USE_OF_DEPRECATED_PROJ_API_H
#error 'To use the proj_api.h you must define the macro ACCEPT_USE_OF_DEPRECATED
_PROJ_API_H'
#endif

[rsbivand]

Yes, PR coming I hope, but same list problems as referenced in OSGeo/PROJ#1306.

[edzer]

I can confirm rgdal now compiles with PROJ6.

[rsbivand]

#1000 and r-spatial/lwgeom#41 made. Needed to bootstrap sf without vignettes, install lwgeom, then build sf and re-install. Haven't fixed lwgeom fully, didn't try to install on PROJ < 6, so maybe missed stuff.

[edzer]

Congrats on hitting #1000!

[rsbivand]

Markus Metz reports on https://lists.osgeo.org/pipermail/grass-dev/2019-March/091539.html

Regarding parsing the traditional PROJ init files, I looked at pj_init.c in
PROJ 5 and earlier and attempted to improve parsing those traditional init
files in GRASS trunk r74224. I guess that R as well as GRASS is trying to
support PROJ versions from 4(.8) onwards.

[edzer]

I created a new branch, PROJH, which uses proj.h instead of proj_api.h. I could get things down to errors like

> demo(meuse, ask = FALSE, echo = FALSE)
> summary(st_as_sf(meuse))
proj_create: init=epsg:/init=IGNF: syntax not supported in non-PROJ4 emulation mode
Error in CPL_proj_is_valid(p4s) : 
  basic_string::_M_construct null not valid
Calls: summary ... valid_proj4string -> structure -> CPL_proj_is_valid -> .Call

indicating that the new PROJ doesn't like the +proj:init=epsg constructs. I don't know whether or how we can set PROJ4 emulation mode.

This branch hard-codes -DPROJH in src/Makevars.in, because autoconf still goes the proj_api.h path.
Tested now with

Linking to GEOS 3.5.1, GDAL 2.4.0, PROJ 6.0.0

[rsbivand]

Does PostGIS use EPSG codes first, rather than using proj4 strings? Can we learn anything from them or GRASS?

[edzer]

PostGIS has its own spatial_ref_sys table

postgis=# select * from spatial_ref_sys limit 3;
 srid | auth_name | auth_srid |                                                                                                                                                                    srtext                                                                                                                                                                    |                                           proj4text                                            
------+-----------+-----------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------------
 3819 | EPSG      |      3819 | GEOGCS["HD1909",DATUM["Hungarian_Datum_1909",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],TOWGS84[595.48,121.69,515.35,4.115,-2.9383,0.853,-3.408],AUTHORITY["EPSG","1024"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","3819"]] | +proj=longlat +ellps=bessel +towgs84=595.48,121.69,515.35,4.115,-2.9383,0.853,-3.408 +no_defs 
 3821 | EPSG      |      3821 | GEOGCS["TWD67",DATUM["Taiwan_Datum_1967",SPHEROID["GRS 1967 Modified",6378160,298.25,AUTHORITY["EPSG","7050"]],AUTHORITY["EPSG","1025"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","3821"]]                                                                 | +proj=longlat +ellps=aust_SA +no_defs 
 3824 | EPSG      |      3824 | GEOGCS["TWD97",DATUM["Taiwan_Datum_1997",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","1026"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","3824"]]                                            | +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
(3 rows)

working with epsg codes through GDAL seems to work, only the +proj:init=epsg breaks.

I haven't looked into liblwgeom, and changes there for dropping proj_api.h, yet.

[rsbivand]

... and the GDAL master IIRC uses PROJ6, so I think must have code to bridge files using proj4 strings. The +init=epsg:# must have raised attention, though in doing the listings in rgdal/src/proj_info6.cpp proj_as_proj_string() going to proj4 strings threw errors which I had to mute by setting a no-op logger.

[edzer]

Given support for proj.h, what should the configure logic be? I'd suggest what PostGIS does:

• if proj_api.h is present, use it
• else if proj.h is present, use it, don't ACCEPT_USE_OF_DEPRECATED_PROJ_API_H
• else break with an error?

[rsbivand]

My safety-first approach for rgdal will be to add a configure flag with-proj-api, defaults "yes" if version < 6, but can be "yes" or "no" if major version == 6 to make checking for changes of behaviour easier. sf may be more adventurous ...

[edzer]

but what will be the default?

[rsbivand]

Now, "yes" where it still works in PROJ 6, (proj.h where it doesn't silently, that is as now) with no change to project(), transform(), etc. List someting like: RGDAL_GetProjectionRef(), RGDAL_checkCRSArgs(), ogrP4S(), p4s_to_wkt(), wkt_to_p4s(), ogrAutoIdentifyEPSG(), project_inv(), RGDAL_project(), transform(), ...

So default "yes" until we know that we are backward-compatible with "no" in PROJ 6, and at that stage flip the default for PROJ 6?

[rsbivand]

Can we tell whether GDAL was built with proj.h from within GDAL? Should we expect both upstream components to use the same headers?

[edzer]

Not that I know of; I'm not sure why this would be needed. I suppose GDAL 2.5 will be proj.h only.

[rsbivand]

Yes, so if need be check across with GDAL version.

[dbaston]

I may be able to help with the API change if you'd like. I'll be at the May OSGeo code sprint with other GDAL, PROJ, and PostGIS devs (others are welcome!) and could devote a day to it then.

[edzer]

Testing against the master branch of GDAL, I now see (both with rgdal and sf)

checking gdal.h usability... yes
checking gdal.h presence... yes
checking for gdal.h... yes
checking GDAL: linking with --libs only... yes
checking GDAL: /usr/local/share/gdal/pcs.csv readable... no
configure: error: pcs.csv not found in GDAL data directory.
ERROR: configuration failed for package 'sf'
* removing '/usr/local/lib/R/site-library/sf'

which was to be expected, right? ("The dreaded ad hoc CSV databases in PROJ_LIB and GDAL_DATA are frustrating for users, pose challenges for developers, and impede interoperability of definitions.")

[edzer]

With these changes, sf passes check in GDAL 2.5.0 (dev) and PROJ 6.

For rgdal, see ac7e4ad#diff-67e997bcfdac55191033d57a16d1408a and this commit. Checking rgdal (after outcommenting the error msg when pcs.csv is missing) gives this

# R CMD check --no-build-vignettes --no-vignettes --no-manuals rgdal_1.4-3.tar.gz 
Warning: unknown option '--no-manuals'
* using log directory '/git/rgdalxx/rgdal.Rcheck'
* using R version 3.4.4 (2018-03-15)
* using platform: x86_64-pc-linux-gnu (64-bit)
* using session charset: ASCII
* using options '--no-vignettes --no-build-vignettes'
* checking for file 'rgdal/DESCRIPTION' ... OK
* this is package 'rgdal' version '1.4-3'
* checking package namespace information ... OK
* checking package dependencies ... OK
* checking if this is a source package ... OK
* checking if there is a namespace ... OK
* checking for executable files ... OK
* checking for hidden files and directories ... OK
* checking for portable file names ... OK
* checking for sufficient/correct file permissions ... OK
* checking whether package 'rgdal' can be installed ... OK
* checking installed package size ... OK
* checking package directory ... OK
* checking DESCRIPTION meta-information ... OK
* checking top-level files ... OK
* checking for left-over files ... OK
* checking index information ... OK
* checking package subdirectories ... OK
* checking R files for non-ASCII characters ... OK
* checking R files for syntax errors ... OK
* checking whether the package can be loaded ... OK
* checking whether the package can be loaded with stated dependencies ... OK
* checking whether the package can be unloaded cleanly ... OK
* checking whether the namespace can be loaded with stated dependencies ... OK
* checking whether the namespace can be unloaded cleanly ... OK
* checking loading without being on the library search path ... OK
* checking dependencies in R code ... OK
* checking S3 generic/method consistency ... OK
* checking replacement functions ... OK
* checking foreign function calls ... OK
* checking R code for possible problems ... OK
* checking Rd files ... OK
* checking Rd metadata ... OK
* checking Rd cross-references ... OK
* checking for missing documentation entries ... OK
* checking for code/documentation mismatches ... OK
* checking Rd \usage sections ... OK
* checking Rd contents ... OK
* checking for unstated dependencies in examples ... OK
* checking contents of 'data' directory ... OK
* checking data for non-ASCII characters ... OK
* checking data for ASCII and uncompressed saves ... OK
* checking line endings in shell scripts ... OK
* checking line endings in C/C++/Fortran sources/headers ... OK
* checking line endings in Makefiles ... OK
* checking compilation flags in Makevars ... OK
* checking for GNU extensions in Makefiles ... OK
* checking for portable use of $(BLAS_LIBS) and $(LAPACK_LIBS) ... OK
* checking compiled code ... OK
* checking files in 'vignettes' ... WARNING
Files in the 'vignettes' directory but no files in 'inst/doc':
  'OGR_shape_encoding.Rnw'
* checking examples ... OK
* checking for unstated dependencies in 'tests' ... OK
* checking tests ...
  Running 'srs_rendering.R'
  Comparing 'srs_rendering.Rout' to 'srs_rendering.Rout.save' ...7c7
< [1] "GDAL 2.5.0dev-4db55e6, released 2019/03/04"
---
> [1] "GDAL 2.4.0, released 2018/12/14"
11,12d10
<                                                                                 kiritimati_primary_roads 
<                                                       "+proj=utm +zone=4 +datum=WGS84 +units=m +no_defs" 
14a13,14
>                                                                                                                                            kiritimati_primary_roads 
>                                                                                                                 "+proj=utm +zone=4 +datum=WGS84 +units=m +no_defs " 
18c18
< "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs" 
---
> "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs " 
22c22
< [1] "+proj=longlat +a=6378137.01 +rf=298.257223563 +towgs84=0,0,0,0,0,0,0 +no_defs"
---
> [1] "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs "
28c28
< [1] "+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.999877420000019 +x_0=600000 +y_0=2200000.00000032 +ellps=clrk80ign +pm=paris +towgs84=-168,-60,320,0,0,0,0 +units=m +no_defs"
---
> [1] "+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.9998774200000193 +x_0=600000 +y_0=2200000.000000325 +a=6378249.2 +b=6356515.000000472 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs "
40,43c40,43
< file:  SP27GTIF.TIF , SRS:  +proj=tmerc +lat_0=36.6666666666667 +lon_0=-88.3333333333333 +k=0.999975 +x_0=152400.30480061 +y_0=0 +ellps=clrk66 +towgs84=-8,160,176,0,0,0,0 +units=us-ft +no_defs 
< file:  cea.tif , SRS:  +proj=cea +lat_ts=33.75 +lon_0=-117.333333333333 +x_0=0 +y_0=0 +datum=NAD27 +units=m +no_defs 
< file:  erdas_spnad83.tif , SRS:  +proj=tmerc +lat_0=30 +lon_0=-82.1666666666667 +k=0.9999 +x_0=200000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs 
< file:  scaleoffset.vrt , SRS:  +proj=utm +zone=18 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
---
> file:  SP27GTIF.TIF , SRS:  +proj=tmerc +lat_0=36.66666666666666 +lon_0=-88.33333333333333 +k=0.9999749999999999 +x_0=152400.3048006096 +y_0=0 +datum=NAD27 +units=us-ft +no_defs  
> file:  cea.tif , SRS:  +proj=cea +lon_0=-117.333333333333 +lat_ts=33.75 +x_0=0 +y_0=0 +datum=NAD27 +units=m +no_defs  
> file:  erdas_spnad83.tif , SRS:  +proj=tmerc +lat_0=30 +lon_0=-82.16666666666667 +k=0.9999 +x_0=200000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs  
> file:  scaleoffset.vrt , SRS:  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs  
  Running 'test_proj.R'
  Comparing 'test_proj.Rout' to 'test_proj.Rout.save' ... OK
  Running 'tests.R'
  Comparing 'tests.Rout' to 'tests.Rout.save' ...7c7
< Error : Null external pointer
---
> Error in getRasterData(x, band = 4) : Null external pointer
21c21
< statistics not supported by this driver 
---
> In GDALinfo(fn) : statistics not supported by this driver
  Running 'tripup.R'
  Comparing 'tripup.Rout' to 'tripup.Rout.save' ... OK
 OK
* checking for unstated dependencies in vignettes ... OK
* checking package vignettes in 'inst/doc' ... WARNING
Package vignette without corresponding PDF/HTML:
   'OGR_shape_encoding.Rnw'

* checking running R code from vignettes ... SKIPPED
* checking re-building of vignette outputs ... SKIPPED
* checking PDF version of manual ... WARNING
LaTeX errors when creating PDF version.
This typically indicates Rd problems.
* checking PDF version of manual without hyperrefs or index ... ERROR
Re-running with no redirection of stdout/stderr.
Hmm ... looks like a package
Error in texi2dvi(file = file, pdf = TRUE, clean = clean, quiet = quiet,  : 
  pdflatex is not available
Error in texi2dvi(file = file, pdf = TRUE, clean = clean, quiet = quiet,  : 
  pdflatex is not available
Error in running tools::texi2pdf()
You may want to clean up by 'rm -rf /tmp/RtmpNsUaqE/Rd2pdf2c582267d012'
* DONE

Status: 1 ERROR, 3 WARNINGs
See
  '/git/rgdalxx/rgdal.Rcheck/00check.log'
for details.

so pretty good (the error comes from pdflatex not installed on the docker image)

[edzer]

@dbaston that is a generous offer, and very very welcome! We're now clearly in "fixing" mode, but it is of course of high interest how also R should deal with e.g. the axis order (hell, or blessing?), reprojections depending on area, and the option (or need?) to use WKT2 rather than PROJ strings. Aligning this with where PostGIS is going sounds like a good idea. When and where is the code sprint?

[dbaston]

@edzer it's in Minneapolis, USA from May 14-17: https://wiki.osgeo.org/wiki/OSGeo_Community_Sprint_2019

[rsbivand]

I'll get back to this next week; I don't have GDAL master installed.

[edzer]

@rsbivand did you look at this, or shall I do so?

[rsbivand]

Please do it, I'm still butchering spdep.

[rsbivand]

Please do it, I'm still butchering spdep. I'm following https://lists.osgeo.org/pipermail/proj/2019-April/008419.html with interest, PROJ 6.1 and proj_normalize_for_visualization() may be a way through ...

[edzer]

Changes here:

rgdal from r-forge now runs against GDAL 2.5.0 (dev), where pcs.csv is missing, with

Step 38/39 : RUN R CMD INSTALL rgdal_*.tar.gz
 ---> Running in 46d728167282
* installing to library '/usr/local/lib/R/site-library'
* installing *source* package 'rgdal' ...
configure: R_HOME: /usr/lib/R
configure: CC: gcc -std=gnu99
configure: CXX: g++
configure: C++11 support available
configure: rgdal: 1.4-3
checking for /usr/bin/svnversion... yes
cat: inst/SVN_VERSION: No such file or directory
configure: svn revision: 
checking for gdal-config... /usr/local/bin/gdal-config
checking gdal-config usability... yes
configure: GDAL: 2.5.0
checking C++11 support for GDAL >= 2.3.0... yes
checking GDAL version >= 1.11.4... yes
checking GDAL version <= 2.5.0... yes
checking gdal: linking with --libs only... yes
checking GDAL: gdal-config data directory readable... yes
checking GDAL: /usr/local/share/gdal/pcs.csv readable... no
configure: pkg-config proj exists, will use it
configure: PROJ version: 6.0.0
configure: Support for PROJ >= 6.0.0 not yet available, deprecated API defined
configure: proj CPP flags: -DACCEPT_USE_OF_DEPRECATED_PROJ_API_H -I/usr/local/include 
checking PROJ header API:... proj_api.h
checking proj_api.h presence and usability... yes
checking PROJ version agreement... yes
checking PROJ version >= 4.8.0... yes
checking PROJ.4: proj.db found and readable... yes
checking PROJ.4: conus found and readable... yes
configure: Package CPP flags:  -I/usr/local/include -I/usr/local/include -DACCEPT_USE_OF_DEPRECATED_PROJ_API_H
configure: Package LIBS:  -L/usr/local/lib -lgdal -lproj
configure: creating ./config.status
config.status: creating src/Makevars

[rsbivand]

OK, thanks!

[GreatEmerald]

I ran into this as well, and I can confirm that the configure of the CRAN version (for R 3.5.3) of sf breaks when reading proj_api.h, the GitHub version works fine. Thanks!

[edzer]

@dbaston are you at the OSGEO sprint? Any concrete plans for working on r-spatial / sf, or any need for inputs?

[dbaston]

Hi @edzer . Yes, I'm at the OSGeo sprint, likely focusing on GEOS for most of the week. Anything you're looking for in sf from the GEOS side?

[edzer]

Maybe you know, or can ask: we now set precision always to 0 (no precision: no rounding taking place), which is quite likely the worst case. Do other OSGEO projects set this by default to some more sensible value, and if yes, to which? Or is this something that (many) GEOS users simply have to learn?

[dbaston]

Do you mean the geometry precision, i.e. as set with GEOSGeom_setPrecision_r? I think most projects do not set this (use the default floating precision model). There is upcoming work to implement snap-rounding in the overlay functions, where a new tolerance setting would be introduced for snapping in the overlays. In that case, a nonzero tolerance should be more robust.

[edzer]

Yes, that one, except that sf doesn't use the GEOS function but uses its own, when converting R doubles to WKB.

What are overlay functions?

[dbaston]

The overlay functions are union / intersection / difference / symdifference. Rounding the inputs can improve robustness, but we don't do it by default because it changes the user's inputs. It's more something that you try if the unrounded operation fails.

[edzer]

OK, thanks - all clear!

[rsbivand]

R CMD check OK for released sf and R-Forge rgdal with PROJ 6.1.0.