如何在带有光栅的列表中循环内部循环？
How to loop inside loop in list with raster?
我有下面的空间数据集。它由 space 中的点组成，这些点构成每个 ID（例如个体动物）和 SUB-IDS（个体动物的路线模拟）所做的轨迹。
library(tidyverse)
library(sf)
library(raster)
library(data.table)
ids <- c("844", "844", "844", "844", "844","844", "844", "844", "844", "844",
"844", "844", "845", "845", "845", "845", "845","845", "845", "845",
"845","845", "845", "845")
sub_ids <- c("2017_844_1", "2017_844_1", "2017_844_1", "2017_844_1", "2017_844_2",
"2017_844_2", "2017_844_2", "2017_844_2", "2017_844_3", "2017_844_3",
"2017_844_3", "2017_844_3", "2017_845_1", "2017_845_1", "2017_845_1",
"2017_845_1", "2017_845_2","2017_845_2", "2017_845_2", "2017_845_2",
"2017_845_3","2017_845_3", "2017_845_3", "2017_845_3")
lat <- c(-30.6456, -29.5648, -27.6667, -31.5587, -30.6934, -29.3147, -23.0538,
-26.5877, -26.6923, -23.40865, -23.1143, -23.28331, -31.6456, -24.5648,
-27.6867, -31.4587, -30.6784, -28.3447, -23.0466, -27.5877, -26.8524,
-23.8855, -24.1143, -23.5874)
long <- c(-50.4879, -49.8715, -51.8716, -50.4456, -50.9842, -51.9787, -41.2343,
-40.2859, -40.19599, -41.64302, -41.58042, -41.55057, -50.4576, -48.8715,
-51.4566, -51.4456, -50.4477, -50.9937, -41.4789, -41.3859, -40.2536,
-41.6502, -40.5442, -41.4057)
df <- tibble(ids, sub_ids, lat, long)
#converting to sf
df.sf <- df %>%
sf::st_as_sf(coords = c("long", "lat"), crs = 4326)
我还有一个由点创建的 sf 网格：
#creating grid
xy <- sf::st_coordinates(df.sf)
grid = sf::st_make_grid(sf::st_bbox(df.sf),
cellsize = .1, square = FALSE) %>%
sf::st_as_sf() %>%
dplyr::mutate(cell = 1:nrow(.))
loop1 inside mylist[[j]]: 这个想法是为每个 SUB ID 生成一个包含每个单元格点数的栅格图层，并将这些图层添加到 rasterbrick 对象中，计算rasterbrick 平均值并产生 1 个 rasterlayer .
loop2 inside mylist: 在每个列表元素中执行此过程，为每个 ID 创建 1 个栅格图层，然后再次将生成的栅格图层添加到栅格砖中以计算平均值并得到单个栅格图层。
#creating rasters objects
r <- raster::raster(grid, res=0.1)
rr <- raster::brick(r)
rr2 <- raster::brick(r)
df.sf$ids <- as.factor(df.sf$ids)
mylist <- split(df.sf, df.sf$ids) #list to each ID
#looping
for(j in 1:length(unique(mylist))){
for(i in 1:length(unique(mylist[[j]]$sub_ids))) {
output1[j]<- stackApply(addLayer(rr,rasterize(subset(mylist[[j]],
sub_ids == unique(mylist[[j]]$sub_ids)[i]),
y = r,
field = 1,
background = 0,
fun="count")), indices = rep(1,nlayers(rr)),
fun = "mean", na.rm = T)+1
output2 <- stackApply(addLayer(rr2, output1[j]), indices = rep(1,nlayers(rr2)),
fun = "mean", na.rm = T) }}
我尝试使用 for() 如上所示，但出现此错误：
Error in h(simpleError(msg, call)) :
erro na avaliação do argumento 'x' na seleção do método para a função 'nlayers': 'operations are possible only for numeric, logical or complex types'
In addition: Warning message:
In x@data@values[i] <- value :
number of items to replace is not a multiple of replacement length
编辑：每个子集由一个 ID 处理
############## ID 844
#subset one ID
id.844 <- subset(df.sf,df.sf$ids=="844")
id.844$sub_ids <- as.factor(id.844$sub_ids)
#subset SUB IDS
id.844.1 <- subset(id.844, id.844$sub_ids=="2017_844_1")
id.844.2 <- subset(id.844, id.844$sub_ids=="2017_844_2")
id.844.3 <- subset(id.844, id.844$sub_ids=="2017_844_3")
#SUB ID 1
output1 <- raster::rasterize(x = id.844.1,
y = r,
field = 1,
background = 0,
fun="count")
#SUB ID 2
output2 <- raster::rasterize(x = id.844.2,
y = r,
field = 1,
background = 0,
fun="count")
#SUB ID 3
output3 <- raster::rasterize(x = id.844.3,
y = r,
field = 1,
background = 0,
fun="count")
#joining rasters SUB IDS
rasterbrick<- raster::brick(output1, output2, output3)
#mean ID 1
result.id1 <- raster::stackApply(rasterbrick,
indices = rep(1,nlayers(rasterbrick)),
fun = "mean", na.rm = T)
考虑使用 user-defined 函数概括您的流程，然后传递给带有嵌套 by 的方法子集，将 object_oriented 包装器传递给 tapply数据框，然后将子集传递给一个函数，在很大程度上等同于 split+lapply。使用 do.call 组合光栅砖。
build_raster_mean <- function(sub_df) {
rasterize_list <- by(
sub_df,
sub_df$sub_ids,
function(sub_id_df) raster::rasterize(
x = sub_id_df, y = r, field = 1, background = 0, fun="count"
),
simplify = FALSE
) |> `attributes<-`(NULL)
raster_brick <- do.call(raster::brick, rasterize_list)
raster_mean <- raster::stackApply(
raster_brick,
indices = rep(1,nlayers(raster_brick)),
fun = "mean",
na.rm = TRUE
)
}
raster_mean_list <- by(
df.sf, df.sf$ids, build_raster_mean, simplify=FALSE
)
最终解决方案
build_raster_mean <- function(sub_df) {
rasterize_list <- by(
sub_df,
sub_df$sub_ids,
function(sub_id_df) raster::rasterize(
x = sub_id_df, y = r, field=1, background = 0, fun='count'
)
)
rasterize_list2 <- unlist(rasterize_list, recursive=F)
raster_brick <- raster::writeRaster(brick(rasterize_list2),
names(rasterize_list2), bylayer=TRUE,
overwrite=TRUE)
raster_mean <- raster::stackApply(
raster_brick,
indices = rep(1,nlayers(raster_brick)),
fun = "mean",
na.rm = TRUE
)
}
#loop in list
raster_mean_list <- purrr::map(mylist, build_raster_mean)
#new raster brick with IDS mean
raster_mean_brick <- raster::writeRaster(brick(raster_mean_list),
names(raster_mean_list), bylayer=TRUE,
overwrite=TRUE)
#rasterlayer with average final
raster_mean_final <- raster::stackApply(raster_mean_brick ,
indices = rep(1,nlayers(raster_mean_brick)),
fun = "mean",
na.rm = TRUE)
我有下面的空间数据集。它由 space 中的点组成，这些点构成每个 ID（例如个体动物）和 SUB-IDS（个体动物的路线模拟）所做的轨迹。
library(tidyverse)
library(sf)
library(raster)
library(data.table)
ids <- c("844", "844", "844", "844", "844","844", "844", "844", "844", "844",
"844", "844", "845", "845", "845", "845", "845","845", "845", "845",
"845","845", "845", "845")
sub_ids <- c("2017_844_1", "2017_844_1", "2017_844_1", "2017_844_1", "2017_844_2",
"2017_844_2", "2017_844_2", "2017_844_2", "2017_844_3", "2017_844_3",
"2017_844_3", "2017_844_3", "2017_845_1", "2017_845_1", "2017_845_1",
"2017_845_1", "2017_845_2","2017_845_2", "2017_845_2", "2017_845_2",
"2017_845_3","2017_845_3", "2017_845_3", "2017_845_3")
lat <- c(-30.6456, -29.5648, -27.6667, -31.5587, -30.6934, -29.3147, -23.0538,
-26.5877, -26.6923, -23.40865, -23.1143, -23.28331, -31.6456, -24.5648,
-27.6867, -31.4587, -30.6784, -28.3447, -23.0466, -27.5877, -26.8524,
-23.8855, -24.1143, -23.5874)
long <- c(-50.4879, -49.8715, -51.8716, -50.4456, -50.9842, -51.9787, -41.2343,
-40.2859, -40.19599, -41.64302, -41.58042, -41.55057, -50.4576, -48.8715,
-51.4566, -51.4456, -50.4477, -50.9937, -41.4789, -41.3859, -40.2536,
-41.6502, -40.5442, -41.4057)
df <- tibble(ids, sub_ids, lat, long)
#converting to sf
df.sf <- df %>%
sf::st_as_sf(coords = c("long", "lat"), crs = 4326)
我还有一个由点创建的 sf 网格：
#creating grid
xy <- sf::st_coordinates(df.sf)
grid = sf::st_make_grid(sf::st_bbox(df.sf),
cellsize = .1, square = FALSE) %>%
sf::st_as_sf() %>%
dplyr::mutate(cell = 1:nrow(.))
loop1 inside mylist[[j]]: 这个想法是为每个 SUB ID 生成一个包含每个单元格点数的栅格图层，并将这些图层添加到 rasterbrick 对象中，计算rasterbrick 平均值并产生 1 个 rasterlayer .
loop2 inside mylist: 在每个列表元素中执行此过程，为每个 ID 创建 1 个栅格图层，然后再次将生成的栅格图层添加到栅格砖中以计算平均值并得到单个栅格图层。
#creating rasters objects
r <- raster::raster(grid, res=0.1)
rr <- raster::brick(r)
rr2 <- raster::brick(r)
df.sf$ids <- as.factor(df.sf$ids)
mylist <- split(df.sf, df.sf$ids) #list to each ID
#looping
for(j in 1:length(unique(mylist))){
for(i in 1:length(unique(mylist[[j]]$sub_ids))) {
output1[j]<- stackApply(addLayer(rr,rasterize(subset(mylist[[j]],
sub_ids == unique(mylist[[j]]$sub_ids)[i]),
y = r,
field = 1,
background = 0,
fun="count")), indices = rep(1,nlayers(rr)),
fun = "mean", na.rm = T)+1
output2 <- stackApply(addLayer(rr2, output1[j]), indices = rep(1,nlayers(rr2)),
fun = "mean", na.rm = T) }}
我尝试使用 for() 如上所示，但出现此错误：
Error in h(simpleError(msg, call)) :
erro na avaliação do argumento 'x' na seleção do método para a função 'nlayers': 'operations are possible only for numeric, logical or complex types'
In addition: Warning message:
In x@data@values[i] <- value :
number of items to replace is not a multiple of replacement length
编辑：每个子集由一个 ID 处理
############## ID 844
#subset one ID
id.844 <- subset(df.sf,df.sf$ids=="844")
id.844$sub_ids <- as.factor(id.844$sub_ids)
#subset SUB IDS
id.844.1 <- subset(id.844, id.844$sub_ids=="2017_844_1")
id.844.2 <- subset(id.844, id.844$sub_ids=="2017_844_2")
id.844.3 <- subset(id.844, id.844$sub_ids=="2017_844_3")
#SUB ID 1
output1 <- raster::rasterize(x = id.844.1,
y = r,
field = 1,
background = 0,
fun="count")
#SUB ID 2
output2 <- raster::rasterize(x = id.844.2,
y = r,
field = 1,
background = 0,
fun="count")
#SUB ID 3
output3 <- raster::rasterize(x = id.844.3,
y = r,
field = 1,
background = 0,
fun="count")
#joining rasters SUB IDS
rasterbrick<- raster::brick(output1, output2, output3)
#mean ID 1
result.id1 <- raster::stackApply(rasterbrick,
indices = rep(1,nlayers(rasterbrick)),
fun = "mean", na.rm = T)
考虑使用 user-defined 函数概括您的流程，然后传递给带有嵌套 by 的方法子集，将 object_oriented 包装器传递给 tapply数据框，然后将子集传递给一个函数，在很大程度上等同于 split+lapply。使用 do.call 组合光栅砖。
build_raster_mean <- function(sub_df) {
rasterize_list <- by(
sub_df,
sub_df$sub_ids,
function(sub_id_df) raster::rasterize(
x = sub_id_df, y = r, field = 1, background = 0, fun="count"
),
simplify = FALSE
) |> `attributes<-`(NULL)
raster_brick <- do.call(raster::brick, rasterize_list)
raster_mean <- raster::stackApply(
raster_brick,
indices = rep(1,nlayers(raster_brick)),
fun = "mean",
na.rm = TRUE
)
}
raster_mean_list <- by(
df.sf, df.sf$ids, build_raster_mean, simplify=FALSE
)
最终解决方案
build_raster_mean <- function(sub_df) {
rasterize_list <- by(
sub_df,
sub_df$sub_ids,
function(sub_id_df) raster::rasterize(
x = sub_id_df, y = r, field=1, background = 0, fun='count'
)
)
rasterize_list2 <- unlist(rasterize_list, recursive=F)
raster_brick <- raster::writeRaster(brick(rasterize_list2),
names(rasterize_list2), bylayer=TRUE,
overwrite=TRUE)
raster_mean <- raster::stackApply(
raster_brick,
indices = rep(1,nlayers(raster_brick)),
fun = "mean",
na.rm = TRUE
)
}
#loop in list
raster_mean_list <- purrr::map(mylist, build_raster_mean)
#new raster brick with IDS mean
raster_mean_brick <- raster::writeRaster(brick(raster_mean_list),
names(raster_mean_list), bylayer=TRUE,
overwrite=TRUE)
#rasterlayer with average final
raster_mean_final <- raster::stackApply(raster_mean_brick ,
indices = rep(1,nlayers(raster_mean_brick)),
fun = "mean",
na.rm = TRUE)