Sf: comparing different coordinate reference systems

[jeffreyhanson]

I'm sorry if I'm posting this in the wrong place, or if this covered somewhere else and I've missed it. What is the correct/recommended approach for comparing different coordinate reference systems (crs and/or CRS) to ensure that they correspond to the same system given PROJ (version 7)? In the past, I have used raster::compareCRS. However, this function works for CRS objects and coercing a crs object to a CRS object can now throw a warning for certain coordinate reference systems.

For example, I could use the following code with sf (0.9.2), rgdal (1.4.6) and PROJ (4):

x1 <- st_crs(3755)
x2  <- x1
raster::compareCRS(as(x1, "CRS"), as(x2, "CRS"))

However, with sf (0.9.2), rgdal (1.5.6) and PROJ (7) this now throws a warning:

x1 <- st_crs(3755)
x2  <- x1
raster::compareCRS(as(x1, "CRS"), as(x2, "CRS"))
# Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
# but +towgs84= values preserved

Is there a way to compare coordinate systems without triggering a warning? I suppose I could compare the WKT data. For example, something like this:

x1 <- st_crs(3755)
x2  <- x1
identical(x1$wkt, x2$wkt)

But could this lead to false negatives (i.e. incorrectly determining that two systems are different when they are in fact the same)? For example, if the order of the parameters were shuffled? I suppose that I could simply wrap the as(x, "CRS") in suppressWarnings but this would not be ideal in case the discarded information is relevant?

[jeffreyhanson]

Looking at the source code for rgdal (version 1.5.6) it would appear that the rgdal::compare_CRS function can compare coordinate reference systems using the WKT data (accessed via comment(CRS)). Using this as a reference, I could have a go at writing a similar compare_crs function for the sf package that works with crs objects, and submit a pull request with this function. Would that be useful?

[Nowosad]

@jeffreyhanson how about this approach:

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
x1 <- st_crs(3755)
x2  <- x1
st_crs(x1) == st_crs(x2)
#> [1] TRUE

^{Created on 2020-04-17 by the reprex package (v0.3.0)}

[rsbivand]

It looks as though raster https://github.com/rspatial/raster is reacting to shifting to WKT2, and terra https://github.com/rspatial/terra. It would be helpful to collaborate on this, helping checking and testing.

[jeffreyhanson]

@Nowosad, thank you very much for the suggestion. I just tried playing around with different representations (sorry, I don't know the correct terminology) of the same coordinate reference system to see how this approach works. It looks like it mostly works -- except when comparing st_crs objects defined using integer vs. character objects -- and it correctly deduces that different character representations correspond to the same system (e.g. it handles shuffled PROJ args):

# load package
library(sf)
# Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0

# define EPSG:3577 projargs
projargs <- "+proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132
+x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

# define EPSG:3577 using different approaches
r1 <- st_crs(3577) # integer 
r2 <- st_crs(projargs) # proj4 args
r3 <- st_crs(paste("+init=epsg:3577", projargs)) # proj args with epsg code
r4 <- st_crs(paste(sample(strsplit(projargs, " ")[[1]]),
                          collapse = " ")) # shuffle proj args

r1 == r2
# FALSE
r1 == r3
# FALSE
r1 == r4
# FALSE
r2 == r3
# TRUE
r2 == r4
# TRUE
r3 == r4
# TRUE

So, moving forward, I'll just use this approach for comparing crs and CRS objects. Thank you!

@rsbivand, that's an excellent point - I didn't think about checking the developmental version of the raster package to see if it was being updated to handle this. It doesn't look like compareCRS has been updated yet (https://github.com/rspatial/raster/blob/master/R/compareCRS.R, unless I'm missing a more recent branch) - but I'll try to remember to keep an eye out in the future. Thanks!

I think you've both addressed my question, so I'll close this issue now. Thanks so much again!

[jeffreyhanson]

Also, sorry I forgot to mention this in my previous post, I'm still happy to try writing a compare_crs function that interfaces with PROJ if it would be useful - just let let me know.

[edzer]

Just for completeness: we're being recommended to stop using proj4strings in general, and definitely stop using +proj=init:xxxx - type proj4strings (and use EPSG:xxxx instead).

[rsbivand]

For the archaeologists among us, https://trac.osgeo.org/gdal/ticket/4880 and https://lists.osgeo.org/pipermail/gdal-dev/2012-November/034558.html were preceded by https://trac.osgeo.org/gdal/ticket/3450 and https://lists.osgeo.org/pipermail/grass-user/2010-March/055146.html. So we have known (for the definition of known - tried to ignore) for ten years that PROJ strings were fragile when datum definitions were involved.