这是我的解决方案。我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)
