3

サンプルデータ

x <- raster(x=matrix(data=1:36, nrow=6), xmn=-1000, xmx=1000, ymn=-100, ymx=900)
x[c(8, 15, 16, 17, 22, 25, 26, 30, 31)] <- NA
plot(x)

ここに画像の説明を入力

問題

ラスター内の穴、つまりセル c(15:17, 22) で囲まれた領域を、穴ではない他のギャップ (つまり、残りの空のセル) と (アルゴリズム的に) 区別するにはどうすればよいですか?

これにより、ラスターの穴/島の領域でのみ操作を実行したり、カスタム値で穴を埋めたりすることが可能になります。

実際のラスターには約 30000 の穴があるため、速度が重要です。R とグラス GIS ソリューションの両方に興味があります。ご協力ありがとうございました。

4

3 に答える 3

3

ここでは、ポリゴン化を使用しないソリューション: (エレガントではありませんが、機能します)。ただし、穴/島を値 (つまり 999) に再分類し、他のすべての非島を NA に再分類する必要があります。このような：

x <- raster(x=matrix(rep(NA,36), nrow=6), xmn=-1000, xmx=1000, ymn=-100, ymx=900)
x[c(8, 15, 16, 17, 22, 25, 26, 30, 31)] <- 999

plot(x)

初期ダミー ラスター

次に、この関数を使用しclump()て島があるかどうかを確認します。その関数の優れた点は、これらの島の ID も返すことです。

#Get Islands with IDs
cl <- clump(x,directions=8)
plot(cl)

ID付きの塊/島

次に、島の周波数からデータフレームを作成します (これは、各島の ID を取得するためのものです)。

freqCl <- as.data.frame(freq(cl))

#remove the (row) which corresponds to the NA values (this is important for the last step)
freqCl <- freqCl[-which(is.na(freqCl$value)),]

島が国境に接しているかどうかを確認します。

#Check if the island touches any border and therefore isn't a "real island" (first and last column or row)

noIslandID <- c()
#First row
if(any(rownames(freqCl) %in% cl[1,])){
  eliminate   <- rownames(freqCl)[rownames(freqCl) %in% cl[1,]]
  noIslandID  <- append(noIslandID, eliminate)
}
#Last row
if(any(rownames(freqCl) %in% cl[nrow(cl),])){
  eliminate  <- rownames(freqCl)[rownames(freqCl) %in% cl[nrow(cl),]]
  noIslandID <- append(noIslandID, eliminate)
}
#First col
if(any(rownames(freqCl) %in% cl[,1])){
  eliminate   <- rownames(freqCl)[rownames(freqCl) %in% cl[,1]]
  noIslandID  <- append(noIslandID, eliminate)
}
#Last col
if(any(rownames(freqCl) %in% cl[,ncol(cl)])){
  eliminate   <- rownames(freqCl)[rownames(freqCl) %in% cl[,ncol(cl)]]
  noIslandID  <- append(noIslandID, eliminate)
}

国境に接する島を排除する:

noIslandID <- unique(noIslandID)
IslandID   <- setdiff(rownames(freqCl), noIslandID)

最後のステップで、最初のラスターからすべての「実際の島」に 1 を割り当てます。 

for(i in 1:length(IslandID)) {
  x[cl[]==as.numeric(IslandID[i])] <- 1
}

plot(x)

ここに画像の説明を入力

于 2016-12-23T23:40:05.627 に答える
2

ポリゴン化はどうですか？速度については、それが何の価値があるかわかりませんが、次のことができます。

x[!is.na(values(x))]<-1
plot(x)
x[is.na(values(x))]<-0

hole <- rasterToPolygons(x, fun=NULL, n=4, na.rm=TRUE, digits=12, dissolve=T)

次に、シングルパート ポリゴンをマルチパートにする必要があります。

hole2 <- SpatialPolygons(lapply(hole@polygons[[1]]@Polygons, function(xx) Polygons(list(xx),ID=round(runif(1,1,100000000)))))

plot(hole2, add=T)

これで、境界に触れていない「本当の」穴が見つかりました。

around <- as(extent(x), "SpatialLines")
touch_border <- gIntersects(around, hole2, byid=T) 
extract(x, hole2[!touch_border,],cellnumbers=T)

これにより、穴ごとのセル数がわかります。穴とは言っていないセル「8」が見つかったので、正しいかどうかはわかりませんが、非常に近いに違いありません!

R で遅い場合は、GRASS で同じアルゴリズムを実行します (主にrasterToPolygons呼び出し)

于 2016-12-22T19:57:48.187 に答える