Mask raster from data in shapefile

Question

I am trying to standardise a raster file by removing the values exceeding the mean +/- 2.5 std calculated for each polygon of a shapefile. I manage to calculate the zonal statistics and feed them back into the shapefile, but am unable to mask my raster file with the data included in the shapefile. I get an error message running the last two lines of code. Extracting the values above mean + 2.5 std gives me only a value (NA), rather than a raster as output and as a result the mask function gives the following error:

Error in (function (classes, fdef, mtable): unable to find an inherited method for function ‘mask’ for signature ‘"RasterBrick", "logical"’

#Read raster file
library(raster)
prec <- getData('worldclim', var='prec', res=10)

#Read shapefile
france <- getData('GADM', country='FRA', level=1)

#Convert shapefile to raster
r_fra <- rasterize(france, prec)

#Calculate mean for each polygon
mean_fra_prec <- zonal(prec, r_fra, FUN='mean')

#Calculate standard deviation for each polygon
std_fra_prec <- zonal(prec, r_fra, FUN = 'sd')

#Calculate mean +/- 2.5 std for each polygon
ul_prec  <- mean_fra_prec + (2.5 * std_fra_prec)
ll_prec <- mean_fra_prec - (2.5*std_fra_prec)

#Merge mean +/- 2.5 std for each polygon
x <- data.frame(zone=1:22)
ul_prec2 <- merge(x, ul_prec, by='zone', all.x=TRUE)
ll_prec2 <- merge(x, ul_prec, by='zone', all.x=TRUE)

#Join new data to shapefile
france@data <- cbind(france@data, ul_prec2[,-1])
france@data <- cbind(france@data, ll_prec2[,-1])

#Remove values above mean + 2.5 std for each polygon
ext_prec <- prec[!(prec > france@data$ul_prec2)] <- NA #Values above mean + 2.5 std
ma_prec <- mask(prec, ext_prec) #Mask values above mean + 2.5 std

What happens? error messages? data values not what you expected? — Carl Witthoft
Stretching the meaning of reproducible a bit here. Why not use rastr = raster(volcano) from raster package? — geotheory
Thanks for your comments. I have changed the code so it is now truly reproducible and states the error message I get. — mb_eco

Robert Hijmans Robert Hijmans · Accepted Answer · 2013-06-13T07:54:26

Here is working script, I think.

library(raster)
# Get raster data
prec <- getData('worldclim', var='prec', res=10)
# Get polygons
france <- getData('GADM', country='FRA', level=1)

# for this example, simplify to one layer, smaller extent
prec <- crop(prec[[1]], france)

# rasterize polygons
r_fra <- rasterize(france, prec)

# Calculate mean for each polygon
mean_fra_prec <- zonal(prec, r_fra, fun='mean')
# order and remove ID
mean_fra_prec <- mean_fra_prec[order(mean_fra_prec[,1]), -1]
#Calculate standard deviation for each polygon
std_fra_prec <- zonal(prec, r_fra, fun = 'sd')
std_fra_prec <- std_fra_prec[order(std_fra_prec[,1]), -1]

# alternative route, may fail on large data set
# v <- extract(prec, france)
# mean_fra_prec <- sapply(v, mean)
# std_fra_prec <- sapply(v, sd)

#Calculate mean +/- 2 std for each polygon
ul  <- mean_fra_prec + 2*std_fra_prec
ll <- mean_fra_prec - 2*std_fra_prec
d <- data.frame(1:length(ll), ll, ul)

v <- subs(r_fra, d, which=2:3)

keep <- prec > v[[1]] & prec < v[[2]]

x <- mask(prec, keep, maskvalue=FALSE)
plot(x)

Mask raster from data in shapefile

1 Answers