4

ポイントとポリゴンの 2 つのシェープファイルがあります。以下のコードではgCentroid()、rgeosパッケージから を使用してポリゴンの重心を計算し、重心の周りにバッファーをプロットします。

各セルから重心の周りの関連付けられたポリゴン バッファー内にある最も近いポイント (赤) までの距離を表すポリゴンからラスター レイヤーを作成したいと考えています。

たとえば、ポリゴン ユニットでは、AI は 2 つの仮想ラスター セルを表示し、計算しようとしている直線距離を示します。

ここに画像の説明を入力

更新 1 : @JMT2080AD のコメントに基づいて実際のバッファーを作成します。leafletコードの交換。

library(raster)
library(rgdal)
library(rgeos)

url <- "https://www.dropbox.com/s/25n9c5avd92b0zu/example.zip?raw=1"
download.file(url, "example.zip")
unzip("example.zip")

myPolygon <- readOGR("myPolygon.shp")
proj4string(myPolygon) <- CRS("+init=epsg:4326")
myPolygon <- spTransform(myPolygon, CRS("+proj=robin +datum=WGS84"))

myPoints <- readOGR("myPoints.shp")
proj4string(myPoints) <- CRS("+init=epsg:4326")
myPoints <- spTransform(myPoints, CRS("+proj=robin +datum=WGS84"))

centroids <- gCentroid(myPolygon, byid=TRUE)
buffer <- gBuffer(centroids, width=5000, byid=TRUE)

plot(myPolygon, col="green")
plot(buffer, col="blue", add = T)
plot(centroids, pch = 20, col = "white", add = T)
plot(myPoints, pch = 20, col = "red", add = T)

gis.stackexchangeでこの質問をしましたが、QGIS のコンテキストで。Rでこれを理解する方が良いと思うので、質問と新しいR MREをここに再投稿しています。質問をSOに移行してMREを変更するより良い方法があるかどうかはわかりません同時に。

4

3 に答える 3

4

sfを使用した別のソリューションを次に示します。私はこれにやや異なる方法でアプローチしています。ベクトルデータ表現を使用してすべての計算を行っており、結果のみをラスタライズしています。これは、ラスター セル内のどこからポイントまでの距離を測定するかが実際に重要であることを強調するために行っています。以下のコードは、各ラスター セルとターゲット ポイントの間の距離を測定する 2 つの方法を提供します。a) 最も近い頂点 (dist_pol) から、b) 重心 (dist_ctr) から。ターゲットの解像度に応じて、これらの違いは非常に大きくなる場合もあれば、無視できる場合もあります。以下のケースでは、約 100m x 100m のセル サイズで、差は平均して、セル エッジの長さに近くなります。

library(sf)
# library(mapview)
library(data.table)
library(raster)
# devtools::install_github("ecohealthalliance/fasterize")
library(fasterize)

url <- "https://www.dropbox.com/s/25n9c5avd92b0zu/example.zip?raw=1"
download.file(url, "/home/ede/Desktop/example.zip")
unzip("/home/ede/Desktop/example.zip")

pls = st_read("/home/ede/Desktop/example/myPolygon.shp")
pts = st_read("/home/ede/Desktop/example/myPoints.shp")
buf = st_read("/home/ede/Desktop/example/myBuffer.shp")

### extract target points within buffers
trgt_pts = st_intersection(pts, buf)

# mapview(pls) + buf + trgt_pts

### make grid and extract only those cells that intersect with the polygons in myPolygon.shp
grd_full = st_make_grid(pls, cellsize = 0.001) # 0.001 degrees is about 100 m longitude in Uganda
grd = grd_full[lengths(st_intersects(grd_full, pls)) > 0]

### do the distance calculations (throughing in some data.rable for the performance & just because)
### dist_pol is distance to nearest polygon vertex
### dist_ctr is distance to polygon centroid
grd = as.data.table(grd)
grd[, pol_id := sapply(st_intersects(grd$geometry, pls$geometry), "[", 1)]
grd[, dist_pol := apply(st_distance(geometry, trgt_pts$geometry[trgt_pts$id.1 %in% pol_id]), 1, min), by = "pol_id"]
grd[, dist_ctr := apply(st_distance(st_centroid(geometry), trgt_pts$geometry[trgt_pts$id.1 %in% pol_id]), 1, min), by = "pol_id"]

### convert data.table back to sf object
grd_sf = st_as_sf(grd)

### finally rasterize sf object using fasterize (again, very fast)
rast = raster(grd_sf, res = 0.001)
rst_pol_dist = fasterize(grd_sf, rast, "dist_pol", fun = "first")
rst_ctr_dist = fasterize(grd_sf, rast, "dist_ctr", fun = "first")

# mapview(rst_ctr_dist)

plot(rst_ctr_dist)
plot(stack(rst_pol_dist, rst_ctr_dist)) # there are no differences visually

### check differences between distances from nearest vertex and centroid
summary(grd_sf$dist_pol - grd_sf$dist_ctr)
于 2018-01-17T21:10:50.383 に答える