Archive

Archive for the ‘Non classé’ Category

Extracting circular buffers from a raster in R

31 January 2014 1 comment

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