Compare CRS

[etiennebr]

What's the best option to compare crs's ? For instance:

epsg <- "+proj=tmerc +lat_0=0 +lon_0=0 +k=0.9999 +x_0=304800 +y_0=0 +datum=NAD83 +units=m +no_defs"
st_crs(epsg)$proj4string; epsg
# [1] "+proj=tmerc +lat_0=0 +lon_0=0 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
# [1] "+proj=tmerc +lat_0=0 +lon_0=0 +k=0.9999 +x_0=304800 +y_0=0 +datum=NAD83 +units=m +no_defs"

Using == like it's currently the case seems too strict.

epsg <- "+proj=tmerc +lat_0=0 +lon_0=0 +k=0.9999 +x_0=304800 +y_0=0 +datum=NAD83 +units=m +no_defs"
y <- x <- st_crs(epsg)
y$proj4string <- epsg  # this is cheating, but it's for illustration; it happened in a "natural" way via st_read
x == y
# FALSE

I see sp::identicalCRS(), sf doesn't seem to have an equivalent, yet. Any thoughts?

[edzer]

If you don't cheat, x == y is TRUE. What is your point? The equivalent is this -- == is a call to Ops.crs.

[etiennebr]

Sorry if the application wasn't clear. I have shapefiles coming from st_read() with x and I need to merge them with shapefiles coming from st_transform() with y (ref. to the example above). But they won't merge because they have different crs strings. I'd be fine merging them because I know they're the same (nad83 uses GRS80). Is there a way to do it automatically? Apparently proj4 knows it's equivalent because it made the switch.

[edzer]

If you st_transform y into the st_crs of x, do they still have different st_crs?

[etiennebr]

Yes, (sorry for giving context piece by piece, I was attempting to simplify the problem). I did something like:

shapefiles <- lapply(lf, st_read)
crs <- st_crs(shapefiles[[1]])
shapefiles_ident_crs <- lapply(shapefiles, function(.) st_transform(., crs))

Those that didn't need reprojection were left with epsg, but the others have st_crs(epsg).

sapply(shapefiles_ident_crs, st_crs)
# [1] "+proj=tmerc +lat_0=0 +lon_0=-0 +k=0.9999 +x_0=304800 +y_0=0 +datum=NAD83 +units=m +no_defs"                       
# [2] "+proj=tmerc +lat_0=0 +lon_0=-0 +k=0.9999 +x_0=304800 +y_0=0 +datum=NAD83 +units=m +no_defs"                       
# [3] "+proj=tmerc +lat_0=0 +lon_0=-0 +k=0.9999 +x_0=304800 +y_0=0 +datum=NAD83 +units=m +no_defs"                       
# [4] "+proj=tmerc +lat_0=0 +lon_0=-0 +k=0.9999 +x_0=304800 +y_0=0 +datum=NAD83 +units=m +no_defs"                       
# [5] "+proj=tmerc +lat_0=0 +lon_0=-0 +k=0.9999 +x_0=304800 +y_0=0 +datum=NAD83 +units=m +no_defs"                       
# [6] "+proj=tmerc +lat_0=0 +lon_0=-0 +k=0.9999 +x_0=304800 +y_0=0 +datum=NAD83 +units=m +no_defs"                       
# [7] "+proj=tmerc +lat_0=0 +lon_0=-0 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
# [8] "+proj=tmerc +lat_0=0 +lon_0=-0 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
# [9] "+proj=tmerc +lat_0=0 +lon_0=-0 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" 

[edzer]

Are these all semantically identical, or rather the best that GDAL could do for you?

st_transform does not do a hard reset on the object's crs when returned after projection: this is kind of what GDAL tells you it did.

@rsbivand any idea?

[etiennebr]

I haven't looked into it, but maybe the fix could be to run crs from st_read through st_crs() to be sure it returns the same crs.

[edzer]

Yes, my question was whether that is the right thing to do (in this case, and in general). Some of these things cannot be fixed.

[etiennebr]

Somewhat related to #85

[edzer]

I'm not sure; #85 is about us not wanting GDAL to decide that a rubbish datum must be WGS84; this is about GDAL not deciding whether GRS80 and NAD83 are the same datum. Wikipedia writes that there is not a simple yes/no answer to this question.

[edzer]

Closing this now; I think sf shouldn't try to be more clever than the most clever of GDAL and Proj.4, that is how I hope it works now.

[edzer]

I now see:

> library(sf)
Linking to GEOS 3.5.1, GDAL 2.2.0, proj.4 4.9.2, lwgeom 2.3.3 r15473
> crs1 = st_crs("+proj=tmerc +lat_0=0 +lon_0=0 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
> crs2 = st_crs("+proj=tmerc +lat_0=0 +lon_0=0 +k=0.9999 +x_0=304800 +y_0=0 +datum=NAD83 +units=m +no_defs")
> crs1 == crs2
[1] TRUE