如何使用terra将栅格数据写入已弃用的CRS？

Question

我正在处理一个使用废弃的crs格式的数据集，并希望在将数据集转换为SpatRaster后避免重新投影。似乎gdal会自动将EPSG:2163替换为EPSG:9311（请参见此处）。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替换，以便在读写文件时使用？

- Joseph

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) 

- Robert Hijmans