3

I am not a GIS expert, and this change caused quite some errors in our analysis. So my question is, why this was changed: From rgdal 1.4.8 to 1.5.16, the +towgs84=0,0,0 disappeard in the output of sp::CRS.

I hope for an answer, which is comprehensive to a non-GIS-expert.

Reproducible example if you have NO sp and NO rgdal package installed:

packageurl <- "http://cran.r-project.org/src/contrib/Archive/sp/sp_1.4-1.tar.gz"
install.packages(packageurl, repos=NULL, type="source")
 
packageurl <- "http://cran.r-project.org/src/contrib/Archive/rgdal/rgdal_1.4-8.tar.gz"
install.packages(packageurl, repos=NULL, type="source")

packageVersion("sp")
# [1] ‘1.4.1’
packageVersion("rgdal")
# [1] ‘1.4.8’

sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
# CRS arguments:
#   +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 

So everything works, that is in our case

spdf <- sp::SpatialPointsDataFrame(
  coords = mission_coords[, c(2,1)],
  data = mission_coords,
  proj4string = sp::CRS(
    "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
df_in <- sp::over(spdf , departement_boundaries , fn = NULL)

are ok. (I have problem to make a small reproducible example - the departement_boundaries from Zurich are open source, but its a huge string from dput.)

Now, update both packages

install.packages("rgdal")
# also installing the dependency ‘sp’
packageVersion("sp")
# [1] ‘1.4.2’
packageVersion("rgdal")
# [1] ‘1.5.16’

and in any case +towgs84=0,0,0 is missing now:

sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs",
        doCheckCRSArgs = TRUE)
# CRS arguments: +proj=longlat +datum=WGS84 +no_defs 
sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs",
        doCheckCRSArgs = FALSE) # ... +ellps=WGS84 is added...(?)
# CRS arguments: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs 

Any ideas for an explanation to a GIS-newbie?


@EdzerPebesma’s answer helped me find R spatial follows GDAL and PROJ development to get information:

x <- sp::CRS("+init=epsg:4326")
cat(comment(x), "\n")

One comment from that site is really unclear to me (perhaps how this relates to +towgs84=?):

  • rgdal generates many warnings about discarded PROJ string keys,
    -> See e.g. Warnung in showSRID(uprojargs, format = "PROJ", multiline = "NO") Discarded datum WGS_1984 in CRS definition, but +towgs84= values preserved from here?
  • intended to alert package maintainers and script authors to the need to review code.
    -> To me it is not clear, what I have to do, see the link above.
  • It is particularly egregious to assign to the CRS object projargs slot directly, and this is unfortunately seem in much code in packages.
    -> Again, what should I do instead?
6

The underlying library, PROJ, has changed dramatically over the past few years, and this is a consequence. The time has come to stop using proj.4 strings. See also this blog post https://www.r-spatial.org/r/2020/03/17/wkt.html which announced upcoming changes to r-spatial.

| improve this answer | |
  • 2
    This is quite exciting and I have been rewriting code to accommodate this change to WKT. But, how mature is this new deployment of PROJ/GDAL across other spatial packages (ie., sp, raster)? Can we just start retooling our packages to incorporate WKT projections with any spatial object and projection operation outside of sf >0.9? Are there helper functions such as sp::is.projected and ways to compare projections for validation that will work with the new format? – Jeffrey Evans Aug 14 at 15:52
  • 1
    Perhaps, you also know the answer here? To me it is really unclear, what I have to do, in order to improve my packages. See my edit after I read your link. (I am struggling with the implications of "stop using proj.4 strings".) – Christoph Aug 14 at 17:16
  • To Jeffrey: sp and rgdal follow this (actually rgdal has pioneered this); the raster author has been informed, but I believe raster still uses on proj4 strings to define CRSs. Not sure what terra does. Please follow r-sig-geo for news on sf, sp and rgdal. – Edzer Pebesma Aug 15 at 16:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.