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 CRS handling changes impact workflows using +datum= (and other fun) #28

Open
rsbivand opened this issue Sep 11, 2019 · 58 comments
Open

Comments

@rsbivand
Copy link
Member

rsbivand commented Sep 11, 2019

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

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

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

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

@MichaelChirico
Copy link

MichaelChirico commented Nov 5, 2019

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?

Screenshot 2019-11-05 at 1 21 48 PM

@Nowosad
Copy link
Contributor

Nowosad commented Nov 5, 2019

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

@Nowosad great! Is this hosted on dockerhub?

@Nowosad
Copy link
Contributor

Nowosad commented Nov 7, 2019

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

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

@edzer
Copy link
Member

edzer commented Nov 12, 2019

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

rsbivand commented Nov 12, 2019 via email

@rsbivand
Copy link
Member Author

rsbivand commented Nov 13, 2019

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

edzer commented Dec 6, 2019

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

edzer commented Dec 6, 2019

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

rsbivand commented Dec 6, 2019

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

edzer commented Dec 6, 2019

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

edzer commented Dec 6, 2019

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

rsbivand commented Dec 6, 2019

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

edzer commented Dec 6, 2019

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

edzer commented Dec 9, 2019

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

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

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

lbusett commented Dec 26, 2019

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

rsbivand commented Feb 5, 2020

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

rsbivand commented Mar 5, 2020

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

Now +towgs84= r-spatial/sf#1328

netbsd-srcmastr pushed a commit to NetBSD/pkgsrc that referenced this issue Jul 23, 2020
Upstream changes:
Changes in version 1.4-1 (2020-xx-yy)

    warn on NULL projargs in CRS(); edzer/sp#74

Changes in version 1.4-0 (2020-02-21)

    prepare for new (>= 1.5.1) rgdal, which creates and listens to a comments() field of a CRS object carrying a WKT representation of a CRS rather than the proj4string; @rsb, edzer/sp#67 and edzer/sp#69 ; for more info see e.g. edzer/sp#68 and r-spatial/discuss#28

Changes in version 1.3-2 (2019-11-07)

    fix length > 1 in coercion to logical error; #54, #60

    add is.na method for CRS objects
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

8 participants