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

Answer 1

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

Check out Mikko's answer below for a more up to date workflow

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.

Answer 2

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

Answer 3

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)

Answer 4

I'll add this as an answer (was originally meant as a comment). Perhaps it helps someone since there is and has been a lot of confusion about this. It expands the answers above. There is another helpful question with answers.

In short, to print the wkt comment and avoid the warnings with PROJ6 CRS objects, use sp::wkt instead of sp::proj4string:

sp::wkt(x)
[1] "GEOGCRS[\"WGS 84 (with axis order normalized for visualization)\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ID[\"EPSG\",6326]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433],\n        ID[\"EPSG\",8901]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]]]"

To assign a CRS to a SpatialPolygons object, which does not have a CRS, use the sp::proj4string function even though you are assigning wkt (= well-known text, PROJ6 CRS format). Note how you can now omit the SRS_string argument within sp::CRS().

library(sp)

x <- SpatialPolygons(
  list(Polygons(
    list(Polygon(cbind(c(1,2,3,4,1), c(3,8,10,12,3)))), 
    ID = 1))
)

proj4string(x) <- CRS("EPSG:4326")

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

wkt(x) # print the PROJ6 wkt comment
#> [1] "GEOGCRS[\"WGS 84 (with axis order normalized for visualization)\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ID[\"EPSG\",6326]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433],\n        ID[\"EPSG\",8901]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]]]"

To borrow the PROJ6 from another SpatialPolygons object, you can either use only proj4string or use the combination of proj4string and wkt:

y <- SpatialPolygons(
  list(Polygons(
    list(Polygon(cbind(c(1,2,3,4,1), c(3,8,10,12,3)))), 
    ID = 1))
)

proj4string(y) <- proj4string(x) # proj4string(y) <- wkt(x) would also work
# but wkt(y) <- wkt(x) would not because wkt() has no assingment.
# proj4string(y) <- CRS(wkt(x)) also works as shown by JonasV

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

wkt(y)
#> [1] "GEOGCRS[\"WGS 84 (with axis order normalized for visualization)\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ID[\"EPSG\",6326]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433],\n        ID[\"EPSG\",8901]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]]]"

The sp::proj4string() should be understood as a function name. It does not necessarily have to do with PROJ4 any longer if it is used in the assignment. The good news is that you only need to change the "inside" from "+init=epsg:4326" to "EPSG:4326" in the sp::CRS() function and things should work as previously as long as you do not use the sp::proj4string() function to retrieve PROJ6 (wkt) information. If you do, you'll get a warning, but the comment is actually not lost anyway in the assignment (as demonstrated by the proj4string(y) <- proj4string(x) example). It is only lost if you use proj4string() to print CRS information:

comment(proj4string(x))
#> NULL

EDIT: What if you want both, the string and comment? Use raster::crs(). This can be useful to keep it PROJ6 through a function/script.

raster::crs(x)
#> CRS arguments: +proj=longlat +datum=WGS84 +no_defs 

comment(raster::crs(x))
[1] "GEOGCRS[\"WGS 84 (with axis order normalized for visualization)\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ID[\"EPSG\",6326]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433],\n        ID[\"EPSG\",8901]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]]]"

^{Created on 2021-05-27 by the reprex package (v2.0.0)}