9

パッケージ：

データ：

  • 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...

サンプルデータがどのように見えるかを示す

4

4 に答える 4

7

編集：ニックのトリックと行列の乗算を使用して更新された回答。

Nickのトリックと行列の乗算を使用して最適化された次の関数を試すことができます。ボトルネックは現在、スタックを個別のレイヤーで埋めていますが、タイミングは今ではかなり問題ないと思います。メモリ使用量は少し少なくなりますが、データと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)

ここに画像の説明を入力してください

于 2011-04-18T14:31:39.640 に答える
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が2桁の場合（「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)
于 2011-04-18T00:59:32.900 に答える
3

私はラスターに精通していませんが、上記から把握すると、基本的に10 * 3000 * 3000の配列がありますよね？

その場合、ラスター内の各位置（2番目と3番目のインデックス、currowとcurcol）について、バイナリを使用して「ゾーン」の一意の識別子を計算できます。「バンド」（最初のインデックス）に対してiを実行し、r[を合計します。 i、currow、curcol] * 2 ^（i-1）。ラスターの内部動作によっては、これをかなり迅速に実装できるはずです。

これにより、各位置の一意の識別子を保持するサイズ3000*3000の新しい「ラスター」が作成されます。そこで一意の値を見つけると、データで実際に発生するゾーンが返され、バイナリロジックを逆にすると、特定のゾーンに属するバンドが得られます。

ラスターの解釈が正しくない場合はご容赦ください。それでは、私の考えを無視してください。いずれにせよ、完全な解決策ではありません。

于 2011-04-18T14:09:50.987 に答える
1

@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)
}
于 2011-04-18T20:26:04.307 に答える