如何读取 R 中 Terra 包中 crs() 的输出

Question

我有一个 SpatRaster，我想找出它的 EPSG 代码。我使用了

crs()

命令，得到了以下输出：

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]]]]

如您所见，EPSG 一词后面有多个不同的数字，我不确定如何读取它并弄清楚我的数据的 EPSG 到底是什么。有人可以帮我吗？以下是一个可重现的示例：

获取 SpatRaster 数据。要使用 R 包在 Linux 上下载它

CropscapeR

，我需要执行以下过程（在 Windows 中更容易）：

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")

我运行的命令给出了上面显示的输出。

crs(ST_CDL)

我想弄清楚 EPSG 代码是什么，以便

terra::project()

获得上面 SpatRaster 的 crs 的 SpatVector。但是，我担心如果 crs(ST_CDL) 是 PROJ.4，

terra::project(my_spatvector, crs(ST_CDL))

会给我错误的输出，因为此函数的 documentation 说“请注意，PROJ.4 符号已被弃用，您只能使用它与 WGS84/NAD83 和 NAD27 基准相匹配。其他基准将被默默忽略。”

Answer 1

使用的 CRS 很可能与任何注册的 EPSG 代码不完全匹配。虽然我在这里没有看到任何真正的问题，但默认情况下

terra::crs()

将 CRS 作为 WKT 字符串返回（正如您在问题中演示的那样），要获取 PROJ 字符串，您必须使用

terra::crs(ST_CDL, proj = TRUE)

来强制它。

示例中的

CropScapeR::GetCDLData()

保存 TIF，返回

RasterLayer

，然后将其转换为

SpatRaster

并保存到另一个 TIF。因此，也许只需跳过

raster

到

terra

的转换，直接使用

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"

^{创建于 2023-08-15，使用 reprex v2.0.2}