r - 计算从多边形的所有部分到最近点的距离

Question

我有两个 shapefile:点和多边形。在下面的代码中，我使用 gCentroid()来自 rgeos包来计算多边形质心，然后我在质心周围绘制一个缓冲区。

我想从多边形创建一个栅格图层，该图层表示从每个像元到位于质心周围相关多边形缓冲区内的最近点(红色)的距离。

例如，在多边形单元 A 中，我显示了两个假设的栅格单元，并指示我要计算的直线距离。



更新 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 MRE 因为我认为我有更好的机会在 R 中解决这个问题。我不知道是否有更好的方法将问题迁移到 SO 并更改 MRE同时。

Answer 1

这是我的解决方案。我正在使用 sf只要有可能。根据我的经验sf与 raster 不完全兼容功能，所以这里有一些不太难看的解决方法。

我使用的基础数据与您提供的不同。

基础数据

library(sf)
library(raster)
library(magrittr)

set.seed(1)

## We will create your polygons from points using a voronoi diagram
x <- runif(10, 640000, 641000)
y <- runif(10, 5200000, 5201000)

myPolyPoints <- data.frame(id = seq(x), x = x, y = y) %>%
    st_as_sf(coords = c("x", "y"))

## Creating the polygons here
myPolygons <- myPolyPoints$geometry %>%
    st_union %>%
    st_voronoi %>%
    st_collection_extract

myPolygons <- st_sf(data.frame(id = seq(x), geometry = myPolygons)) %>%
    st_intersection(y = st_convex_hull(st_union(myPolyPoints)))

## Creating points to query with buffers then calculate distances to
polygonExt <- extent(myPolygons)
x <- runif(50, polygonExt@xmin, polygonExt@xmax)
y <- runif(50, polygonExt@ymin, polygonExt@ymax)

myPoints <- data.frame(id = seq(x), x = x, y = y) %>%
    st_as_sf(coords = c("x", "y"))

## Set projection info
st_crs(myPoints) <- 26910
st_crs(myPolygons) <- 26910

## View base data
plot(myPolygons$geometry)
plot(myPoints$geometry, add = T, col = 'blue')

## write out data
saveRDS(list(myPolygons = myPolygons,
             myPoints = myPoints),
        "./basedata.rds")

我生成的基础数据如下所示:



距离处理

library(sf)
library(raster)
library(magrittr)

## read in basedata
dat <- readRDS("./basedata.rds")

## makeing a grid of points at a resolution using the myPolygons extent
rast <- raster(extent(dat$myPolygons), resolution = 1, vals = 0, crs = st_crs(dat$myPoints))

## define a function that masks out the raster with each polygon, then
## generate a distance grid to each point with the masked raster
rastPolyInterDist <- function(maskPolygon, buffDist){
    maskPolygon <- st_sf(st_sfc(maskPolygon), crs = st_crs(dat$myPoints))
    mRas <- mask(rast, maskPolygon)
    cent <- st_centroid(maskPolygon)
    buff <- st_buffer(cent, buffDist)
    pSel <- st_intersection(dat$myPoints$geometry, buff)

    if(length(pSel) > 0){
        dRas <- distanceFromPoints(mRas, as(pSel, "Spatial"))
        return(dRas + mRas)
    }
    return(mRas)
}

dat$distRasts <- lapply(dat$myPolygons$geometry,
                        rastPolyInterDist,
                        buffDist = 100)

## merge all rasters back into a single raster
outRast <- dat$distRasts[[1]]

mergeFun <- function(mRast){
     outRast <<- merge(outRast, mRast)
}

lapply(dat$distRasts[2:length(dat$distRasts)], mergeFun)

## view output
plot(outRast)
plot(dat$myPoints$geometry, add = T)
dat$myPolygons$geometry %>%
    st_centroid %>%
    st_buffer(dist = 100) %>%
    plot(add = T)

结果如下所示。您可以看到，当缓冲质心不与其多边形中找到的任何位置相交时，会处理一个条件。



使用您的基础数据，我对如何在 R 中读取和处理数据进行了以下编辑。

OP基础数据

library(raster)
library(sf)
library(magrittr)

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

myPolygons <- st_read("myPolygon.shp") %>%
    st_transform(st_crs("+proj=robin +datum=WGS84"))

myPoints <- st_read("myPoints.shp") %>%
    st_transform(st_crs("+proj=robin +datum=WGS84"))

centroids <- st_centroid(myPolygons)
buffer <- st_buffer(centroids, 5000)

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

saveRDS(list(myPoints = myPoints, myPolygons = myPolygons), "op_basedata.rds")

使用 OP 数据进行距离处理

要使用我建议的计算程序，您只需要修改起始栅格的分辨率和缓冲区距离输入。否则，一旦您将数据读入 R 后，它的行为应该相同，正如我上面概述的那样。

library(sf)
library(raster)
library(magrittr)

## read in basedata
dat <- readRDS("./op_basedata.rds")

## makeing a grid of points at a resolution using the myPolygons extent
rast <- raster(extent(dat$myPolygons), resolution = 100, vals = 0, crs = st_crs(dat$myPoints))

## define a function that masks out the raster with each polygon, then
## generate a distance grid to each point with the masked raster
rastPolyInterDist <- function(maskPolygon, buffDist){
    maskPolygon <- st_sf(st_sfc(maskPolygon), crs = st_crs(dat$myPoints))
    mRas <- mask(rast, maskPolygon)
    cent <- st_centroid(maskPolygon)
    buff <- st_buffer(cent, buffDist)
    pSel <- st_intersection(dat$myPoints$geometry, buff)

    if(length(pSel) > 0){
        dRas <- distanceFromPoints(mRas, as(pSel, "Spatial"))
        return(dRas + mRas)
    }
    return(mRas)
}

dat$distRasts <- lapply(dat$myPolygons$geometry,
                        rastPolyInterDist,
                        buffDist = 5000)

## merge all rasters back into a single raster
outRast <- dat$distRasts[[1]]

mergeFun <- function(mRast){
     outRast <<- merge(outRast, mRast)
}

lapply(dat$distRasts[2:length(dat$distRasts)], mergeFun)

## view output
plot(outRast)
plot(dat$myPoints$geometry, add = T)
dat$myPolygons$geometry %>%
    st_centroid %>%
    st_buffer(dist = 5000) %>%
    plot(add = T)