如何在 r 中切换到 PROJ6
How to switch to PROJ6 in r
我正在尝试放弃 rgdal 和 PROJ4，因为它们都已弃用。我尝试使用 sf 分配一个坐标参考系统，但继续收到关于 crs 仍然是 PROJ4 的错误。
pacman::p_load(sf, raster)
cell.coords <- data.frame(lon = rep(c(2, 3), 2), lat =
rep(c(1, 2), each = 2), cell.id = 10:13)
blank.grid <- st_as_sf(rasterToPolygons(rasterFromXYZ(
cell.coords[ , c("lon", "lat", "cell.id")])))
# Try, in turn, each of the following.
# First three successfully show the crs when blank.grid is called,
# but crs(blank.grid) includes "Deprecated Proj.4 representation"
wgs84 <- 4326
# next is based on https://gis.stackexchange.com/questions/372692/how-should-i-handle-crs-properly-after-the-major-change-in-proj-library
wgs84 <- "OGC:CRS84"
# next is based on https://gis.stackexchange.com/questions/385114/how-to-update-crs-r-object-from-proj4-to-proj6
wgs84 <- "+proj=lcc +lat_1=49 +lat_2=77 +lat_0=0 +lon_0=-95 +x_0=0 +y_0=0 +ellps=GRS80 +datum=WGS84 +units=m +no_defs"
# For the rest, blank.grid returns CRS: NA
wgs84 <- "4326+datum=WGS84"
wgs84 <- "+init=epsg:4326+datum=WGS84"
wgs84 <- "EPSG:4326"
# Apply and check each option with
st_crs(blank.grid) <- wgs84
blank.grid
crs(blank.grid)
我意识到我们被引用了（https://gis.stackexchange.com/questions/372692/how-should-i-handle-crs-properly-after-the-major-change-in-proj-library) to the technical references https://cran.r-project.org/web/packages/rgdal/vignettes/CRS_projections_transformations.html and/or http://rgdal.r-forge.r-project.org/articles/PROJ6_GDAL3.html. 两个都是我不理解的长文档，并且涉及到一定程度的细节，我希望使用简单的栅格不需要。
最佳解决方案是创建 blank.grid 的单个命令，并将 crs 作为参数。使用 terra 而不是栅格是最好的，因为我想我读到栅格对于这个世界来说同样不长。
您可以使用 rast 创建光栅模板（空白光栅）。有不同的方式来指定 crs（authority:code、PROJ 或 WKT）。例如：
library(terra)
# authority:code
r1 <- rast(crs="EPSG:4326")
r1
#class : SpatRaster
#dimensions : 180, 360, 1 (nrow, ncol, nlyr)
#resolution : 1, 1 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (EPSG:4326)
# PROJ
r2 <- rast(crs="+proj=longlat")
r2
#class : SpatRaster
#dimensions : 180, 360, 1 (nrow, ncol, nlyr)
#resolution : 1, 1 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs
r1 和 r2 对于大多数意图和目的（longitude/latitude，WGS84）具有相同的坐标参考系统，但请注意它们在上面的打印方式上的差异。对于 EPSG 代码，通常可以使用像“lon/lat”这样的简短描述。
但是您尝试的这些字符串：“4326+datum=WGS84”和“+init=epsg:4326+datum=WGS84”无效，因为它们混合了 EPSG 和 PROJ。您可以使用单个 EPSG 代码，或通过 PROJ 表示法指定参数。一个例外是像这样使用 PROJ 表示法：“+init=epsg:4326”（但没有理由再使用该表示法，因为您可以改用“epsg:4326”）。
用sf也可以把EPSG码指定为数字“4326”，不带权限前缀，但是terra不接受（因为太不透明）
您也可以使用“WKT”（可能包括对 EPSG 代码的引用），但这非常冗长，通常不在脚本中用于指定 crs。但这是它们在内部的表示方式，您可以这样检查它们：
crs(r1) |> cat("\n")
#GEOGCRS["WGS 84",
# DATUM["World Geodetic System 1984",
# ELLIPSOID["WGS 84",6378137,298.257223563,
# LENGTHUNIT["metre",1]]],
# PRIMEM["Greenwich",0,
# ANGLEUNIT["degree",0.0174532925199433]],
# CS[ellipsoidal,2],
# AXIS["geodetic latitude (Lat)",north,
# ORDER[1],
# ANGLEUNIT["degree",0.0174532925199433]],
# AXIS["geodetic longitude (Lon)",east,
# ORDER[2],
# ANGLEUNIT["degree",0.0174532925199433]],
# USAGE[
# SCOPE["Horizontal component of 3D system."],
# AREA["World."],
# BBOX[-90,-180,90,180]],
# ID["EPSG",4326]]
crs(r2) |> cat("\n")
#GEOGCRS["unknown",
# DATUM["World Geodetic System 1984",
# ELLIPSOID["WGS 84",6378137,298.257223563,
# LENGTHUNIT["metre",1]],
# ID["EPSG",6326]],
# PRIMEM["Greenwich",0,
# ANGLEUNIT["degree",0.0174532925199433],
# ID["EPSG",8901]],
# CS[ellipsoidal,2],
# AXIS["longitude",east,
# ORDER[1],
# ANGLEUNIT["degree",0.0174532925199433,
# ID["EPSG",9122]]],
# AXIS["latitude",north,
# ORDER[2],
# ANGLEUNIT["degree",0.0174532925199433,
# ID["EPSG",9122]]]]
再次表明 r1 比 r2 更花哨。
也不是说您可以使用 rast 的附加参数来设置所需的范围、分辨率和层数。但是用你的cell.coords你可以做到：
cell.coords <- data.frame(lon = rep(c(2, 3), 2), lat =
rep(c(1, 2), each = 2), cell.id = 10:13)
x <- rast(cell.coords, crs="+proj=longlat")
x
#class : SpatRaster
#dimensions : 2, 2, 1 (nrow, ncol, nlyr)
#resolution : 1, 1 (x, y)
#extent : 1.5, 3.5, 0.5, 2.5 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs
#source : memory
#name : cell.id
#min value : 10
#max value : 13
我不确定你为什么在你的例子中使用 rasterToPolygons 但如果你想要矩形多边形或栅格你可以做
p <- as.polygons(x)
#plot(p, "cell.id")
我会注意到，尽管您从 rgdal 收到了这些讨厌的消息，但对于大多数意图和目的而言，PROJ.4 表示法工作正常，并且是最透明的 (human-readable) 代码。
你是对的，前进的道路是忘记 raster/sp 和朋友 (rgdal, rgeos)，而是使用 terra and/or sf。
我正在尝试放弃 rgdal 和 PROJ4，因为它们都已弃用。我尝试使用 sf 分配一个坐标参考系统，但继续收到关于 crs 仍然是 PROJ4 的错误。
pacman::p_load(sf, raster)
cell.coords <- data.frame(lon = rep(c(2, 3), 2), lat =
rep(c(1, 2), each = 2), cell.id = 10:13)
blank.grid <- st_as_sf(rasterToPolygons(rasterFromXYZ(
cell.coords[ , c("lon", "lat", "cell.id")])))
# Try, in turn, each of the following.
# First three successfully show the crs when blank.grid is called,
# but crs(blank.grid) includes "Deprecated Proj.4 representation"
wgs84 <- 4326
# next is based on https://gis.stackexchange.com/questions/372692/how-should-i-handle-crs-properly-after-the-major-change-in-proj-library
wgs84 <- "OGC:CRS84"
# next is based on https://gis.stackexchange.com/questions/385114/how-to-update-crs-r-object-from-proj4-to-proj6
wgs84 <- "+proj=lcc +lat_1=49 +lat_2=77 +lat_0=0 +lon_0=-95 +x_0=0 +y_0=0 +ellps=GRS80 +datum=WGS84 +units=m +no_defs"
# For the rest, blank.grid returns CRS: NA
wgs84 <- "4326+datum=WGS84"
wgs84 <- "+init=epsg:4326+datum=WGS84"
wgs84 <- "EPSG:4326"
# Apply and check each option with
st_crs(blank.grid) <- wgs84
blank.grid
crs(blank.grid)
我意识到我们被引用了（https://gis.stackexchange.com/questions/372692/how-should-i-handle-crs-properly-after-the-major-change-in-proj-library) to the technical references https://cran.r-project.org/web/packages/rgdal/vignettes/CRS_projections_transformations.html and/or http://rgdal.r-forge.r-project.org/articles/PROJ6_GDAL3.html. 两个都是我不理解的长文档，并且涉及到一定程度的细节，我希望使用简单的栅格不需要。
最佳解决方案是创建 blank.grid 的单个命令，并将 crs 作为参数。使用 terra 而不是栅格是最好的，因为我想我读到栅格对于这个世界来说同样不长。
您可以使用 rast 创建光栅模板（空白光栅）。有不同的方式来指定 crs（authority:code、PROJ 或 WKT）。例如：
library(terra)
# authority:code
r1 <- rast(crs="EPSG:4326")
r1
#class : SpatRaster
#dimensions : 180, 360, 1 (nrow, ncol, nlyr)
#resolution : 1, 1 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (EPSG:4326)
# PROJ
r2 <- rast(crs="+proj=longlat")
r2
#class : SpatRaster
#dimensions : 180, 360, 1 (nrow, ncol, nlyr)
#resolution : 1, 1 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs
r1 和 r2 对于大多数意图和目的（longitude/latitude，WGS84）具有相同的坐标参考系统，但请注意它们在上面的打印方式上的差异。对于 EPSG 代码，通常可以使用像“lon/lat”这样的简短描述。
但是您尝试的这些字符串：“4326+datum=WGS84”和“+init=epsg:4326+datum=WGS84”无效，因为它们混合了 EPSG 和 PROJ。您可以使用单个 EPSG 代码，或通过 PROJ 表示法指定参数。一个例外是像这样使用 PROJ 表示法：“+init=epsg:4326”（但没有理由再使用该表示法，因为您可以改用“epsg:4326”）。
用sf也可以把EPSG码指定为数字“4326”，不带权限前缀，但是terra不接受（因为太不透明）
您也可以使用“WKT”（可能包括对 EPSG 代码的引用），但这非常冗长，通常不在脚本中用于指定 crs。但这是它们在内部的表示方式，您可以这样检查它们：
crs(r1) |> cat("\n")
#GEOGCRS["WGS 84",
# DATUM["World Geodetic System 1984",
# ELLIPSOID["WGS 84",6378137,298.257223563,
# LENGTHUNIT["metre",1]]],
# PRIMEM["Greenwich",0,
# ANGLEUNIT["degree",0.0174532925199433]],
# CS[ellipsoidal,2],
# AXIS["geodetic latitude (Lat)",north,
# ORDER[1],
# ANGLEUNIT["degree",0.0174532925199433]],
# AXIS["geodetic longitude (Lon)",east,
# ORDER[2],
# ANGLEUNIT["degree",0.0174532925199433]],
# USAGE[
# SCOPE["Horizontal component of 3D system."],
# AREA["World."],
# BBOX[-90,-180,90,180]],
# ID["EPSG",4326]]
crs(r2) |> cat("\n")
#GEOGCRS["unknown",
# DATUM["World Geodetic System 1984",
# ELLIPSOID["WGS 84",6378137,298.257223563,
# LENGTHUNIT["metre",1]],
# ID["EPSG",6326]],
# PRIMEM["Greenwich",0,
# ANGLEUNIT["degree",0.0174532925199433],
# ID["EPSG",8901]],
# CS[ellipsoidal,2],
# AXIS["longitude",east,
# ORDER[1],
# ANGLEUNIT["degree",0.0174532925199433,
# ID["EPSG",9122]]],
# AXIS["latitude",north,
# ORDER[2],
# ANGLEUNIT["degree",0.0174532925199433,
# ID["EPSG",9122]]]]
再次表明 r1 比 r2 更花哨。
也不是说您可以使用 rast 的附加参数来设置所需的范围、分辨率和层数。但是用你的cell.coords你可以做到：
cell.coords <- data.frame(lon = rep(c(2, 3), 2), lat =
rep(c(1, 2), each = 2), cell.id = 10:13)
x <- rast(cell.coords, crs="+proj=longlat")
x
#class : SpatRaster
#dimensions : 2, 2, 1 (nrow, ncol, nlyr)
#resolution : 1, 1 (x, y)
#extent : 1.5, 3.5, 0.5, 2.5 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs
#source : memory
#name : cell.id
#min value : 10
#max value : 13
我不确定你为什么在你的例子中使用 rasterToPolygons 但如果你想要矩形多边形或栅格你可以做
p <- as.polygons(x)
#plot(p, "cell.id")
我会注意到，尽管您从 rgdal 收到了这些讨厌的消息，但对于大多数意图和目的而言，PROJ.4 表示法工作正常，并且是最透明的 (human-readable) 代码。
你是对的，前进的道路是忘记 raster/sp 和朋友 (rgdal, rgeos)，而是使用 terra and/or sf。