sp::proj4string equivalent in the PROJ6 framework?

Question

I am confused about the shift from PROJ4 to PROJ6 in the (relatively) recent versions of the spatial packages of R (and gdal for that matter). I used to use sp::proj4string() to get the projection code / CRS argument from an sp object and use that code further in my scripts. Now I get a warning (which tells me that my projection might be inaccurate). I am obviously doing something wrong. Which function should I use to get the projection code without a warning? I would personally prefer the EPSG code if that level of simplicity was possible...

library(sp)

x <- SpatialPolygons(
            list(Polygons(
                 list(Polygon(cbind(c(1,2,3,4,1), c(3,8,10,12,3)))), 
                 ID = 1)), 
            proj4string = CRS(SRS_string = "EPSG:4326"))

proj4string(x)
#> Warning in proj4string(x): CRS object has comment, which is lost in output
#> [1] "+proj=longlat +datum=WGS84 +no_defs"

Created on 2020-09-18 by the reprex package (v0.3.0)
See also the references and comments in this question

JonasV · Accepted Answer

I was also quite confused about the move away from PROJ Strings and what to use instead until I watched this video of a talk held last year at FOSS4G.
Here's my main takeaways:
Why PROJ and GDAL are moving away from PROJ-strings
With PROJ4 every transformation is done using the "hub-approach". So everything will first be converted to WGS84 and only then transformed to the target projection.
This approach can lead to larger errors in the transformation process than direct transformation between the two projections. At the worst case the needed transformation isn't even possible with this approach.
Since PROJ strings can be ambiguous, direct transformations are often not possible.
An example from the talk:
The old Australian datum GDA94 and the new GDA2020 are both based on WGS84. However one is using the reference frame ITRF2014 at Epoch 2020.0 and the other ITRF92 at Epoch 1994.0. Those differences can't be captured in a PROJ-string (as far as I understood) and thus all projections have the almost exact same PROJ String as WGS84. Converting between both with the hub-approach would lead to no change in the coordinates. While actually Australia has moved 1.8 meters between 1994 and 2020 and the datum conversion should reflect that.
With a direct transformation erros like that would not happen.
What you should use instead
Instead you should ideally use non-ambiguous definitions for projections like EPSG Codes or definitions in WKT2 Format.
However if you do not intend to do datum transformations (or then don't care about a few meters or more of inaccuracy) you can still use PROJ-Strings, but do be aware that best-practice is to use EPSG Codes or WKT2.
How is this implemented in sp
CRS-objects in sp now have an additional field comment in which a WKT2 representation of the CRS is saved.
So the best practice workflow with sp now looks like this (taken and expanded a bit from here):
# Define the CRS using an EPSG Code
x <- CRS(SRS_string='EPSG:4326')

# Display the stored CRS using comment()
cat(comment(x), "n")

# Store the wkt in a variable
wkt <- comment(x)

# Use this to assign the CRS of another sp-object
y <- CRS(SRS_string = wkt)

All in all sf has a lot more convenient functions for accessing and handling crs than sp. You can read more about how sf handles it in the same link given above.

Jeffrey Evans · Answer

Proj4strings are being phased out. Thanks to the sf package and its developers there is a simple trick for converting a proj4string to the WKT format.
Just ignore the datum warnings coming from rgdal and, in fact, you can just turn them off because they do become irksome. Either add options("rgdal_show_exportToProj4_warnings"="none") to your Rprofile.site file or issue it before adding any libraries.
For example with the geographic WGS84 projection, using sf::st_crs("+proj=longlat +datum=WGS84 +no_defs") will result in
Coordinate Reference System:
  User input: +proj=longlat +datum=WGS84 +no_defs 
  wkt:
GEOGCRS["unknown",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6326]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]

From the R-spatial blog R spatial follows GDAL and PROJ development sf developers state:

Use of so-called PROJ4-strings (like +proj=longlat +datum=WGS84) are
discouraged, they no longer offer sufficient description of coordinate
reference systems; use of +init=epsg:XXXX leads to warnings

César Arquero · Answer

That warning doesn't affect your output at all. As the other folks suggested, sf package can be useful for cleaner and simple EPSG coding. Here is your code adapted:
library(sf)
library(mapview)
coords <- list(rbind(c(1,3), c(2,8), c(3,10), c(4,12), c(1,3)))
x2 <- st_sf(st_geometry(st_polygon(x = coords)), crs="EPSG:4326")
x2$ID <- 1
mapview(x2)

sp::proj4string equivalent in the PROJ6 framework?

3 Answers

Why PROJ and GDAL are moving away from PROJ-strings

What you should use instead

How is this implemented in sp

Add your own answers!

Ask a Question