如何迭代 SpatRaster 的像元值？

Question

我有 Python 背景，但最近一直在使用 R 的 terra 包，以便更好地处理地理空间栅格数据。

我一直想创建一个 for 循环，允许迭代 SpatRaster 的单元格值，以便填充另一个具有相同大小、分辨率和坐标的 SpatRaster。参考号取决于某些条件。

我的目标是创建一个类似于下面的循环；这适用于 Python 列表列表，我习惯了在 Python 中具有这种结构的数据框，但我无法在 R 中创建等效的东西：

for row in raster_1:
 for column in row:
  if value  value_1 <= x < value_1:
   raster_2[row, column] = 5

Answer 1

通常应尽可能避免循环遍历栅格单元，并使用子集化、栅格 algbera、

terra::ifel()

和其他

{terra}

方法 - https://rspatial.github.io/terra/reference/terra-package.html .

这不仅仅是 1 .. 5% 的性能增益，即使对于循环遍历每个单元格的微小栅格，应用一些逻辑并一次读/写单个单元格值也很容易比更合适的方法慢 1000 倍以上。

例如考虑这样的事情只移动特定值并将范围设置为固定值：

library(terra)
#> terra 1.7.39
par(mfrow = c(2, 2))

# input matrix 
m <- matrix(1:25, ncol = 5, byrow = TRUE)
m
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    1    2    3    4    5
#> [2,]    6    7    8    9   10
#> [3,]   11   12   13   14   15
#> [4,]   16   17   18   19   20
#> [5,]   21   22   23   24   25

# raster 1
rast_1 <- rast(m) |> `names<-`("raster_1")
rast_1
#> class       : SpatRaster 
#> dimensions  : 5, 5, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 0, 5, 0, 5  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> source(s)   : memory
#> name        : raster_1 
#> min value   :        1 
#> max value   :       25
plot(rast_1, main = "rast_1")

# raster 2, filled with NA values
rast_2 <- rast(rast_1, names = "raster_2", vals = NA)

# creating masks (TRUE / FALSE values):
(rast_1 %% 2 == 0) |> plot(main = "rast_1 %% 2 == 0")
(rast_1 > 10 & rast_1 <= 15) |> plot(main = "rast_1 > 10 & rast_1 <= 15")

## using those for submitting
# copy certain values from one raster to another:
rast_2[rast_1 %% 2 == 0] <- rast_1[rast_1 %% 2 == 0]

# set certain target values based on source values:
rast_2[rast_1 > 10 & rast_1 <= 15] <- 25
plot(rast_2, main = "rast_2[...] <- ...")

如果您确实需要循环单元格，这应该可以帮助您开始：

rast_3 <- rast(rast_1, names = "raster_3", vals = NA)
for (cid in 1:ncell(rast_1)){
  if(rast_1[cid] > 10 && rast_1[cid] <= 15){
    rast_3[cid] <- 5
  }
}

bench::mark()

子集 vs 具有简单条件的循环 vs

ifel()

虽然这是一个简单的基准测试，带有一个微小的 100x100 栅格，只是为了说明您应该做好准备，但所有三种方法的结果都是可比较的（

identical(rast_x[], rast_y[])

是

TRUE

）。

# 100x100 input rater for benchmark, 
# random values in range of 1..10, uniform distribution, 
# each value is used for ~10% of cells
set.seed(123)
rast_uniform <- 
  sample(1:10, size = 10000, replace = TRUE) |> 
  matrix(nrow = 100) |>
  rast()
rast_uniform
#> class       : SpatRaster 
#> dimensions  : 100, 100, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> source(s)   : memory
#> name        : lyr.1 
#> min value   :     1 
#> max value   :    10
summary(rast_uniform)
#>      lyr.1       
#>  Min.   : 1.000  
#>  1st Qu.: 3.000  
#>  Median : 6.000  
#>  Mean   : 5.507  
#>  3rd Qu.: 8.000  
#>  Max.   :10.000
table(rast_uniform[])
#> 
#>    1    2    3    4    5    6    7    8    9   10 
#> 1034 1003 1021  941  987  966  980 1028 1009 1031

# NA rasters for benchmark
rast_na_1 <- rast_na_2 <- rast_na_3 <- rast(rast_uniform, vals = NA)

# banchmark sets ~10% of cell values, subset vs loop vs ifel()
bnch <- bench::mark(
  subset = { rast_na_1[rast_uniform == 5] <- 0; rast_na_1 },
  loop   = { for (cid in 1:ncell(rast_uniform)) {if (rast_uniform[cid] == 5) rast_na_2[cid] <- 0}; rast_na_2 },
  ifel   = { rast_na_3 <- ifel(rast_uniform == 5, 0, NA); rast_na_3 },
  filter_gc = FALSE
)

注意单位，毫秒与秒，与取子集相比，测试循环慢了约 1700 倍（！）。

bnch
#> # A tibble: 3 × 6
#>   expression      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 subset       4.68ms   4.85ms  191.        6.26KB     5.97
#> 2 loop         10.26s   10.26s    0.0975    3.23MB     6.34
#> 3 ifel         7.33ms   8.71ms   94.9       8.84KB     5.93

^{创建于 2023-08-19，使用 reprex v2.0.2}