Here is my solution. I am using sf whenever possible. From my experience sf is not completely compatible with the raster functions yet, so there are a couple of workarounds here that aren't too ugly.
I am using different base data than what you provided.
Base Data
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")
The base data I generated looks like so:
Distance Processing
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)
The results can be seen below. You can see that there is a condition handled when the buffered centroid does not intersect any locations found in its polygon.
Using your base data, I have made the following edits to how your data is read and processed in R.
OP Base Data
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")
Distance Processing Using OP Data
To use the calculation routine I have suggested, you just need to modify the resolution of the starting raster, and the buffer distance input. Othewise it should behave the same once you read your data into R as I have outlined above.
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)
