Sf: CRS with equivalent EPSG codes not seen as equivalent

[dbaston]

Would it make sense for the following two CRS objects to be seen as equivalent?

> st_crs(r)
Coordinate Reference System:
  EPSG: 4326 
  proj4string: "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs"
> st_crs(all)
Coordinate Reference System:
  EPSG: 4326 
  proj4string: "+proj=longlat +datum=WGS84 +no_defs"
> st_crs(all) == st_crs(r)
[1] FALSE

Related: #180

[edzer]

Maybe yes, and that's why sf pushes this problem all the way down to GDAL's OGRSpatialReference::IsSame(). But it doesn't try to be any smarter than that.

[dbaston]

Maybe I can help with a doc clarification. I think the description of equality is incorrect:
https://github.com/r-spatial/sf/blob/master/R/crs.R#L46
(the EPSG code is ignored in checking equality, and the equivalent-but-different proj strings are in fact considered equivalent)

Would it be accurate to say that the proj4 string is the "source of truth" in sf and that, while sf tracks an EPSG code, it is not used?

[edzer]

If the alternative to "document it" would be "fix it", what would you propose?

[edzer]

The wkt2 branch spits this out:

> st_crs("+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")
Coordinate Reference System:
  EPSG: 4326 
  proj4string: "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs"
  wkt2:
BOUNDCRS[
    SOURCECRS[
        GEOGCRS["unknown",
            DATUM["Unknown based on WGS84 ellipsoid",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["longitude",east,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["latitude",north,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["latitude",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["longitude",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],
    ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
        METHOD["Position Vector transformation (geog2D domain)",
            ID["EPSG",9606]],
        PARAMETER["X-axis translation",0,
            ID["EPSG",8605]],
        PARAMETER["Y-axis translation",0,
            ID["EPSG",8606]],
        PARAMETER["Z-axis translation",0,
            ID["EPSG",8607]],
        PARAMETER["X-axis rotation",0,
            ID["EPSG",8608]],
        PARAMETER["Y-axis rotation",0,
            ID["EPSG",8609]],
        PARAMETER["Z-axis rotation",0,
            ID["EPSG",8610]],
        PARAMETER["Scale difference",1,
            ID["EPSG",8611]]]]
> st_crs("+proj=longlat +datum=WGS84 +no_defs")
Coordinate Reference System:
  EPSG: 4326 
  proj4string: "+proj=longlat +datum=WGS84 +no_defs"
  wkt2:
GEOGCRS["unknown",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

suggesting (to me at least) that a different thing is meant with the two proj4 strings.

[rsbivand]

Not to mention the insecurity of the lack of epoch definition, which 4326 doesn't help with. It isn't obvious how to move forward from WGS84, really.

[edzer]

But I agree with @dbaston that the current approach needs better documentation.

[dbaston]

If the alternative to "document it" would be "fix it", what would you propose?

The current approach makes sense to me, I think it just needs to be documented more clearly. From my non-expert perspective it is odd that GDAL returns 4326 from both proj strings yet says they are not the same. But that's on GDAL, not sf.

[dbaston]

Should this remain open until the docs are changed? The docs say that

As a consequence, equivalent but different proj4strings, e.g. \code{ "+proj=longlat +datum=WGS84" } and \code{ "+datum=WGS84 +proj=longlat" } are considered different.

and yet

st_crs("+proj=longlat +datum=WGS84") == st_crs( "+datum=WGS84 +proj=longlat" )

returns TRUE