Маска растра из данных в шейп-файле

Question

Я пытаюсь стандартизировать растровый файл, удаляя значения, превышающие средний +/- 2.5 std, рассчитанный для каждого многоугольника шейп-файла. Мне удается рассчитать статистику по зонам и вернуть их обратно в шейп-файл, но я не могу замаскировать свой растровый файл данными, содержащимися в шейп-файле. Я получаю сообщение об ошибке, выполняющее последние две строки кода. Извлечение значений выше среднего + 2,5 std дает мне только значение (NA), а не растр в качестве результата, и в результате функция маски дает следующую ошибку:

Ошибка в (function (classes, fdef, mtable) : не удалось найти наследуемый метод для функции «маски» для подписания «RasterBrick", „логический“»

#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 

Answer 1

Здесь работает сценарий, я думаю

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)