PROJ 6 CRS handling changes impact workflows using +datum= (and other fun)

[rsbivand]

Ongoing changes in the representation of coordinate reference systems in the PROJ, external software used by sf, sp and other R packages, will impact R packages and workflows using these packages. In particular, the changes affect transformation between coordinate reference systems, as for example shown in this presentation at the OpenGeoHub Summer School 2019.

• Invite users/developers to contribute cases

• Add other examples and links, also from other geospatial software.

GDAL RFC 73: https://gdal.org/development/rfc/rfc73_proj6_wkt2_srsbarn.html

Notes from Spatialite

PROJ wiki proj.h adoption status

GRASS-dev thread and PR

PostGIS blog posting

cartography CRS issues

geojsonio geojsonio seems vulnerable to PROJ 6

MODIS getTile() sp subset or raster::crop issues with PROJ 6.1 and GDAL 3.0

MODIStsp PROJ 6 failures in MODIStsp

rangeMapper PROJ 6 failures in rangeMapper

rdefra PROJ 6 failures in rdefra

rnrfa PROJ 6 failures in rnfra

rshapemapper PROJ 6 problems?

Incomplete list of packages with CRAN issues after PROJ 6.2.0 installed on some test machines:

https://cran.r-project.org/web/checks/check_results_foieGras.html
https://cran.r-project.org/web/checks/check_results_plotKML.html
https://cran.r-project.org/web/checks/check_results_rangeMapper.html
https://cran.r-project.org/web/checks/check_results_rnrfa.html
https://cran.r-project.org/web/checks/check_results_sf.html
https://cran.r-project.org/web/checks/check_results_vapour.html

• Provide recommendations for checking/testing PROJ versions in use and expected transformation outcomes.

See also r-spatial/sf#1146 for sf and rgdal reprex with PROJ 4.9.3/GDAL 2.4.2 output from the same file contrasted with PROJ 6.2.0/GDAL 3.0.1.

This issue is based on https://github.com/rsbivand/geostat19_talk; https://rsbivand.github.io/geostat19_talk/Discuss_issue_drafts.html

[rsbivand]

Relevant CDN for serving grid files discussion on https://twitter.com/howardbutler/status/1171778886646022145?s=20. On the client side related to thread: https://lists.osgeo.org/pipermail/proj/2019-September/008825.html.

[rsbivand]

Adding files (scripts, output) map_package_analyses.zip for packages with reverse dependencies (depends, imports, suggests) on rgdal using pkgapi map_package() to try to see which packages use spTransform(), project(), and raster projectExtent() or projectRaster(). About 80 packages of almost 300 in total seem to use projection mechanisms on the sp side. The scripts may be used to set up similar reports on the sf side. Other columns report CRS() and other vulnerable points.

[rsbivand]

See also suggestions on https://github.com/pyproj4/pyproj; https://pyproj4.github.io/pyproj/stable/index.html.

[MichaelChirico]

Hi Roger, much appreciate your reaching out and pointing to all the resources.

However, I feel a bit overwhelmed by it all and am struggling to turn this into action for how to get my package (geohashTools) working with PROJ>=6.

I see this entry in map_package_analyses/df_all.csv for my package:

   raster_projectExtent raster_projectRaster rgdal_checkCRSargs rgdal_CRSargs rgdal_make_EPSG
1:                    0                    0                  0             0               0
   rgdal_project rgdal_spTransform sp_CRS sp_proj4string sp_spTransform sum
1:             0                 0      1              2              1   4

How can I use this to proceed?

Also, is it safe to use rocker/geospatial for testing (i.e., does that image have the proper system libraries in place already)?

Seems rocker/geospatial is still using PROJ 4.9.3, can you recommend any Docker image to use for testing?

[Nowosad]

@MichaelChirico for the purpose of testing PROJ6 in the Geocomputation with R book, I've created a new docker image based on rocker:

FROM rocker/rstudio:3.6.1

RUN apt-get update \
  && apt-get install -y --no-install-recommends \
	gdb \
	git \
	libcairo2-dev \
	libcurl4-openssl-dev \
	libexpat1-dev \
	libpq-dev \
	libsqlite3-dev \
	libudunits2-dev \
	make \
	pandoc \
	qpdf \
	r-base-dev \
        sqlite3 \
	subversion \
	valgrind \
	vim \
	wget

RUN apt-get install -y --no-install-recommends \
	libv8-3.14-dev \
	libjq-dev \
	libprotobuf-dev \
	libxml2-dev \
	libprotobuf-dev \
	protobuf-compiler \
	unixodbc-dev \
	libssh2-1-dev \
	libgit2-dev \
	libnetcdf-dev \
	locales \
	libssl-dev

RUN locale-gen en_US.UTF-8

