如何切换到r中的proj6

Question

我正在尝试从RGDAL和PROJ4切换出来，因为两者都被弃用。＆nbsp;我尝试使用SF分配一个坐标参考系统，但继续遇到有关CRS仍然是ProJ4的错误。 ＆nbsp;

    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)

＆nbsp; 我意识到我们被推荐了（ https://gis.stackexchange.com/questions/372692/how-should-ihold-i handle-crs-properly-after-after-the-major-change-change-change-in-proj-library ） 和/或 http://rgdal.r-forge.r-project.org/articles/proj6_gdal3.html.  这两个都是我不了解的长文档，并且涉及一个细节i级别i使用简单的射手不需要希望。

最佳解决方案将是一个单个命令，以创建空白。grid，以CRS作为参数。使用Terra而不是栅格是最好的，因为我认为我读到栅格也不是这个世界的时间。

Answer 1

您可以使用 rast 创建栅格模板（空白栅格）。有不同的方法来指定CRS（权威：代码，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 具有相同的坐标参考系统，可用于大多数目的和目的（经度/纬度，wgs84），但请注意上面打印的方式的差异。对于EPSG代码，通常可以使用诸如“ LON/LAT”之类的简短描述。

但是您尝试过的这些字符串：“ 4326+datum = wgs84”和“+init = epsg：4326+datum = wgs84”无效，因为它们混合了EPSG和ProJ。您可以使用单个EPSG代码，也可以通过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符号效果很好，并为最透明的（人类可读）代码提供了良好的作用。

您是对的，前进的道路是忘记栅格/SP和朋友（RGDAL，RGEO），并使用 Terra 和/或 sf 而不是。