识别 R 栅格包中的重叠区域

Posted 2023-03-23

【中文标题】识别 R 栅格包中的重叠区域

【英文标题】：Identifying overlap zones in R raster package

【发布时间】：2011-08-05 19:57:27

【问题描述】：

包装：

raster

数据：

具有 10 个波段的 rasterStack。 每个波段都包含一个被 NA 包围的图像区域 波段是合乎逻辑的，即图像数据为“1”，周围区域为“0”/NA 每个波段的“图像区域”彼此不完全对齐，尽管大多数有部分重叠

目标：

编写一个快速函数，该函数可以返回每个“区域”的 rasterLayer 或像元编号，例如，仅包含波段 1 和 2 的数据的像素位于区域 1，仅包含波段 3 和 4 的数据的像素位于在区域 2 等。如果返回 rasterLayer，我需要能够稍后将区域值与波段编号匹配。

第一次尝试：

# Possible band combinations
values = integer(0)
for(i in 1:nlayers(myraster))
 combs = combn(1:nlayers(myraster), i)
 for(j in 1:ncol(combs))
  values = c(values, list(combs[,j]))
 


# Define the zone finding function
find_zones = function(bands)

 # The intersection of the bands of interest
 a = subset(myraster, 1)
 values(a) = TRUE
 for(i in bands)
  a = a & myraster[[i]]
 

 # Union of the remaining bands
 b = subset(myraster, 1)
 values(b) = FALSE
 for(i in seq(1:nlayers(myraster))[-bands])
  b = b | myraster[[i]]
 

 #plot(a & !b)
 cells = Which(a & !b, cells=TRUE)
 return(cells)


# Applying the function
results = lapply(values, find_zones)

我当前的函数需要很长时间才能执行。你能想出更好的方法吗？请注意，我不只是想知道每个像素有多少波段有数据，我还需要知道哪些波段。这样做的目的是在之后以不同的方式处理不同的区域。

另请注意，实际场景是 3000 x 3000 或更大的栅格，可能有超过 10 个波段。

---

编辑

由10个偏移图像区域组成的一些样本数据：

# Sample data
library(raster)    
for(i in 1:10) 
  start_line = i*10*1000
  end_line = 1000000 - 800*1000 - start_line
  offset = i * 10
  data = c(rep(0,start_line), rep(c(rep(0,offset), rep(1,800), rep(0,200-offset)), 800), rep(0, end_line))
  current_layer = raster(nrows=1000, ncols=1000)
  values(current_layer) = data
  if(i == 1) 
    myraster = stack(current_layer)
   else 
    myraster = addLayer(myraster, current_layer)
  

NAvalue(myraster) = 0  # You may not want to do this depending on your solution...

【问题讨论】：

你能详细说明什么是“区域”吗？

我会定义一个“区域”，一组单元格具有相同波段的数据（并且只有那些共同的波段）。例如，如果您有两个图层，每个图层都是正方形，但有一个偏移 100 像素，那么您将有 3 个区域，一个只有波段 1，一个只有波段 2，一个有两个。我需要在 rasterLayer 中对它们进行编号，使用数据框来链接波段编号和区域编号，或者使用可以返回哪些单元格编号属于每个区域的函数。最后，需要将至少 1 个波段中的数据所在的每个像素分配给这样的“区域”。

有点像你对多边形特征进行联合，但需要知道子区域共有哪个原始多边形。

@Benjamin：想添加一些样本数据作为测试用例吗？

@Joris Meys：完成。好主意。

【参考方案1】：

编辑：使用尼克的技巧和矩阵乘法更新答案。

---

您可以尝试使用尼克的技巧和矩阵乘法优化的以下函数。现在的瓶颈是用单独的层填充堆栈，但我想现在时间还不错。内存使用量要少一些，但考虑到您的数据和 R 的性质，我不知道您是否可以在不影响性能的情况下吃一点。

> system.time(T1 <- FindBands(myraster,return.stack=T))
   user  system elapsed 
   6.32    2.17    8.48 
> system.time(T2 <- FindBands(myraster,return.stack=F))
   user  system elapsed 
   1.58    0.02    1.59 
> system.time(results <- lapply(values, find_zones))
  Timing stopped at: 182.27 35.13 217.71

该函数返回具有图中存在的不同级别组合的 rasterStack（这不是所有可能的级别组合，因此您已经获得了一些收益），或者带有级别编号和级别名称的矩阵。这使您可以执行以下操作：