ENV PROJ_VERSION=6.2.0
# ENV LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH

RUN wget http://download.osgeo.org/proj/proj-${PROJ_VERSION}.tar.gz \
  && tar zxf proj-*tar.gz \
  && cd proj* \
  && ./configure \
  && make \
  && make install \
  && cd .. \
  && ldconfig

# install proj-datumgrid:
RUN cd /usr/local/share/proj \
  && wget http://download.osgeo.org/proj/proj-datumgrid-latest.zip \
  && unzip -o proj-datumgrid*zip \
  && rm proj-datumgrid*zip \
  && wget https://download.osgeo.org/proj/proj-datumgrid-europe-latest.zip \
  && unzip -o proj-datumgrid*zip \
  && rm proj-datumgrid*zip \
  && cd -

# GDAL:

ENV GDAL_VERSION=3.0.1
ENV GDAL_VERSION_NAME=3.0.1

RUN wget http://download.osgeo.org/gdal/${GDAL_VERSION}/gdal-${GDAL_VERSION_NAME}.tar.gz \
  && tar -xf gdal-${GDAL_VERSION_NAME}.tar.gz \
  && cd gdal* \
  && ./configure \
  && make \
  && make install \
  && cd .. \
  && ldconfig

#RUN git clone --depth 1 https://github.com/OSGeo/gdal.git
#RUN cd gdal/gdal \
#  && ls -l \
#  && ./configure \
#  && make \
#  && make install \
#  && cd .. \
#  && ldconfig

# GEOS:
ENV GEOS_VERSION=3.7.2
#
RUN wget http://download.osgeo.org/geos/geos-${GEOS_VERSION}.tar.bz2 \
  && bzip2 -d geos-*bz2 \
  && tar xf geos*tar \
  && cd geos* \
  && ./configure \
  && make \
  && make install \
  && cd .. \
  && ldconfig

RUN install2.r --error \
    remotes

RUN R -e "remotes::install_github('geocompr/geocompkg')"

[MichaelChirico]

@Nowosad great! Is this hosted on dockerhub?

[Nowosad]

It was not... until today;)
There are two versions:

1. https://hub.docker.com/r/jakubnowosad/rspatial_proj6 - with all external libraries (PROJ6, GDAL3, GEOS)
2. https://hub.docker.com/r/jakubnowosad/geocompr_proj6 - the same as above + many R packages for geocomputation (https://github.com/geocompr/geocompkg/blob/master/DESCRIPTION)

[rsbivand]

This is the WIP-vignette from R-Forge/rgdal 1.5-1 rev 888. Comments welcome!
PROJ6_GDAL3.zip

[edzer]

DO make the WKT representations readable using cat and pretty, as in

> library(sf)
Linking to GEOS 3.7.1, GDAL 2.4.2, PROJ 5.2.0
> st_as_text(st_crs(4326))
[1] "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]"
> st_as_text(st_crs(4326), pretty = TRUE)
[1] "GEOGCS[\"WGS 84\",\n    DATUM[\"WGS_1984\",\n        SPHEROID[\"WGS 84\",6378137,298.257223563,\n            AUTHORITY[\"EPSG\",\"7030\"]],\n        AUTHORITY[\"EPSG\",\"6326\"]],\n    PRIMEM[\"Greenwich\",0,\n        AUTHORITY[\"EPSG\",\"8901\"]],\n    UNIT[\"degree\",0.0174532925199433,\n        AUTHORITY[\"EPSG\",\"9122\"]],\n    AUTHORITY[\"EPSG\",\"4326\"]]"
> cat(st_as_text(st_crs(4326), pretty = TRUE), "\n")
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]
>

(misses the final newline, but nevertheless)

[rsbivand]

Suggestion: do it yourself and commit to svn. I'm busy, and don't think appearances matter enough. Roger Bivand Falsensvei 32 5063 Bergen tir. 12. nov. 2019, 19.22 skrev Edzer Pebesma <notifications@github.com>:

…

DO make the WKT representations readable using cat and pretty, as in > library(sf) Linking to GEOS 3.7.1, GDAL 2.4.2, PROJ 5.2.0 > st_as_text(st_crs(4326)) [1] "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]" > st_as_text(st_crs(4326), pretty = TRUE) [1] "GEOGCS[\"WGS 84\",\n DATUM[\"WGS_1984\",\n SPHEROID[\"WGS 84\",6378137,298.257223563,\n AUTHORITY[\"EPSG\",\"7030\"]],\n AUTHORITY[\"EPSG\",\"6326\"]],\n PRIMEM[\"Greenwich\",0,\n AUTHORITY[\"EPSG\",\"8901\"]],\n UNIT[\"degree\",0.0174532925199433,\n AUTHORITY[\"EPSG\",\"9122\"]],\n AUTHORITY[\"EPSG\",\"4326\"]]" > cat(st_as_text(st_crs(4326), pretty = TRUE), "\n") GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4326"]] > (misses the final newline, but nevertheless) — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#28?email_source=notifications&email_token=ACNZ3BELABVSHPA6HMBLTRDQTLX6JA5CNFSM4IVRDCLKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOED3JLTY#issuecomment-553031119>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACNZ3BFJH27CCQ6XG6TOD3TQTLX6JANCNFSM4IVRDCLA> .

