[R][r] 栅格操作,从缓冲区中提取值而不重叠

问题描述

我和我的研究伙伴正在撰写关于巴西蝙蝠的适宜性/种群数量的论文 我们正在努力操纵栅格以实现我们的方法, 模型中将使用大约200个栅格(不同物种的不同场景) 我们为这个问题选择并裁剪了一个

我们正在努力做的是:

  • 基于值高于阈值(在本例中为 0.8)的像素周围的缓冲区,从每个栅格中提取信息
  • 理想情况下,重叠的缓冲区只会选择一个像素来保留。创建一组生成非重叠缓冲区的像素。
  • 此外,对于下一步,我们希望根据当前值和线性模型更改栅格的值(这部分很好)
  • 值更改后,我们希望合并来自同一场景的栅格以寻找物种之间的相似性

我们一直在努力寻找一种不会产生重叠缓冲区或以更规范的方式聚合像素的方法

我们的栅格分辨率为 30 弧秒(大约 1 平方公里),我们只能聚合像素,但我们需要特定区域,并且不可能仅聚合整个像素生成值,例如 6 平方公里 另一种可能性是聚合部分像素以访问例如 6 平方公里的区域。可能吗?

到目前为止,我已经尝试了两种方法提取和聚合),但都没有返回所需的结果。以下是目前使用的代码

library(raster)
# library(ggplot2)
library(dplyr)

raster_values <- c(0.7151692,0.7234125,0.7242436,0.8838134,0.9855102,0.9921246,0.9679778,0.8245632,0.8965716,NA,0.721549,0.6988058,0.7333487,0.8138089,0.829727,0.8689544,0.8607966,0.9794912,0.9012381,0.7118917,0.7103891,0.7527786,0.7792872,0.8320968,0.8082156,0.920545,0.9788723,0.7859345,0.6824703,0.7039984,0.7589136,0.7905939,0.9024587,0.9848859,0.969553,0.8404503,0.9922802,0.6883243,0.731391,0.7682831,0.8586601,0.9850862,0.9705435,0.9888299,0.8164479,0.666971,0.757661,0.7628939,0.7868114,0.7910978,0.7650495,0.8227689,0.8148086,0.8691386,0.7376462,0.9176324,0.998813,0.7585487,0.7721,0.6508481,0.6098195,0.7708853,0.7119401,0.625409,0.6886432,0.8641906,0.9991203,0.7227083,0.6550816,0.5863016,0.7050957,0.6629267,0.6550342,0.6217013,0.8762864,0.9989462,0.6901366,0.7041199,0.6307223,0.6305411,0.7033732,0.7092581,0.7340803,0.7865254,0.9964261,0.7940102,0.6661926,0.7381653,0.6544684,0.6170949,0.641997,0.7506128,0.9248958,0.9903375,0.7662657,0.7847621,0.6582187,0.7361452,0.6488761,0.6309077,0.6542051,0.781707,0.940975,0.9350743,0.7824667,0.7215696,0.7388574,0.6753907,0.6716958,0.7162136,0.7920918,0.8702987,0.9929227,0.9775091,NA)

adeq_raster <- raster(nrows=12,ncols=12,xmn=-49,xmx=-48.5,ymn=-28,ymx=-27.5,crs = "+proj=longlat +datum=wgs84 +no_defs ",vals=raster_values)
adeq_raster_df <- as.data.frame(adeq_raster,xy = T)

# ggplot()+
#   geom_raster(data = adeq_raster_df,aes(x = x,y = y,fill = layer))+
#   scale_fill_viridis_c()+
#   coord_quickmap()


# - Extract method
#this method works fine,but does not avoid overlapping buffer.
adeq_raster_df_points <- filter(adeq_raster_df,layer >= 0.80) %>% 
  dplyr::select(x,y) #seleciona apenas as duas colunas de localização
#transfor the points df into a sp format to be used in the extract function and maintain the coordenates values.
adeq_raster_df_points_sp <- sp::SpatialPoints(coordinates(adeq_raster_df_points))

adeq_raster_extract <- extract(adeq_raster,adeq_raster_df_points_sp,buffer = 5000,fun = mean,na.rm=TRUE,sp = TRUE)
adeq_raster_extract_df <- as.data.frame(adeq_raster_extract)


# - Aggregate method
# - This method have not been so promissing,beucase it restrain the agregation to a certain number of pixels and we need buffers based in distance (3,5 or 10 km).
adeq_raster_5km <- aggregate(adeq_raster,2,fun=mean,expand=F,na.rm=TRUE)#,filename='output/adeq_raster_5km')
adeq_raster_5km_df <- as.data.frame(adeq_raster_5km,xy = T)

我们欢迎任何关于如何处理缓冲区的建议。

解决方法

请在 代码 中提供示例数据 --- 应该没有理由下载数据(请参阅 R 帮助文件和这些页面以获取 100 多个示例)。提出一个重点问题也很好。

我认为这是使用 terra 包获取您想要的缓冲区。你可以对raster和朋友做类似的事情。

示例数据:

library(terra)
r <- rast(nrow=20,ncol=20,xmin=0,xmax=1,ymin=0,ymax=1,crs="+proj=utm +zone=1 +datum=WGS84")
set.seed(123)
values(r) <- runif(ncell(r))
r <- ifel(r < .97,NA,r)

制作点和缓冲

p <- as.points(r)
b <- buffer(p,0.1)

找到不相交的缓冲区

x <- relate(b,relation="intersects")
i <- rowSums(x,na.rm=TRUE) == 1
bb <- b[i,]

或者你想要合并重叠缓冲区?在这种情况下,您可以使用聚合(可能后跟 disaggregate 以创建单个多边形)

aa <- aggregate(b)

插图

plot(r,xlim=c(-0.1,1.1),ylim=c(-0.1,1.1))
lines(b,lty=2)
lines(aa,col="blue",lwd=2)
lines(bb,col="red",lwd=3)

enter image description here

如果您现在想为每个聚合缓冲区采样一个单元

a <-disaggregate(aa)
set.seed(3) 
e <- extract(r,a,function(i) sample(na.omit(i),1))

==== 使用您的新示例数据 ====

你也可以这样想

library(terra)

v <- c(0.72,0.72,0.88,0.99,0.97,0.82,0.9,0.7,0.73,0.81,0.83,0.87,0.86,0.98,0.71,0.75,0.78,0.92,0.79,0.68,0.76,0.84,0.69,0.77,0.67,0.74,1,0.65,0.61,0.63,0.66,0.59,0.62,0.64,0.94,NA)
r <- rast(nrows=12,ncols=12,xmin=-49,xmax=-48.5,ymin=-28,ymax=-27.5,vals=v)
x <- r > 0.8

x <- ifel(r > 0.8,NA)
p <- disaggregate(as.polygons(x))
b <- buffer(p,.01)
plot(x)
lines(b)