levelnames <- attr(T2,"levels")[T2]

获取每个单元格点的级别名称。如下所示，您可以轻松地将该矩阵放入 rasterLayer 对象中。

功能：

 FindBands <- function(x,return.stack=F)
    dims <- dim(x)
    Values <- getValues(x)
    nn <- colnames(Values)

    vec <- 2^((1:dims[3])-1)
    #Get all combinations and the names
    id <- unlist(
                lapply(1:10,function(x) combn(1:10,x,simplify=F))
              ,recursive=F)

    nameid <- sapply(id,function(i)
      x <- sum(vec[i])
      names(x) <- paste(i,collapse="-")
      x
    )
    # Nicks approach
    layers <- Values %*% vec
    # Find out which levels we need
    LayerLevels <- unique(sort(layers))
    LayerNames <- c("No Layer",names(nameid[nameid %in% LayerLevels]))

    if(return.stack)
        myStack <- lapply(LayerLevels,function(i)
          r <- raster(nr=dims[1],nc=dims[2])
          r[] <- as.numeric(layers == i)
          r
           )
        myStack <- stack(myStack)
        layerNames(myStack) <- LayerNames
        return(myStack)

     else 

      LayerNumber <- match(layers,LayerLevels)
      LayerNumber <- matrix(LayerNumber,ncol=dims[2],byrow=T)
      attr(LayerNumber,"levels") <- LayerNames
      return(LayerNumber)
        

概念证明，使用 RobertH 的数据：

r <- raster(nr=10, nc=10)
r[]=0
r[c(20:60,90:93)] <- 1
s <- list(r)
r[]=0
r[c(40:70,93:98)] <- 1
s <- c(s, r)
r[]=0
r[50:95] <- 1
s <- (c(s, r))
aRaster <- stack(s)


> X <- FindBands(aRaster,return.stack=T)
> plot(X)

> X <- FindBands(aRaster,return.stack=F)
> X
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    1    1    1    1    1    1    1    1    1     1
 [2,]    1    1    1    1    1    1    1    1    1     2
 [3,]    2    2    2    2    2    2    2    2    2     2
 [4,]    2    2    2    2    2    2    2    2    2     4
 [5,]    4    4    4    4    4    4    4    4    4     8
 [6,]    8    8    8    8    8    8    8    8    8     8
 [7,]    7    7    7    7    7    7    7    7    7     7
 [8,]    5    5    5    5    5    5    5    5    5     5
 [9,]    5    5    5    5    5    5    5    5    5     6
[10,]    6    6    8    7    7    3    3    3    1     1
attr(,"levels")
[1] "No Layer" "1"        "2"        "3"        "1-2"      "1-3"
       "2-3"      "1-2-3"   

> XX <- raster(ncol=10,nrow=10)
> XX[] <- X
> plot(XX)

【讨论】：

不错！这正是我的想法。我想现在唯一的问题是是否可以更快地做到这一点，和/或使用更少的内存。感谢您的回答。

@Benjamin：使用尼克斯方法和矩阵乘法进行了更新。如果这是您正在寻找的，谢谢尼克。

嘿，我不是让你睡过头的人 :-) 荣誉归你所有。

谢谢，这非常快，而且在内存上也相当不错（无论如何，return.stack=F）。

@Benjamin：显然。返回堆栈时，您在测试用例中创建了 34 个 1000x1000 的矩阵，这确实很好地填满了内存。该选项仅用于提取单独的图层。

【参考方案2】：

我不熟悉光栅，但从我上面的理解来看，你基本上有一个 10*3000*3000 的数组，对吧？

如果是这样，对于栅格中的每个位置（第二个和第三个索引，currow 和 curcol），您可以使用二进制计算其“区域”的唯一标识符：在“波段”（第一个索引）上运行 i 并总和 r[i,currow, curcol]*2^(i-1)。根据 raster 的内部工作原理，应该可以快速实现这一点。

这会生成一个大小为 3000*3000 的新“栅格”，其中包含每个位置的唯一标识符。找到其中的唯一值会为您返回数据中实际出现的区域，并且反转二进制逻辑应该为您提供属于给定区域的波段。

如果我对光栅的解释不正确，请原谅我：那么请忽略我的想法。无论哪种方式都不是一个完整的解决方案。