[rsbivand]

All I need is an indication if someone else is attacking this problem and can see obviously better resolutions, or indications (thanks @Nowosad !) of checks on the current state of play. I do not need comments about aesthetics, they do not matter for the task in hand. I'll wrap the output, but I will not change to multiline out from exportToWkt(), because the idea is to present a string that is machine readable. Unless it can be shown that WKT needs to be human-readable (worse editable), I'll keep then as single line strings.

[edzer]

It is good to know that this swapping doesn't happen when a CRS is initialised with a PROJ string. As this is up to now the only acceptable way, this problem shouldn't affect any r-spatial work done so far.

[edzer]

Inserting, in rgdal, pkg/src/proj_info6.cpp, line 358, the line

pj_transform = proj_normalize_for_visualization(ctx, pj_transform); // EJP

gives

coordinates(sp::spTransform(mypoint, CRS(SRS_string = "EPSG:4326")))
#      coords.x1 coords.x2
# [1,]  7.199383  3.617083

but the more important discussion is whether we want this, or want to be OGC compliant. As of: what is wrong with people wanting to do this to themselves. Or default to one and allow for the other?

[rsbivand]

Yes, that is the problem. How can we adopt a default (visualization axis order) but still permit users who need OGC compliant output to achieve that?

[edzer]

One option would be a global setting. It feels more elegant though to keep data in authority (wkt2/ogc) order, and do an axis swap on the fly when it is needed, e.g.

• before plotting
• before (and after) doing geometry operations

I will run some experiments in sf with setting SetAxisMappingStrategy(OAMS_AUTHORITY_COMPLIANT); by default.

[edzer]

Mmm for that to work, datasets like nc would need to be rewritten to have lat long coordinates first, or undergo an axis swap at read time. But based on what? Then, the bounding box logic breaks (what was xmin?), but after that projecting works fine! Feels like a problem that will not easily go away.

[rsbivand]

Perhaps restrict sp workflows to be GIS/visualization mode, because that is how the classes were designed? We can't rewrite legacy data sets that themselves may not encode/report axis order because we need to keep things running for both GDAL2/PROJ5 and GDAL3/PROJ6 at the same time (think CentOs).

[edzer]

Sounds fair enough. R CMD check sf_*gz has btw a test where proj_normalize_for_visualization as used above returned NULL, which I didn't understand, but which caused a segfault when left unchecked.

[edzer]

In case others want to experiment with sf and handling objects with authority compliant axes order, see r-spatial/sf#1146 : branch wkt2 in the sf github repo.

[rsbivand]

Video of talk on this problem on: https://www.youtube.com/playlist?list=PLXUoTpMa_9s10NVk4dBQljNOaOXAOhcE0, third in playlist; presentation at https://rsbivand.github.io/ECS530_h19/ECS530_III.html.

[rsbivand]

@ranghetti , @lbusett : please see rgdal SVN revision >= 903 (905 current) and scroll down in http://rgdal.r-forge.r-project.org/articles/PROJ6_GDAL3.html to Axis swapping to see a version of the Rpubs round-trip issue (when the served version updates). I've added code to address this, but using list_coordOps() already permitted the creation of a transformation object that would output visualization order objects. In PROJ, it is the transformation object that loses/gains a step to handle coordinate swapping, so setting the GDAL SRS axis policy on a CRS may or may not work (when using PROJ code to transform). So far, I hope, so good ... and thanks for the nice present!

[lbusett]

@rsbivand Looks good, thanks! A bit short on time right now, but I'll try to have a better look in the next days.

PS: sorry for the present... ;-(

[rsbivand]

New fun: see https://lists.osgeo.org/pipermail/proj/2020-February/009325.html. Briefly, when a grid is dowloaded and used in a coordinate operation pipeline, coordinates outside the grid's area of interest may be returned as missing.

[rsbivand]

Video https://www.youtube.com/watch?v=D4-roPsMz48 and slides https://github.com/rsbivand/celebRation20_files from a talk at celebRation 2020 in Copenhagen published. The PROJ part starts about 37:20; enjoy (if that is the right expression).

[rsbivand]

Now +towgs84= r-spatial/sf#1328