1

I have raster of the following features:

library(raster)
library(rgeos)
test <- raster(nrow=225, ncols=478, xmn=-15.8, xmx=32, ymn=-9.4, ymx=13.1)

I want to mask in this raster the cells that are within a given distance of a point. I create the spatial points as followed:

p2=readWKT("POINT(31.55 -1.05)")

Then I create a spatial polygon object by adding a 0.5 buffer:

p2_Buffered <- gBuffer(p2, width = 0.5)

mask(test, mask=p2_Buffered,inverse=T)

When I mask my raster given this spatial object, I have the following error message:

Error in .polygonsToRaster(x, y, field = field, fun = fun, background = background, : number of items to replace is not a multiple of replacement length

I do not understand because this is script I have been running many many times with different point and different buffer width without any problem.

What is strange is that when I change the width of the buffer, it works fine:

p2_Buffered <- gBuffer(p2, width = 0.4)
mask(test, mask=p2_Buffered,inverse=T)

This is also true for a different focal point:

p2=readWKT("POINT(32.55 -1)")
p2_Buffered <- gBuffer(p2, width = 0.5)
mask(test, mask=p2_Buffered,inverse=T)

I would like to identify the specific problem I have for that point because this is a script I should run in a routine (I have been doing it without any problem so far).

Thanks a lot

dauby gilles
  • 57
  • 1
  • 5
  • From where come functions `raster`, `readWKT`, `gBuffer`, `mask`? –  Apr 08 '15 at 08:59
  • packages `rgeos` and `raster` I assume. @dauby gilles you should edit your post and mention the packages you used in your code by including the `library("...")` command – rmuc8 Apr 08 '15 at 09:19
  • Maybe a bug. As you noticed, it works with "31.54" and "31.56", but not with "31.55" –  Apr 08 '15 at 09:58

2 Answers2

3

This is indeed a bug with polygons that go over the edge of a raster. It has been fixed in version 2.3-40 (now on CRAN), so it should go away if you update the raster package.

Here is a workaround (removing the part of the polygon that goes over the edge). 

library(raster)
library(rgeos)
r <- raster(nrow=225, ncols=478, xmn=-15.8, xmx=32, ymn=-9.4, ymx=13.1)
e <- as(extent(r), 'SpatialPolygons')

p <- readWKT("POINT(31.55 -1.05)")   
pb <- gBuffer(p, width = 0.5)
pbe <- intersect(pb, e)

values(r) 
x <- mask(r, mask=pbe, inverse=TRUE)
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
1

You usually need to set some values to the raster layer. For a mask layer its always best to set values to 1. 

    library(raster)
    library(rgeos)

# make sample raster    
test <- raster(nrow=225, ncols=478, xmn=-15.8, xmx=32, ymn=-9.4, ymx=13.1)
# set values of raster for mask
test <- setValues(test, 1)

# make point buffer
p2=readWKT("POINT(15 5)")
p2_Buffered <- gBuffer(p2, width = 1.5)
# name projection of buffer (assume its the same as raster)
projection(p2_Buffered) <- projection(test)

# visual check 
plot(test); plot(p2_Buffered, add=T)

If you want to trim down your raster layer to the just the single polygon then try this workflow. 

step1 <- crop(test, p2_Buffered) # crop to same extent
step2 <- rasterize(p2_Buffered, step1) # rasterize polygon 
final <- step1*step2 # make your final product
plot(final)

If you just want to poke a hole in your raster layer then use the mask function 

# rasterize your polygon
p2_Buffered <- rasterize(p2_Buffered, test, fun='sum')

# now mask it 
my_mask <- mask(test, mask=p2_Buffered,inverse=T) # try changing the inverse argument 
plot(my_mask)
Matthew Bayly
  • 556
  • 5
  • 7