2017-06-27 13 views
1

次の問題があります。私はCH(Fribourg)の州のデータを持っており、この地域で作成されたグリッドセルと3つの円を交差させたいと考えています。ここで添付の地図やポリゴンを生成するためにいくつかのコードされています。私は基本的に何をしたいのか円で交差する格子

data file

library(raster) 

centroid.x <- c(566052, 576768, 560652) 
centroid.y <- c(155501, 185501, 165325)` 

centroids <- data.frame(centroid.x = centroid.x, centroid.y = centroid.y) 

radius <- 1000 
n = 1000 
pts = seq(0, 2 * pi, length.out = n) 
ZH = cbind(centroid.x[1] + radius * sin(pts), centroid.y[1] + radius * cos(pts)) 
WH <- cbind(centroid.x[2] + radius * sin(pts), centroid.y[2] + radius * cos(pts)) 
GL <- cbind(centroid.x[3] + radius * sin(pts), centroid.y[3] + radius * cos(pts)) 
sl = SpatialPolygons(list(Polygons(list(Polygon(ZH)), 1), Polygons(list(Polygon(WH)), 2) 
          Polygons(list(Polygon(GL)), 3))) ` 

CRS_CH <- CRS("+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs") 
crs(sl) <- CRS_CH 
GRID <- 1000 
grid <- raster(extent(pol)) 

res(grid) <- GRID 
proj4string(grid)<-proj4string(pol) 

gridpolygon <- rasterToPolygons(grid) 
plot(gridpolygon) 

dry.grid1000 <- intersect(pol, gridpolygon) 
plot(dry.grid1000) 
plot(sl, add = TRUE, border = "red") 

enter image description here

は、円の外側のすべてのグリッドセルの値1を取得し、そしてあります2つの円の内部のグリッドセル。少し難しいのは、サークル全体にないグリッドセルです。これらのグリッドセルでは、重み付けされた平均値1と2をそれぞれ円の外側と内側の領域で重み付けしたいと考えています。例えば

は、私は、このセルが撮りたい

enter image description here

(赤が内側円と青の外側の部分である)我々は、以下の画像のようにサークルとグリッドセルの重なりを持っていることを言うことができます値(1 * | A | + 2 * | B |)/ 2

誰でも構いません。

答えて

1

ライブラリrasterに関数rasterizeを使用できます。パラメータはgetCover=TRUEです。

TRUEの場合、ポリゴンで覆われている各グリッドセルの割合が返されます(field、fun、mask、およびupdateの値は無視されます)。サブセルとあなたがより正確に何かをしたい場合は、より正確なラスター(100よりも高く、に自分でラスタをdisaggregateできる各々の中心にあるポリゴンの存在/不在は

サブセル決定、それ以外の場合は同じですgetCover)、元のラスタと分解されたラスタの対応を取得します。