[R][r] 栅格操作，从缓冲区中提取值而不重叠
[R][r] raster manipulation, extract values from buffer area without overlap
我和我的研究伙伴正在撰写一篇关于巴西蝙蝠的 suitability/population 论文
我们正在努力操纵栅格来实现我们的方法，
模型中将使用大约200个栅格（不同物种的不同场景）
我们select针对这个问题编辑并裁剪了一个
我们正在尝试做的是：
- 根据值高于阈值（在本例中为 0.8）的像素周围的缓冲区从每个栅格中提取信息
- 理想情况下，重叠的缓冲区将 select 仅保留一个像素。创建一组生成非重叠缓冲区的像素。
- 此外，对于下一步，我们想根据当前值和线性模型更改栅格的值（这部分很好）
- 值更改后，我们想合并来自同一场景的栅格以寻找物种之间的相似性
我们一直在努力寻找一种不会产生重叠缓冲区或以更规范的方式聚合像素的方法。
我们的栅格分辨率为 30 弧秒（大约 1 平方公里），我们只能聚合像素，但我们需要特定区域，不可能只聚合整个像素生成值，例如 6 平方公里
另一种可能性是聚合部分像素以访问 6 平方公里等区域。可能吗？
到目前为止，我已经尝试了两种方法（提取和聚合）并且 none 返回了所需的结果。以下是目前使用的代码：
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, NA, NA, 0.721549, 0.6988058, 0.7333487, 0.8138089, 0.829727, 0.8689544, 0.8607966, 0.9794912, 0.9012381, NA, NA, NA, 0.7118917, 0.7103891, 0.7527786, 0.7792872, 0.8320968, 0.8082156, 0.920545, 0.9788723, 0.7859345, NA, NA, NA, 0.6824703, 0.7039984, 0.7589136, 0.7905939, 0.9024587, 0.9848859, 0.969553, 0.8404503, NA, NA, NA, 0.9922802, 0.6883243, 0.731391, 0.7682831, 0.8586601, 0.9850862, 0.9705435, 0.9888299, 0.8164479, 0.666971, NA, NA, 0.757661, 0.7628939, 0.7868114, 0.7910978, 0.7650495, 0.8227689, 0.8148086, 0.8691386, 0.7376462, 0.9176324, NA, NA, 0.998813, 0.7585487, 0.7721, 0.6508481, 0.6098195, 0.7708853, 0.7119401, 0.625409, 0.6886432, 0.8641906, NA, NA, 0.9991203, 0.7227083, 0.6550816, 0.5863016, 0.7050957, 0.6629267, 0.6550342, 0.6217013, 0.8762864, 0.9989462, NA, NA, NA, 0.6901366, 0.7041199, 0.6307223, 0.6305411, 0.7033732, 0.7092581, 0.7340803, 0.7865254, 0.9964261, 0.7940102, NA, NA, 0.6661926, 0.7381653, 0.6544684, 0.6170949, 0.641997, 0.7506128, 0.9248958, 0.9903375, 0.7662657, 0.7847621, NA, NA, 0.6582187, 0.7361452, 0.6488761, 0.6309077, 0.6542051, 0.781707, 0.940975, 0.9350743, 0.7824667, NA, NA, NA, 0.7215696, 0.7388574, 0.6753907, 0.6716958, 0.7162136, 0.7920918, 0.8702987, 0.9929227, 0.9775091, NA, NA, 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)
我们欢迎任何有关如何处理缓冲区的建议。
请在 code 中提供示例数据 --- 应该没有理由下载数据（请参阅 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)
如果您现在想为每个聚合缓冲区采样一个单元格
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.72, 0.88, 0.99, 0.99, 0.97, 0.82, 0.9, NA, NA, NA, 0.72, 0.7, 0.73, 0.81, 0.83, 0.87, 0.86, 0.98,0.9,NA,NA,NA,0.71,0.71,0.75,0.78,0.83,0.81,0.92,0.98,0.79,NA,NA,NA,0.68,0.7,0.76,0.79,0.9,0.98,0.97,0.84,NA,NA,NA,0.99,0.69,0.73,0.77,0.86,0.99,0.97,0.99,0.82,0.67,NA,NA,0.76,0.76,0.79,0.79,0.77,0.82,0.81,0.87,0.74,0.92,NA,NA,1,0.76,0.77,0.65,0.61,0.77,0.71,0.63,0.69,0.86,NA,NA,1,0.72,0.66,0.59,0.71,0.66,0.66,0.62,0.88,1,NA,NA,NA,0.69,0.7,0.63,0.63,0.7,0.71,0.73,0.79,1,0.79,NA,NA,0.67,0.74,0.65,0.62,0.64,0.75,0.92,0.99,0.77,0.78,NA,NA,0.66,0.74,0.65,0.63,0.65,0.78,0.94,0.94,0.78,NA,NA,NA,0.72,0.74,0.68,0.67,0.72,0.79,0.87,0.99,0.98,NA,NA,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, 1, NA)
p <- disaggregate(as.polygons(x))
b <- buffer(p, .01)
plot(x)
lines(b)
我和我的研究伙伴正在撰写一篇关于巴西蝙蝠的 suitability/population 论文 我们正在努力操纵栅格来实现我们的方法， 模型中将使用大约200个栅格（不同物种的不同场景） 我们select针对这个问题编辑并裁剪了一个
我们正在尝试做的是：
- 根据值高于阈值（在本例中为 0.8）的像素周围的缓冲区从每个栅格中提取信息
- 理想情况下，重叠的缓冲区将 select 仅保留一个像素。创建一组生成非重叠缓冲区的像素。
- 此外，对于下一步，我们想根据当前值和线性模型更改栅格的值（这部分很好）
- 值更改后，我们想合并来自同一场景的栅格以寻找物种之间的相似性
我们一直在努力寻找一种不会产生重叠缓冲区或以更规范的方式聚合像素的方法。
我们的栅格分辨率为 30 弧秒（大约 1 平方公里），我们只能聚合像素，但我们需要特定区域，不可能只聚合整个像素生成值，例如 6 平方公里 另一种可能性是聚合部分像素以访问 6 平方公里等区域。可能吗？
到目前为止，我已经尝试了两种方法（提取和聚合）并且 none 返回了所需的结果。以下是目前使用的代码：
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, NA, NA, 0.721549, 0.6988058, 0.7333487, 0.8138089, 0.829727, 0.8689544, 0.8607966, 0.9794912, 0.9012381, NA, NA, NA, 0.7118917, 0.7103891, 0.7527786, 0.7792872, 0.8320968, 0.8082156, 0.920545, 0.9788723, 0.7859345, NA, NA, NA, 0.6824703, 0.7039984, 0.7589136, 0.7905939, 0.9024587, 0.9848859, 0.969553, 0.8404503, NA, NA, NA, 0.9922802, 0.6883243, 0.731391, 0.7682831, 0.8586601, 0.9850862, 0.9705435, 0.9888299, 0.8164479, 0.666971, NA, NA, 0.757661, 0.7628939, 0.7868114, 0.7910978, 0.7650495, 0.8227689, 0.8148086, 0.8691386, 0.7376462, 0.9176324, NA, NA, 0.998813, 0.7585487, 0.7721, 0.6508481, 0.6098195, 0.7708853, 0.7119401, 0.625409, 0.6886432, 0.8641906, NA, NA, 0.9991203, 0.7227083, 0.6550816, 0.5863016, 0.7050957, 0.6629267, 0.6550342, 0.6217013, 0.8762864, 0.9989462, NA, NA, NA, 0.6901366, 0.7041199, 0.6307223, 0.6305411, 0.7033732, 0.7092581, 0.7340803, 0.7865254, 0.9964261, 0.7940102, NA, NA, 0.6661926, 0.7381653, 0.6544684, 0.6170949, 0.641997, 0.7506128, 0.9248958, 0.9903375, 0.7662657, 0.7847621, NA, NA, 0.6582187, 0.7361452, 0.6488761, 0.6309077, 0.6542051, 0.781707, 0.940975, 0.9350743, 0.7824667, NA, NA, NA, 0.7215696, 0.7388574, 0.6753907, 0.6716958, 0.7162136, 0.7920918, 0.8702987, 0.9929227, 0.9775091, NA, NA, 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)
我们欢迎任何有关如何处理缓冲区的建议。
请在 code 中提供示例数据 --- 应该没有理由下载数据（请参阅 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)
如果您现在想为每个聚合缓冲区采样一个单元格
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.72, 0.88, 0.99, 0.99, 0.97, 0.82, 0.9, NA, NA, NA, 0.72, 0.7, 0.73, 0.81, 0.83, 0.87, 0.86, 0.98,0.9,NA,NA,NA,0.71,0.71,0.75,0.78,0.83,0.81,0.92,0.98,0.79,NA,NA,NA,0.68,0.7,0.76,0.79,0.9,0.98,0.97,0.84,NA,NA,NA,0.99,0.69,0.73,0.77,0.86,0.99,0.97,0.99,0.82,0.67,NA,NA,0.76,0.76,0.79,0.79,0.77,0.82,0.81,0.87,0.74,0.92,NA,NA,1,0.76,0.77,0.65,0.61,0.77,0.71,0.63,0.69,0.86,NA,NA,1,0.72,0.66,0.59,0.71,0.66,0.66,0.62,0.88,1,NA,NA,NA,0.69,0.7,0.63,0.63,0.7,0.71,0.73,0.79,1,0.79,NA,NA,0.67,0.74,0.65,0.62,0.64,0.75,0.92,0.99,0.77,0.78,NA,NA,0.66,0.74,0.65,0.63,0.65,0.78,0.94,0.94,0.78,NA,NA,NA,0.72,0.74,0.68,0.67,0.72,0.79,0.87,0.99,0.98,NA,NA,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, 1, NA)
p <- disaggregate(as.polygons(x))
b <- buffer(p, .01)
plot(x)
lines(b)