3

For example, I would like to build a 20-m buffer (represented in yellow in the image below) outside polygons having a value of 2 (represented in green) and to assign the value of 3 to buffer cells.

enter image description here

I tested the function buffer (package raster) but all buffer cells (including those of original polygons) have the same value.

Here is a reproducible example to build the raster with cell polygons.

pol1 <- rbind(c(-159,-19), c(-158,3), c(-81, 3), c(-148,-75), c(-180,-32))
pol2 <- rbind(c(40,0), c(121,52), c(151,0), c(121,-58), c(81,0))
polys <- SpatialPolygons(list(Polygons(list(Polygon(pol1)), 1),  
                              Polygons(list(Polygon(pol2)), 2)))
proj4string(polys) = CRS("+init=epsg:32188") ## NAD83 / MTM zone 8
plot(polys)
r <- raster(ncol=180, nrow=180)
extent(r) <- extent(polys) + 100
projection(r) <- CRS("+init=epsg:32188") 
res(r) <- 10
rp <- rasterize(polys, r)
rp[is.na(rp)] <- 0
plot(rp)
Improve this question
1
  • Can you show us your buffer function call? I suspect all you really need to do is to combine the buffered raster with the original raster.
    – Spacedman
    Apr 7, 2016 at 7:09

2 Answers 2

Reset to default
6

Some simple raster arithmetic should sort this out for you.

First make a raster where its NA anywhere except where the original was equal to 2:

> rp2 = rp ; rp2[rp2[]!=2]=NA
> plot(rp2)

Now we can buffer that to 20m:

> rp2b = buffer(rp2, 20)
> plot(rp2b)

Now the ring of the buffer is where rp2 is NA and rp2b is not NA:

> rbuff = is.na(rp2) & !is.na(rp2b))
> plot(rbuff)

enter image description here

This should give you a raster which is 1 in the buffer border and zero inside and outside (see screenshot). You could then add this to your original raster or whatever you want to do with it. I think this does what you asked for, although the colours are slightly different because of the cells with value 1 as well:

> plot(rbuff*3 + rp)

enter image description here

If you want you can also do this for vector data by using the gDifference function in package:rgeos to subtract the original polygon from the buffered polygon.

Improve this answer
2

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:

Figure

I believe that's an order of magnitude faster than alternatives such as the buffer operation.

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.

Improve this answer
2
  • Thank you very much whuber for your answer. Your codes are very useful because I have a raster with 600000000 cells. I tested the two codes but the given value (for example, value of 3) is not assigned to the buffer cells. In my case, the cells around polygons have different values and not only a value of 0 (as in my example). So, how can I modify the condition ifelse(y > 0, y, ifelse(x > 0, 3, 0))) to do this ?
    – Pierre
    Apr 11, 2016 at 15:18
  • As you can see, the code worked in the test cases I constructed. The reason for its failure in your example therefore needs investigation and debugging. That the surrounding cells have nonzero values should not matter, so I suspect you might not have constructed the indicator grid correctly (in which all values outside the target polygons are converted to NA). Work with very small grids until you get this issue sorted out and working, and then experiment with larger grids: it will save you time.
    – whuber
    Apr 11, 2016 at 15:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.