如何使用 terra 将光栅写入折旧的 crs？

Question

我正在处理折旧 crs 中的数据，并希望在将数据集转换为 SpatRaster 后避免重新投影。似乎 gdal 会自动将 EPSG:2163 替换为 EPSG:9311（请参阅 here）。 Gdal 明显显示警告，但 terra 没有显示警告。

我可以使用完整的 crs 字符串强制 SpatRaster 进入正确的 crs。后续操作（例如掩码）似乎根据需要使用折旧的 crs。

r <- rast(ncol = 1, nrow = 1, vals = 1,
     xmin = -1694992.5, xmax = -1694137.5, 
     ymin = -430492.5, ymax = -429367.5, 
     crs = 'EPSG:2163') # creates raster in EPSG:9311

crs(r) <- 'EPSG:2163' # doesn't work

# however, this successfully forces crs to EPSG:2163:
crs(r) <-"PROJCRS[\"US National Atlas Equal Area\",\n    BASEGEOGCRS[\"Unspecified datum based upon the Clarke 1866 Authalic Sphere\",\n        DATUM[\"Not specified (based on Clarke 1866 Authalic Sphere)\",\n            ELLIPSOID[\"Clarke 1866 Authalic Sphere\",6370997,0,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4052]],\n    CONVERSION[\"US National Atlas Equal Area\",\n        METHOD[\"Lambert Azimuthal Equal Area (Spherical)\",\n            ID[\"EPSG\",1027]],\n        PARAMETER[\"Latitude of natural origin\",45,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-100,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"False easting\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"easting (X)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"northing (Y)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Statistical analysis.\"],\n        AREA[\"United States (USA) - onshore and offshore.\"],\n        BBOX[15.56,167.65,74.71,-65.69]],\n    ID[\"EPSG\",2163]]"

但是，当我将栅格写入 GeoTIFF 时，我还没有想出如何覆盖 gdal 来替换 crs。

tf <- tempfile(fileext = '.tif')
writeRaster(r, tf) # this doesn't work

writeRaster(r, tf, gdal = 'OSR_USE_NON_DEPRECATED=NO') # neither does this

笨拙的解决方案是在将 GeoTiff 加载到 R 后仅使用完整的 crs 字符串重新定义 crs。是否有更好的方法来覆盖写入和读取文件的 crs 替换？

Answer 1

你的例子数据

library(terra)
r <- rast(ncol = 1, nrow = 1, vals = 1,
     xmin = -1694992.5, xmax = -1694137.5, 
     ymin = -430492.5, ymax = -429367.5, 
     crs = 'EPSG:2163') # creates raster in EPSG:9311

crs(r) <-"PROJCRS[\"US National Atlas Equal Area\", BASEGEOGCRS[\"Unspecified datum based upon the Clarke 1866 Authalic Sphere\", DATUM[\"Not specified (based on Clarke 1866 Authalic Sphere)\", ELLIPSOID[\"Clarke 1866 Authalic Sphere\",6370997,0, LENGTHUNIT[\"metre\",1]]], PRIMEM[\"Greenwich\",0, ANGLEUNIT[\"degree\",0.0174532925199433]], ID[\"EPSG\",4052]], CONVERSION[\"US National Atlas Equal Area\", METHOD[\"Lambert Azimuthal Equal Area (Spherical)\", ID[\"EPSG\",1027]], PARAMETER[\"Latitude of natural origin\",45, ANGLEUNIT[\"degree\",0.0174532925199433], ID[\"EPSG\",8801]], PARAMETER[\"Longitude of natural origin\",-100, ANGLEUNIT[\"degree\",0.0174532925199433], ID[\"EPSG\",8802]], PARAMETER[\"False easting\",0, LENGTHUNIT[\"metre\",1], ID[\"EPSG\",8806]], PARAMETER[\"False northing\",0, LENGTHUNIT[\"metre\",1], ID[\"EPSG\",8807]]], CS[Cartesian,2], AXIS[\"easting (X)\",east, ORDER[1], LENGTHUNIT[\"metre\",1]], AXIS[\"northing (Y)\",north, ORDER[2], LENGTHUNIT[\"metre\",1]], USAGE[ SCOPE[\"Statistical analysis.\"], AREA[\"United States (USA) - onshore and offshore.\"], BBOX[15.56,167.65,74.71,-65.69]], ID[\"EPSG\",2163]]"

默认值：

tf <- tempfile(fileext = '.tif')
(writeRaster(r, tf, overwrite=TRUE))
# coord. ref. : NAD27 / US National Atlas Equal Area (EPSG:9311)

解决方案 1，使用 PROJ.4 表示法

## get the proj.4 value 
prj <- crs(r, proj=TRUE)
prj 
#[1] "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +ellps=sphere +units=m +no_defs"
crs(r) <- prj
(writeRaster(r, tf, overwrite=TRUE))
#coord. ref. : +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +ellps=sphere +units=m +no_defs

方案二，使用setGDALconfig。这需要 terra 版本 >= 1.5-27（目前是开发版本，在 windows 上你可以用 install.packages('terra', repos='https://rspatial.r-universe.dev') 安装它）

setGDALconfig("OSR_USE_NON_DEPRECATED", value="NO")
(writeRaster(r, tf, overwrite=TRUE))
# coord. ref. : US National Atlas Equal Area (EPSG:2163)

如何使用 terra 将光栅写入折旧的 crs？

How to write raster to depreciated crs with terra?

r

terra