R 语言，不合理的 (?) 光栅分辨率答案

【问题标题】：R language, unreasonable (?) raster resolutionsR 语言，不合理的 (?) 光栅分辨率
【发布时间】：2022-06-10 17:24:10
【问题描述】：

很抱歉这个非常愚蠢的问题，但我真的被困在这里......我需要为我的研究区域创建一个数字高程模型。为此，我下载了 SRTM（1 arc-seg 分辨率，可从网络免费获得）图像，其中包含比我感兴趣的区域更宽的区域。原始栅格具有以下特征：

class      : RasterLayer 
dimensions : 3601, 3601, 12967201  (nrow, ncol, ncell)
resolution : 0.0002777778, 0.0002777778  (x, y)
extent     : -45.00014, -43.99986, -22.00014, -20.99986  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : s22_w045_1arc_v3.tif 
names      : s22_w045_1arc_v3 
values     : -32768, 32767  (min, max)

我需要 (1) 将分辨率（最初为 30.75662 * 28.68392 m）提高到 1 * 1 m（也就是说，我真的不关心海拔的精确度）和 (2) 裁剪一个平方区域2000 * 2000 m 以给定坐标为中心。所以，我要做的第一步是重新投影到 UTM：

projection(r) <- "+proj=utm +zone=23 +datum=WGS84"

但之后分辨率单位不会改变：

class      : RasterLayer 
dimensions : 3601, 3601, 12967201  (nrow, ncol, ncell)
resolution : 0.0002777778, 0.0002777778  (x, y)
extent     : -45.00014, -43.99986, -22.00014, -20.99986  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=23 +datum=WGS84 +units=m +no_defs 
source     : s22_w045_1arc_v3.tif 
names      : s22_w045_1arc_v3 
values     : -32768, 32767  (min, max)

如果我尝试手动设置以米为单位的分辨率，则会生成一个空栅格。谁能在这里给我一些启发？

【问题讨论】：

标签： r raster resolution

【解决方案1】：

您正在更改（即覆盖）CRS，而不是投影栅格。通常，建议使用 CRS 和您需要的分辨率创建模板栅格，然后使用此模板重新投影栅格。

这里有一个例子，我正在切换到terra 进行分析，因为它是一个更新更快的包，但我也会展示如何将它转换回raster 格式：

library(raster)
#> Loading required package: sp

# Faking your data
r <- raster(
  nrows = 3601, ncols = 3601,
  ext = extent(c(-45.00014, -43.99986, -22.00014, -20.99986))
)



values(r) <- seq(-32768, 32767, length.out = ncell(r))
r
#> class      : RasterLayer 
#> dimensions : 3601, 3601, 12967201  (nrow, ncol, ncell)
#> resolution : 0.0002777784, 0.0002777784  (x, y)
#> extent     : -45.00014, -43.99986, -22.00014, -20.99986  (xmin, xmax, ymin, ymax)
#> crs        : +proj=longlat +datum=WGS84 +no_defs 
#> source     : memory
#> names      : layer 
#> values     : -32768, 32767  (min, max)
plot(r)

# End of faking data

# Change to terra, much faster
library(terra)
#> terra 1.5.21
r_terra <- terra::rast(r)
template <- terra::project(r_terra, "+proj=utm +zone=23 +datum=WGS84")

# Change to the desired res
res(template) <- c(2000, 2000)


# Reproject
r_terra_reproj <- terra::project(r_terra, template)

r_terra_reproj
#> class       : SpatRaster 
#> dimensions  : 56, 52, 1  (nrow, ncol, nlyr)
#> resolution  : 2000, 2000  (x, y)
#> extent      : 499985.4, 603985.4, -2433195, -2321195  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=utm +zone=23 +datum=WGS84 +units=m +no_defs 
#> source      : memory 
#> name        :     layer 
#> min value   : -32371.95 
#> max value   :  32182.81


terra::plot(r_terra_reproj)


# Back to RasterLayer

r_reproj <- raster(r_terra_reproj)

r_reproj
#> class      : RasterLayer 
#> dimensions : 56, 52, 2912  (nrow, ncol, ncell)
#> resolution : 2000, 2000  (x, y)
#> extent     : 499985.4, 603985.4, -2433195, -2321195  (xmin, xmax, ymin, ymax)
#> crs        : +proj=utm +zone=23 +datum=WGS84 +units=m +no_defs 
#> source     : memory
#> names      : layer 
#> values     : -32371.95, 32182.81  (min, max)

^{由reprex package (v2.0.1) 于 2022-06-10 创建}

【讨论】：