【讨论】：

这给了我一个像素上存在的波段数量的标识符，而不是哪些波段的唯一标识符。例如，当（仅）波段 3 和 10 与（仅）波段 4 和 5 存在时，我应该得到不同的标识符，但在你的情况下，我得到相同的数字。

我不这么认为：当波段 3 和 10 存在时，你会得到 2^2+2^9，而对于波段 4 和 5，你会得到 2^3+2^4。但再一次，我可能完全没有抓住重点。

抱歉，实际上效果很好，而且速度非常快。只需要将每一层重新赋值为 2^i，然后取所有波段的总和。现在我只需要弄清楚如何从那个数字中识别出哪些波段。酷！

【参考方案3】：

这个怎么样？

library(raster)
#setting up some data

r <- raster(nr=10, nc=10)
r[]=0
r[c(20:60,90:93)] <- 1
s <- list(r)
r[]=0
r[c(40:70,93:98)] <- 1
s <- c(s, r)
r[]=0
r[50:95] <- 1
s <- (c(s, r))
plot(stack(s))

# write a vectorized function that classifies the data
# 
fun=function(x,y,z)cbind(x+y+z==0, x==1&y+z==0, y==1&x+z==0, z==1&x+y==0, x==0&y+z==2, y==0&x+z==2, z==0&x+y==2,x+y+z==3)

z <- overlay(s[[1]], s[[2]], s[[3]], fun=fun)
# equivalent to
#s <- stack(s)
#z <- overlay(s[[1]], s[[2]], s[[3]], fun=fun)

ln <- c("x+y+z==0", "x==1&y+z==0", "y==1&x+z==0", "z==1&x+y==0", "x==0&y+z==2", "y==0&x+z==2", "z==0&x+y==2", "x+y+z==3")
layerNames(z) <- ln
x11()
plot(z)

更通用：

s <- stack(s)
fun=function(x)as.numeric(paste(which(x==1), collapse=""))
x <- calc(s,fun)

当 nlayers(s) 有两位数（“1”，“2”与“12”相同时，这并不好，在这些情况下，您可以改用下面的函数 (fun2)：

fun2=function(x)as.numeric(paste(c(9, x), collapse=""))
x2 <- calc(s,fun2)

unique(x)
# [1]   1   2   3  12  13  23 123

unique(x2)
# [1] 9000 9001 9010 9011 9100 9101 9110 9111

仅用于玩具示例：

plot(x)
text(x)
p=rasterToPolygons(x)
plot(p, add=T)

【讨论】：

感谢您的回答。当然这很快，它是一个只有 3 层的 10x10 栅格！有没有一种简单的方法可以扩展它以处理 n 层？我并不是要苛刻，只是这并不能完全解决问题（尽管它适用于这种特定情况）。

这个例子不是为了速度，而是为了说明你应该如何用光栅来解决这个问题。这个想法是编写自己的函数，如“fun”，然后将其与“raster”函数（如叠加或计算）一起使用。避免像 Joris 示例中的 'findBounds' 之类的函数，这些函数一次对整个栅格进行操作，并使用 getValues() 会导致内存问题。

【参考方案4】：

我为@Nick Sabbe 的建议编写了代码，我认为它非常简洁且相对较快。这假设输入 rasterStack 已经有逻辑 1 或 0 数据：

# Set the channels to 2^i instead of 1
bands = nlayers(myraster)
a = stack()
for (i in 1:bands) 
  a = addLayer(a, myraster[[i]] * 2^i)

coded = sum(a)
#plot(coded)
values = unique(coded)[-1]
remove(a, myraster)

# Function to retrieve which coded value means which channels
which_bands = function(value) 
  single = numeric()
  for (i in bands:1) 
    if ((0 < value) & (value >= 2^i)) 
     value = value - 2^i
      single = c(single, i)
    
  
  return(single)

【讨论】：

@Joris Meys：我不明白为什么这是必要的。有没有什么不同的情况？我以为他只是这样做，因为他认为有一个乐队 0...

对您的情况没有太大影响。但是，如果带的数量会变得更大，那么您很快就会遇到 R 中整数的限制，因此我建议将其保持在尽可能低的水平。由于（假设）没有波段 0，我可以将所有索引向下移动。

您可以使用矩阵乘法来加快速度。请参阅我的更新答案。

你也可以这样做：v