子集栅格以分隔两个类

Question

我有一个太大的光栅，无法在此处发布，但可以从此处下载：https://login.filesanywhere.com/fs/v.aspx?v=8c6c688b586472bcab6a。当我绘制它时，我看到水和土地的混合体。我想用以下标准将其分开：如果值大于 1，则为土地，如果值小于 1，则为水。

然后，将水类转换为多边形特征。然后用它来裁剪原始光栅以仅显示水。非常感谢任何帮助。

Here is what I tried:

  library(terra)
  library(mapview)
  library(raster)
  
  xx <- rast("MyRaster.asc")
xx
ncell(xx)
plot(xx)

#Convert to raster first to vizualise with mapview
xxx <- raster(xx)
mapview(xxx)

Answer 1

您可以使用 terra 轻松做到这一点，避免转换和潜在的错误。请参阅https://gis.stackexchange.com/questions/421821/how-to-subset-a-spatraster-by-value-in-rterra开发人员对此进行了进一步解释。

模拟文件示例：

library(terra)
#> terra 1.7.23

# Mock data
f <- system.file("extdata/asia.tif", package = "tidyterra")
xx <- rast(f)
# End of Mock data

# xx <- rast("MyRaster.asc")

xx
#> class       : SpatRaster 
#> dimensions  : 232, 432, 1  (nrow, ncol, nlyr)
#> resolution  : 22550.66, 22512.94  (x, y)
#> extent      : 7619120, 17361007, -1304745, 3918256  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
#> source      : asia.tif 
#> name        : file44bc291153f2 
#> min value   :        -10071.50 
#> max value   :          6064.73
ncell(xx)
#> [1] 100224
plot(xx)

# Clamp it

# Water
xx_water <- clamp(xx, upper = 1, value = FALSE)
xx_water
#> class       : SpatRaster 
#> dimensions  : 232, 432, 1  (nrow, ncol, nlyr)
#> resolution  : 22550.66, 22512.94  (x, y)
#> extent      : 7619120, 17361007, -1304745, 3918256  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
#> source(s)   : memory
#> name        : file44bc291153f2 
#> min value   :    -1.007150e+04 
#> max value   :     9.916311e-01

plot(xx_water, col = hcl.colors(10, "Blues2"))

# Land
xx_land <- clamp(xx, lower = 1, value = FALSE)
xx_land
#> class       : SpatRaster 
#> dimensions  : 232, 432, 1  (nrow, ncol, nlyr)
#> resolution  : 22550.66, 22512.94  (x, y)
#> extent      : 7619120, 17361007, -1304745, 3918256  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
#> source(s)   : memory
#> name        : file44bc291153f2 
#> min value   :         1.009122 
#> max value   :      6064.730469

plot(xx_land, col = hcl.colors(10, "BuGn"))

^{创建于 2023-04-21 与 reprex v2.0.2}

我试过你的文件（~1.3Gb），这个过程可能需要一些时间，但上面的代码应该仍然有效。仅查看水的结果：

library(terra)
#> terra 1.7.23

xx <- rast("del/sfbaydeltadem10m2016.asc")

xx

# class       : SpatRaster 
# dimensions  : 15795, 13588, 1  (nrow, ncol, nlyr)
# resolution  : 10, 10  (x, y)
# extent      : 520545, 656425, 4136705, 4294655  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=utm +zone=10 +datum=NAD83 +units=m +no_defs 
# source      : sfbaydeltadem10m2016.asc 
# name        : sfbaydeltadem10m2016 
# min value   :            -114.8981 
# max value   :             432.1640 
ncell(xx)
#> [1] 214622460


# Clamp it

# Water
xx_water <- clamp(xx, upper = 1, value = FALSE)
plot(xx_water, col = hcl.colors(10, "Blues2"))

Answer 2

读取数据：

library(terra)
xx = rast("sfbaydeltadem10m2016.asc")

创建一个二进制 0/1 栅格，其中陆地为 1，海洋为 0：

landsea = xx > 1

转换为多边形（需要一段时间）：

pol = as.polygons(landsea)

显示（也需要一段时间）：

plot(pol, col=c("red","blue")

注意这里的两个“多边形”是非常复杂的多边形，每行

pol

代表该值的所有小区域。

要完成将栅格裁剪到水域范围的最后一步，您实际上不必经历这个多边形化过程，只需将所有土地值设置为缺失值并从各个方向切割完整的行和列 - 可能有

terra

中的一个函数来执行此操作 - 是的，它称为

trim()

，您可以向它传递一个值以计为缺失。