Home > Non classé, R > Extracting circular buffers from a raster in R

Extracting circular buffers from a raster in R

Cropping circular buffers centered on a point is pretty easy in R, at least when you’re dealing with rasters. Two packages allow to do this task, if not more. Let’s see how we can do with the package raster.

Note that in the following example, we assume that the package adehabitat is properly installed.


library(raster)
f<-system.file("ascfiles/elevation.asc", package="adehabitat")
raster<-rasterbis<-raster(f)# read example file from adehabitat folder
plot(raster)# have a look
rasterbis[]<-NA # empty one of the rasters
points<-data.frame(704000, 3162000);names(points)<-c("x","y")# create coordinates of a sample point
cell<-cellFromXY(object=raster,xy=points[1,])# get the cell number corresponding to that sample point
radius<-1000 # set the buffer radius
rasterbis[cell]<-1 # set the value of the corresponding cell to 1 (the others remaining NAs)
plot(raster)
points(points$x,points$y,pch=3)# one more look...
buf <- buffer(rasterbis, width=radius) # buffer extraction per se
plot(buf)
p<-which(buf[]==1) # identify which cells belong to the buffer
buf[p]<-raster[p]# pick up original values
# Now we crop the buffer :
xmin<-points$x[1] - 2*radius
xmax<-points$x[1] + 2*radius
ymin<-points$y[1] - 2*radius
ymax<-points$y[1] + 2*radius
ev<-c(xmin, xmax, ymin, ymax)
e<-extent(ev)
buffer<-crop(buf,e) # that's it :
plot(buffer)
points(points$x,points$y,pch=3)

I hope this might be helpful !
Cheers.
Next post : same thing with the package adehabitat

————————————–
EDIT: Robert Hijmans send me two sets of R commands which are more compact and very efficient ! Thank you Robert ! I add this code below so interested readers can make use of it. Thanks also to Curlew for the contribution (see comments below).

————————————–

library(raster)
library(rgeos)
points <- SpatialPoints(cbind(704000, 3162000))
pbuf <- gBuffer(points, widt=1000)
buf <- mask(raster, pbuf)
buffer <- trim(buf, pad=2)

library(raster)
raster <- raster(system.file(“ascfiles/elevation.asc”, package=”adehabitat”))
r1 <- raster(f)
points <- cbind(704000, 3162000)
r2 <- rasterize(points, r1, field=1)
pbuf <- buffer(r2, width=1000)
buf <- mask(raster, pbuf)
buffer <- trim(buf, pad=2)

 

Advertisements
Categories: Non classé, R
  1. 2 February 2014 at 2109 23

    Nice!
    Another way which is maybe more straight forward is to generate the points buffer in R (using the rgeos package) and then use the buffer to mask the raster’s values to a new one:

    library(rgeos)
    sp <- SpatialPoints(points) # make spatialpoints object of your points file
    buf <- gBuffer(sp,width=1000,quadsegs=50) # create your buffer
    raster2 <- crop(raster, buf, snap="out") # First crop the raster
    crop<- setValues(raster2, NA) # Create an empty dummy raster for the rasterization
    bufr <- rasterize(buf, crop) # rasterize the generated buffer
    out <- mask(x=raster2, mask=bufr) # Create your mask by putting NA in all cells outside the boundary
    plot(mask)

  1. No trackbacks yet.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: