How to read output of crs() in the Terra package in R

Question

I have a SpatRaster and I want to figure out the EPSG code of it. I have used the crs() command and I got the following output:

PROJCRS["unknown",
        BASEGEOGCRS["NAD83",
                    DATUM["North American Datum 1983",
                          ELLIPSOID["GRS 1980",6378137,298.257222101004,
                                    LENGTHUNIT["metre",1]]],
                    PRIMEM["Greenwich",0,
                           ANGLEUNIT["degree",0.0174532925199433]],
                    ID["EPSG",4269]],
        CONVERSION["Albers Equal Area",
                   METHOD["Albers Equal Area",
                          ID["EPSG",9822]],
                   PARAMETER["Latitude of false origin",23,
                             ANGLEUNIT["degree",0.0174532925199433],
                             ID["EPSG",8821]],
                   PARAMETER["Longitude of false origin",-96,
                             ANGLEUNIT["degree",0.0174532925199433],
                             ID["EPSG",8822]],
                   PARAMETER["Latitude of 1st standard parallel",29.5,
                             ANGLEUNIT["degree",0.0174532925199433],
                             ID["EPSG",8823]],
                   PARAMETER["Latitude of 2nd standard parallel",45.5,
                             ANGLEUNIT["degree",0.0174532925199433],
                             ID["EPSG",8824]],
                   PARAMETER["Easting at false origin",0,
                             LENGTHUNIT["metre",1],
                             ID["EPSG",8826]],
                   PARAMETER["Northing at false origin",0,
                             LENGTHUNIT["metre",1],
                             ID["EPSG",8827]]],
        CS[Cartesian,2],
        AXIS["easting",east,
             ORDER[1],
             LENGTHUNIT["metre",1,
                        ID["EPSG",9001]]],
        AXIS["northing",north,
             ORDER[2],
             LENGTHUNIT["metre",1,
                        ID["EPSG",9001]]]]

As you can see, there are multiple different numbers after the word EPSG and I am unsure how to read it and figure out what actually is the EPSG of my data. Can someone help me with that? Below is a reproducible example:

Get the SpatRaster data. To download it on Linux using the R package CropscapeR I need to do the process below (in Windows it is easier):

library(CropScapeR)
library(terra)
# Skip the SSL check
httr::set_config(httr::config(ssl_verifypeer = 0L))
# Automatically generate a temporary path to save the data
tif_file <- tempfile(fileext = '.tif')
# Download the raster TIF file into specified path, also read into R
# This is a county in Iowa, US
ST_CDL <- GetCDLData(aoi = '19001', year = 2020, type = 'f', save_path = tif_file) 
terra::writeRaster(ST_CDL, "ST_CDL.tif", overwrite=TRUE)
ST_CDL = terra::rast("ST_CDL.tif")

The command I run that gives me the output shown above.

crs(ST_CDL)

I want to figure out what the EPSG code is in order to terra::project() a SpatVector I have to the crs of my SpatRaster above. However, I was afraid terra::project(my_spatvector, crs(ST_CDL)) would give me the wrong output if crs(ST_CDL) is PROJ.4, since the documentation for this functions says "note that the PROJ.4 notation has been deprecated, and you can only use it with the WGS84/NAD83 and NAD27 datums. Other datums are silently ignored."

Answer 1

There's a good chance that the used CRS does not perfectly match any of the registered EPSG codes. Though I don't see any real issues here, by default terra::crs() returns CRS as a WKT string (as you demonstrate in your question), to get PROJ-string you'd have to force it by using terra::crs(ST_CDL, proj = TRUE) .

CropScapeR::GetCDLData() in your example saves TIF, returns RasterLayer, which you then convert to SpatRaster and save to another TIF. So perhaps just skip that raster to terra conversion and open the downloaded file directly with terra:

tif_file <- "ST_CDL.tif"
# set readr = FALSE to skip reading downloaded file with  
# raster::raster(save_path)
CropScapeR::GetCDLData(aoi = '19001', year = 2020, type = 'f', save_path = tif_file, readr = FALSE) 
#> Data is saved at:ST_CDL.tif
#> NULL

# read with terra
ST_CDL <- terra::rast(tif_file)
ST_CDL
#> class       : SpatRaster 
#> dimensions  : 1312, 1292, 1  (nrow, ncol, nlyr)
#> resolution  : 30, 30  (x, y)
#> extent      : 107655, 146415, 2017395, 2056755  (xmin, xmax, ymin, ymax)
#> coord. ref. : NAD83_Conus_Albers 
#> source      : ST_CDL.tif 
#> color table : 1 
#> name        : Layer_1

# using terra::crs() result for vector crs is fine, 
# crs() returns PROJ-string only if you explicitly set proj to TRUE:
terra::crs(ST_CDL, proj = TRUE)
#> [1] "+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"

^{Created on 2023-08-15 with reprex v2.0.2}