How to switch to PROJ6 in r

user_5892805 0
Viewed 94

I am trying to switch away from rgdal and PROJ4 since both are deprecated.  I try to assign a coordinate reference system using sf but continue to get an error about the crs still being PROJ4.  

    pacman::p_load(sf, raster)
    cell.coords <- data.frame(lon = rep(c(2, 3), 2), lat =
       rep(c(1, 2), each = 2), cell.id = 10:13)
    blank.grid <- st_as_sf(rasterToPolygons(rasterFromXYZ(
       cell.coords[ , c("lon", "lat", "cell.id")])))
 
    # Try, in turn, each of the following.
    # First three successfully show the crs when blank.grid is called,
    # but crs(blank.grid) includes "Deprecated Proj.4 representation"
    wgs84 <- 4326
    # next is based on https://gis.stackexchange.com/questions/372692/how-should-i-handle-crs-properly-after-the-major-change-in-proj-library
    wgs84 <- "OGC:CRS84"
    # next is based on https://gis.stackexchange.com/questions/385114/how-to-update-crs-r-object-from-proj4-to-proj6
    wgs84 <- "+proj=lcc +lat_1=49 +lat_2=77 +lat_0=0 +lon_0=-95 +x_0=0 +y_0=0 +ellps=GRS80 +datum=WGS84 +units=m +no_defs"
 
    # For the rest, blank.grid returns CRS:  NA
    wgs84 <- "4326+datum=WGS84"
    wgs84 <- "+init=epsg:4326+datum=WGS84"
    wgs84 <- "EPSG:4326"
 
    # Apply and check each option with
    st_crs(blank.grid) <- wgs84
    blank.grid
    crs(blank.grid)

  I realize we are referred (https://gis.stackexchange.com/questions/372692/how-should-i-handle-crs-properly-after-the-major-change-in-proj-library) to the technical references https://cran.r-project.org/web/packages/rgdal/vignettes/CRS_projections_transformations.html and/or http://rgdal.r-forge.r-project.org/articles/PROJ6_GDAL3.html.  Both are long documents I don’t understand, and involve a level of detail I hope is not required for using simple rasters.

An optimal solution would be a single command to create blank.grid, with crs as an argument. Using terra rather than raster would be best, since I think I read that raster is likewise not long for this world.

proj terra r-raster sf r
1 Answers

You can use rast to create a raster template (blank raster). There are different ways to specify the crs (authority:code, PROJ, or WKT). For example:

library(terra)
# authority:code
r1 <- rast(crs="EPSG:4326")
r1
#class       : SpatRaster 
#dimensions  : 180, 360, 1  (nrow, ncol, nlyr)
#resolution  : 1, 1  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (EPSG:4326) 

# PROJ
r2 <- rast(crs="+proj=longlat")
r2
#class       : SpatRaster 
#dimensions  : 180, 360, 1  (nrow, ncol, nlyr)
#resolution  : 1, 1  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs

r1 and r2 have the same coordinate reference system for most intents and purposes (longitude/latitude, WGS84), but note the difference in how they are printed above. For an EPSG code a short description like "lon/lat" is typically available.

But these strings you tried: "4326+datum=WGS84" and "+init=epsg:4326+datum=WGS84" are not valid because they mix EPSG and PROJ. You can either use a single EPSG code, or specify parameters via the PROJ notation. The one exception is using the PROJ notation like this: "+init=epsg:4326" (but there is no reason to use that notation anymore as you can instead use "epsg:4326").

With sf, you can also specify the EPSG code as the number only "4326", without the authority prefix, but terra does not accept that (because it is too opaque).

You can also use "WKT" (that may include references to EPSG codes) but that is very verbose and not typically used in scripts to specify the crs. But it is how they are represented internally, and you can inspect them like this:

crs(r1) |> cat("\n")
#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["geodetic latitude (Lat)",north,
#            ORDER[1],
#            ANGLEUNIT["degree",0.0174532925199433]],
#        AXIS["geodetic longitude (Lon)",east,
#            ORDER[2],
#            ANGLEUNIT["degree",0.0174532925199433]],
#    USAGE[
#        SCOPE["Horizontal component of 3D system."],
#        AREA["World."],
#        BBOX[-90,-180,90,180]],
#    ID["EPSG",4326]] 

crs(r2) |> cat("\n")
#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]]]]

Again showing that r1 is a bit more fancy than r2.

Also not that you can use additional argument to rast, to set the desired extent, resolution, and number of layers. But with your cell.coords you can do:

cell.coords <- data.frame(lon = rep(c(2, 3), 2), lat =
   rep(c(1, 2), each = 2), cell.id = 10:13)
x <- rast(cell.coords, crs="+proj=longlat")
x
#class       : SpatRaster 
#dimensions  : 2, 2, 1  (nrow, ncol, nlyr)
#resolution  : 1, 1  (x, y)
#extent      : 1.5, 3.5, 0.5, 2.5  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : memory 
#name        : cell.id 
#min value   :      10 
#max value   :      13

I am not sure why you used rasterToPolygons in your example but if you wanted rectangular polygons instead or a raster you could do

p <- as.polygons(x)
#plot(p, "cell.id")

I would note that despite these pesky messages you get from rgdal, for most intents and purposes the PROJ.4 notation works fine, and makes for the most transparent (human-readable) code.

And you are right that the way forward is to forget about raster/sp and friends (rgdal, rgeos), and use terra and/or sf instead.

user_6352450
Related