R 栅格函数将不接受 WKT 格式的 crs

Posted 2023-03-23

【中文标题】R 栅格函数将不接受 WKT 格式的 crs

【英文标题】：R raster function will not accept crs in WKT format

【发布时间】：2021-08-07 07:01:14

【问题描述】：

我正在尝试生成栅格并为其分配 CRS 投影。但是，我的 CRS 采用新的 WKT 格式，并且 raster() 函数要求我提供 proj4string。到目前为止，这是我的代码：

library(sf)
library(raster)

crs_dem <- st_crs(
  'PROJCS["NAD_1983_2011_StatePlane_California_II_FIPS_0402",
    GEOGCS["GCS_NAD_1983_2011",
      DATUM["D_NAD_1983_2011",
         SPHEROID["GRS_1980",6378137.0,298.257222101]],
       PRIMEM["Greenwich",0.0],
       UNIT["Degree",0.0174532925199433]],
     PROJECTION["Lambert_Conformal_Conic"],
     PARAMETER["False_Easting",2000000.0],
     PARAMETER["False_Northing",500000.0],
     PARAMETER["Central_Meridian",-122.0],
     PARAMETER["Standard_Parallel_1",38.33333333333334],
     PARAMETER["Standard_Parallel_2",39.83333333333334],
     PARAMETER["Latitude_Of_Origin",37.66666666666666],
     UNIT["Meter",1.0]]')

ext <- extent(1895000, 1935000, 579500, 616500)
grid <- raster(ext, resolution = c(40,40), crs = crs(dem))

上面的代码生成了一个 crs = NA 的栅格。我也试过用 crs(grid)

【问题讨论】：

您的对象属于crs 类（sf 包），而raster 需要CRS 类（例如sp 包）。他们是两个不同的东西。例如，尝试test <- sp::CRS(crs_dem$input) 并提供它；它会起作用的。但是，我对这些软件包没有经验，所以我不知道这是否是一个有效的解决方案。这只是说明出了什么问题。

【参考方案1】：

@slamballais 的回答成功了！昨晚我还通过反复试验找到了另一种（不太干净）的方法，所以这里有两种解决方案。

选项 1：

test <- sp::CRS(crs_dem$input)
grid <- raster(ext, resolution = c(40,40), crs = test)

选项 2：

library(dplyr)
aoi <- ext %>% 
  as('SpatialPolygons') %>% 
  st_as_sf %>% 
  st_set_crs(crs_dem)
grid <- raster(ext, resolution = c(40,40), crs = projection(aoi))

【参考方案2】：

raster 需要文本，所以如果你有 wkt 格式 crs，你可以直接使用它。无需创建更复杂的对象。

crs_dem <- 'PROJCS["NAD_1983_2011_StatePlane_California_II_FIPS_0402",
    GEOGCS["GCS_NAD_1983_2011", DATUM["D_NAD_1983_2011", SPHEROID["GRS_1980",6378137.0,298.257222101]],
     PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]],
     PROJECTION["Lambert_Conformal_Conic"], PARAMETER["False_Easting",2000000.0],
     PARAMETER["False_Northing",500000.0], PARAMETER["Central_Meridian",-122.0],
     PARAMETER["Standard_Parallel_1",38.33333333333334], PARAMETER["Standard_Parallel_2",39.83333333333334],
     PARAMETER["Latitude_Of_Origin",37.66666666666666], UNIT["Meter",1.0]]'
     
library(raster)
r <- raster(crs=crs_dem)

或者，如果您以 sf 对象开头

r <- raster(crs=crs_dem$input)