如何使用`terra`通过R中的条件语句对栅格进行子集化？答案

【问题标题】：How can I subset a raster by conditional statement in R using `terra`?如何使用`terra`通过R中的条件语句对栅格进行子集化？
【发布时间】：2022-01-18 08:39:49
【问题描述】：

我正在尝试仅从我正在使用的分类土地覆盖栅格中绘制某些值。我已经使用 terra 包将它加载到 R 中，并且绘制得很好。但是，由于原始数据没有附带图例，我试图找出哪个栅格值对应于地图上的值。

与此处提供的答案类似：How to subset a raster based on grid cell values

我已尝试使用以下行：

> landcover
class       : SpatRaster 
dimensions  : 20057, 63988, 1  (nrow, ncol, nlyr)
resolution  : 0.0005253954, 0.0005253954  (x, y)
extent      : -135.619, -102, 59.99989, 70.53775  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : spat_n5WpgzBuVAV3Ijm.tif 
name        : CAN_LC_2015_CAL_wgs 
min value   :                   1 
max value   :                  18 

> plot(landcover[landcover == 18])
                                      
Error: cannot allocate vector of size 9.6 Gb

但是，此行需要很长时间才能运行并产生向量内存错误。该对象在全局环境中为 1.3 kb，原始 tif 约为 300 mb。

【问题讨论】：

标签： r raster terra

【解决方案1】：

您可以使用cats 找出哪些值对应于哪些类别。

library(terra)
set.seed(0)
r <- rast(nrows=10, ncols=10)
values(r) <- sample(3, ncell(r), replace=TRUE) - 1
cls <- c("forest", "water", "urban")
levels(r) <- cls
names(r) <- "land cover"

cats(r)[[1]]
#  ID category
#1  0   forest
#2  1    water
#3  2    urban

要为一个类别绘制逻辑（布尔）层，您可以这样做

plot(r == "water")

从上面你可以看出，在这种情况下，它相当于

plot(r == 1)

【讨论】：

【解决方案2】：

我想我找到了在绘图函数中编写条件的解决方案，如下所示：

plot(landcover == 18)

对于那些寻找可复制示例的人，只需加载 rlogo：

s <- rast(system.file("ex/logo.tif", package="terra"))
s <- s$red
plot(s == 255)

【讨论】：