1

我试图通过删除超过为 shapefile 的每个多边形计算的平均值 +/- 2.5 std 的值来标准化光栅文件。我设法计算了区域统计数据并将它们反馈回 shapefile，但无法使用 shapefile 中包含的数据屏蔽我的栅格文件。我收到一条运行最后两行代码的错误消息。提取上面的值 mean + 2.5 std 只给我一个值（NA），而不是一个栅格作为输出，因此 mask 函数给出以下错误：

（函数（类，fdef，mtable）中的错误：无法为签名'“RasterBrick”，“逻辑”'找到函数“掩码”的继承方法</p> 

#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
4

1 回答 1

2

我认为这是工作脚本。

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)
于 2013-06-13T07:54:26.200 回答