Simply overlay a reclassified focal mean or distance grid of the polygon indicator.
The focal mean requires a circular neighborhood
w. Here is a way to create it in terms of the radius, 20. It starts with constant values (line 4). Values beyond the desired radius are zeroed out (line 5). The result is normalized to sum to unity (line 6).
radius <- ceiling(20 / min(res(r)))
diameter <- 2*radius + 1
i <- outer(seq(-radius,radius), seq(-radius,radius), function(a,b) a^2+b^2) > radius^2
w <- matrix(1, diameter, diameter)
w[i] <- 0
w <- w / sum(w)
With this in hand, the raster operations are fast:
r.focal <- focal(rp==2, w)
result <- overlay(r.focal, rp, fun=function(x, y) ifelse(y > 0, y, ifelse(x > 0, 3, 0)))
The first line computes the focal mean of the indicator. It will be nonzero at any cell for which the polygon is within the radius. The second line reclassifies the nonzero values of the focal mean to 3 and overlays those on the cells outside the polygon.
For tiny test grids as in the question, this method doesn't look superior to others that might be suggested because of the overhead of creating the neighborhood matrix
w. For any substantial grid, though, it's likely to be pretty fast. For instance, changing the resolution of
r from 10 to 1 makes the grid 100 times larger (it now has almost 100,000 cells). The total elapsed time for the entire operation is 0.4 to 0.5 seconds:
I believe that's an order of magnitude faster than alternatives such as the
focal operation in the
raster package does not scale optimally. On a grid with ten times as many cells (almost a million), this focal sum solution takes almost half a minute: about 60 times faster. Here's an alternative that appears to scale better for this particular package. It is a solution based on computing Euclidean distances. To do this, all the cells outside the buffered polygon must first be set to NA. Then the distance grid is computed, limited to values less than the buffer radius, and overlaid as before.
r.0 <- rp
r.0[r.0 != 2] <- NA
result <- overlay(distance(r.0) <= 20, rp, fun=function(x,y) ifelse(y > 0, y, 3*x))
The map looks identical to the previous one. This calculation took 0.6 seconds on the 100,000 cell grid and 12 seconds on the 1,000,000 cell grid. Still not great scaling, but obviously better than the previous solution.
Expert users of the
raster package may be able to suggest slightly faster ways to carry out these individual operations.