如何使用 R 在 3D 数组（沿维度时间）中找到周期性出现的 NA 值
How can I find periodically appearing NA values in an 3D array (along dimension time) with R
我在一个数组中有空间数据（最初是 ncdf）的时间序列（几年内的月度值）。如果连续出现超过 2 个，例如带有 NA 的 1 月，我想通过在所有时间步长中将其放入 NA 来完全禁止进一步研究这个像素（现在是一个时间步长矩阵中的单元格）。
就我而言，"time.series" 仅对向量或矩阵（二维的最大值）有效。
我可以看到（但也无法实施）的一个解决方法是：
以不再纯粹按时间顺序排列而是按月排序的方式对数组进行排序（2001 年 1 月、2002 年 1 月、2003 年 1 月、2001 年 2 月、2002 年 2 月、2003 年 2 月，...）已经有很大帮助了。但是，如果例如，它会留下像素得到 NA 的情况。 2002 年 1 月、2003 年 1 月和 2001 年 2 月是 NA。
任何帮助将不胜感激。请询问我的问题是否不清楚 - 这是我的第一个问题 - 我已尽力而为。
编辑：
我的实际数据集是基于全球卫星的辐射数据集。由于例如周期性出现的云（在每年同月的雨季期间），不应进一步考虑这些像素。我还有一些其他标准可以消除像素。只缺少一个标准。
# create any array with scattered NAs
set.seed (10)
array <- replicate(48, replicate(10, rnorm(20)))
na_pixels <- array((sample(c(1, NA), size = 7200, replace = TRUE, prob = c(0.95, 0.05))), dim = c(20,10,48))
na_array <- array * na_pixels
dimnames(na_array) <- list(NULL, NULL, as.character(seq(as.Date("2001-01-01"), as.Date("2004-12-01"), "month")))
#I want to test several conditions that would make a pixel not usable for me
#in the end I want to retrieve a mask of usable "pixels".
#what I am doing already is:
mask <- apply(na_array, MARGIN = c(1,2), FUN=function(x){
#check if more than 10% of a pixel are NA over time
if (sum(is.na(x)) > (length(x)*0.05)){
mask_val <- 0
}
#check if more than 5 pixel are missing consecutively
else if (max(with(rle(is.na(a)), lengths[values])) > 5){
mask_val <- 0
}
#this is the missing part
else if (...more than 2 januaries or 2 feburaries or... are NA){#check for periodically appearing NAs
mask_val <- 0
}
else {
mask_val <- 1
}
return(mask_val)
})
在 'long' "data.frame":
中更改 3D 数组可能更方便（如果存在必要的内存）
as.data.frame(as.table(na_array))
# Var1 Var2 Var3 Freq
#1 A A 2001-01-01 0.01874617
#2 B A 2001-01-01 -0.18425254
#3 C A 2001-01-01 -1.37133055
# ...........................
#9598 R J 2004-12-01 NA
#9599 S J 2004-12-01 -1.11411416
#9600 T J 2004-12-01 0.01435433
与其依赖 as.table 和 as.data.frame 强制转换，不如手动完成，效率更高：
dat = data.frame(i = rep_len(seq_len(dim(na_array)[1]), prod(dim(na_array))),
j = rep_len(rep(seq_len(dim(na_array)[2]), each = dim(na_array)[1]), prod(dim(na_array))),
date = rep(as.Date(dimnames(na_array)[[3]]), each = prod(dim(na_array)[1:2])) ,
month = rep(format(as.Date(dimnames(na_array)[[3]]), "%b"), each = prod(dim(na_array)[1:2])),
isNA = c(is.na(na_array)))
dat
# i j date month isNA
#1 1 1 2001-01-01 Jan FALSE
#2 2 1 2001-01-01 Jan FALSE
#3 3 1 2001-01-01 Jan FALSE
#4 4 1 2001-01-01 Jan TRUE
# ..............
#9597 17 10 2004-12-01 Dec FALSE
#9598 18 10 2004-12-01 Dec TRUE
#9599 19 10 2004-12-01 Dec FALSE
#9600 20 10 2004-12-01 Dec FALSE
其中 i：na_array 中的行，j：na_array 中的列，date：na_array 的第 3 个维度，month：date列的月份（后面会用到），isNA：na_array的值是否为NA。
并建立三个条件：
cond1 = aggregate(isNA ~ i + j, dat, function(x) sum(x) > (dim(na_array)[3] * 0.05))
（创建 cond1 更有效的方法是 rowSums(is.na(na_array), dims = 2) > (dim(na_array)[3] * 0.05)）。
cond2 = aggregate(isNA ~ i + j, dat, function(x) any(with(rle(x), lengths[values]) > 5))
并计算 cond3，首先找到每个 "month" 每个 'cell'（即 [i, j]）的缺失值数量（"month" 是一个在开始创建 'long' "data.frame" dat 时来自 dimnames(na_array)[[3]] 的变量 created/extracted:
NA_per_month = aggregate(isNA ~ i + j + month, dat, function(x) sum(x))
对于每个 [i, j]，每个 "month" 有 NAs 的数量，我们通过检查每个 [i, j] 是否包含 any 来构建 cond3 ] "month" 超过 2 NAs:
cond3 = aggregate(isNA ~ i + j, NA_per_month, function(x) any(x > 2))
（将上述 'group-by' 操作中的 aggregate 替换为任何其他可用的操作很简单）。
也许我们可以避免创建 'long' "data.frame" 并直接对 na_array 进行操作。例如，使用 rowSums 版本计算 cond1 更加高效和直接。 cond2 也可以通过 na_array 上的 apply 保存。但是 cond3 使用 'long' "data.frame" 比使用 3D 数组要简单得多。因此，考虑到效率，尝试使用数据中存在的结构总是更好，如果它变得足够麻烦，那么我们可能应该更改一次数据结构并在另一个脚手架中计算任何东西。
要得到最终结果，分配一个合适大小的"matrix"：
ans = matrix(NA, dim(na_array)[1], dim(na_array)[2])
并在OR条件后填写：
ans[cbind(cond1$i, cond1$j)] = cond1$isNA | cond2$isNA | cond3$isNA
ans
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
# [1,] TRUE TRUE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE
# [2,] TRUE FALSE FALSE FALSE TRUE TRUE FALSE TRUE FALSE FALSE
# [3,] FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE
# [4,] FALSE FALSE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE
# [5,] FALSE FALSE TRUE FALSE TRUE FALSE FALSE TRUE FALSE FALSE
# [6,] FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
# [7,] FALSE FALSE TRUE TRUE TRUE FALSE FALSE TRUE TRUE FALSE
# [8,] TRUE TRUE TRUE TRUE FALSE FALSE TRUE FALSE TRUE FALSE
# [9,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
#[10,] TRUE FALSE TRUE TRUE FALSE FALSE FALSE TRUE FALSE FALSE
#[11,] FALSE TRUE TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE
#[12,] TRUE TRUE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
#[13,] FALSE TRUE TRUE FALSE TRUE FALSE FALSE TRUE FALSE TRUE
#[14,] FALSE FALSE TRUE FALSE TRUE FALSE FALSE TRUE FALSE TRUE
#[15,] TRUE TRUE TRUE TRUE FALSE TRUE FALSE FALSE TRUE FALSE
#[16,] FALSE FALSE FALSE TRUE TRUE FALSE TRUE TRUE TRUE TRUE
#[17,] TRUE FALSE TRUE TRUE FALSE FALSE TRUE FALSE TRUE FALSE
#[18,] FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE TRUE TRUE
#[19,] FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE FALSE
#[20,] TRUE FALSE TRUE TRUE FALSE TRUE TRUE FALSE FALSE TRUE
@ alexis_laz：是的，现在可以了。不幸的是，我意识到 ans[cbind(cond1$i, cond1$j)] = cond1$isNA | cond2$isNA | cond3$isNA 不起作用。我收到错误消息：要替换的项目数不是替换长度的倍数。我认为它只需要 cond1 进行更换。 （我很抱歉我的示例数据集在所有情况下都为 cond2 和 cond3 提供 'FALSE'，但它仍然应该检查 code.Even 中的 'OR'，尽管结果看起来与 cond1 相同) 我想出了以下代码，它可以工作，但肯定不是很好或高效，因为我不太熟悉布尔值。也许你可以优化我的代码或编辑你的行（因为我的真实数据集很大，我会很高兴能进行任何优化）。在远端，我需要所有 True 条件（即 NA）为 0，所有 FALSE 条件为 1。这就是为什么我已经在我的代码中这样做了。
ans = matrix(NA, dim(na_array)[1], dim(na_array)[2])
cond1_bool <- ans
cond1_bool[cbind(cond1$i, cond1$j)] = cond1$isNA
cond2_bool <- ans
cond2_bool[cbind(cond2$i, cond2$j)] = cond2$isNA
cond3_bool <- ans
cond3_bool[cbind(cond3$i, cond3$j)] = cond3$isNA
ans_bool <- ans
ans_bool[which(cond1_bool == T|cond2_bool == T|cond3_bool == T)] <- 0
ans_bool[which(is.na(ans_bool))] <- 1
我在一个数组中有空间数据（最初是 ncdf）的时间序列（几年内的月度值）。如果连续出现超过 2 个，例如带有 NA 的 1 月，我想通过在所有时间步长中将其放入 NA 来完全禁止进一步研究这个像素（现在是一个时间步长矩阵中的单元格）。
就我而言，"time.series" 仅对向量或矩阵（二维的最大值）有效。
我可以看到（但也无法实施）的一个解决方法是： 以不再纯粹按时间顺序排列而是按月排序的方式对数组进行排序（2001 年 1 月、2002 年 1 月、2003 年 1 月、2001 年 2 月、2002 年 2 月、2003 年 2 月，...）已经有很大帮助了。但是，如果例如，它会留下像素得到 NA 的情况。 2002 年 1 月、2003 年 1 月和 2001 年 2 月是 NA。
任何帮助将不胜感激。请询问我的问题是否不清楚 - 这是我的第一个问题 - 我已尽力而为。
编辑： 我的实际数据集是基于全球卫星的辐射数据集。由于例如周期性出现的云（在每年同月的雨季期间），不应进一步考虑这些像素。我还有一些其他标准可以消除像素。只缺少一个标准。
# create any array with scattered NAs
set.seed (10)
array <- replicate(48, replicate(10, rnorm(20)))
na_pixels <- array((sample(c(1, NA), size = 7200, replace = TRUE, prob = c(0.95, 0.05))), dim = c(20,10,48))
na_array <- array * na_pixels
dimnames(na_array) <- list(NULL, NULL, as.character(seq(as.Date("2001-01-01"), as.Date("2004-12-01"), "month")))
#I want to test several conditions that would make a pixel not usable for me
#in the end I want to retrieve a mask of usable "pixels".
#what I am doing already is:
mask <- apply(na_array, MARGIN = c(1,2), FUN=function(x){
#check if more than 10% of a pixel are NA over time
if (sum(is.na(x)) > (length(x)*0.05)){
mask_val <- 0
}
#check if more than 5 pixel are missing consecutively
else if (max(with(rle(is.na(a)), lengths[values])) > 5){
mask_val <- 0
}
#this is the missing part
else if (...more than 2 januaries or 2 feburaries or... are NA){#check for periodically appearing NAs
mask_val <- 0
}
else {
mask_val <- 1
}
return(mask_val)
})
在 'long' "data.frame":
中更改 3D 数组可能更方便（如果存在必要的内存）as.data.frame(as.table(na_array))
# Var1 Var2 Var3 Freq
#1 A A 2001-01-01 0.01874617
#2 B A 2001-01-01 -0.18425254
#3 C A 2001-01-01 -1.37133055
# ...........................
#9598 R J 2004-12-01 NA
#9599 S J 2004-12-01 -1.11411416
#9600 T J 2004-12-01 0.01435433
与其依赖 as.table 和 as.data.frame 强制转换，不如手动完成，效率更高：
dat = data.frame(i = rep_len(seq_len(dim(na_array)[1]), prod(dim(na_array))),
j = rep_len(rep(seq_len(dim(na_array)[2]), each = dim(na_array)[1]), prod(dim(na_array))),
date = rep(as.Date(dimnames(na_array)[[3]]), each = prod(dim(na_array)[1:2])) ,
month = rep(format(as.Date(dimnames(na_array)[[3]]), "%b"), each = prod(dim(na_array)[1:2])),
isNA = c(is.na(na_array)))
dat
# i j date month isNA
#1 1 1 2001-01-01 Jan FALSE
#2 2 1 2001-01-01 Jan FALSE
#3 3 1 2001-01-01 Jan FALSE
#4 4 1 2001-01-01 Jan TRUE
# ..............
#9597 17 10 2004-12-01 Dec FALSE
#9598 18 10 2004-12-01 Dec TRUE
#9599 19 10 2004-12-01 Dec FALSE
#9600 20 10 2004-12-01 Dec FALSE
其中 i：na_array 中的行，j：na_array 中的列，date：na_array 的第 3 个维度，month：date列的月份（后面会用到），isNA：na_array的值是否为NA。
并建立三个条件：
cond1 = aggregate(isNA ~ i + j, dat, function(x) sum(x) > (dim(na_array)[3] * 0.05))
（创建 cond1 更有效的方法是 rowSums(is.na(na_array), dims = 2) > (dim(na_array)[3] * 0.05)）。
cond2 = aggregate(isNA ~ i + j, dat, function(x) any(with(rle(x), lengths[values]) > 5))
并计算 cond3，首先找到每个 "month" 每个 'cell'（即 [i, j]）的缺失值数量（"month" 是一个在开始创建 'long' "data.frame" dat 时来自 dimnames(na_array)[[3]] 的变量 created/extracted:
NA_per_month = aggregate(isNA ~ i + j + month, dat, function(x) sum(x))
对于每个 [i, j]，每个 "month" 有 NAs 的数量，我们通过检查每个 [i, j] 是否包含 any 来构建 cond3 ] "month" 超过 2 NAs:
cond3 = aggregate(isNA ~ i + j, NA_per_month, function(x) any(x > 2))
（将上述 'group-by' 操作中的 aggregate 替换为任何其他可用的操作很简单）。
也许我们可以避免创建 'long' "data.frame" 并直接对 na_array 进行操作。例如，使用 rowSums 版本计算 cond1 更加高效和直接。 cond2 也可以通过 na_array 上的 apply 保存。但是 cond3 使用 'long' "data.frame" 比使用 3D 数组要简单得多。因此，考虑到效率，尝试使用数据中存在的结构总是更好，如果它变得足够麻烦，那么我们可能应该更改一次数据结构并在另一个脚手架中计算任何东西。
要得到最终结果，分配一个合适大小的"matrix"：
ans = matrix(NA, dim(na_array)[1], dim(na_array)[2])
并在OR条件后填写：
ans[cbind(cond1$i, cond1$j)] = cond1$isNA | cond2$isNA | cond3$isNA
ans
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
# [1,] TRUE TRUE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE
# [2,] TRUE FALSE FALSE FALSE TRUE TRUE FALSE TRUE FALSE FALSE
# [3,] FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE
# [4,] FALSE FALSE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE
# [5,] FALSE FALSE TRUE FALSE TRUE FALSE FALSE TRUE FALSE FALSE
# [6,] FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
# [7,] FALSE FALSE TRUE TRUE TRUE FALSE FALSE TRUE TRUE FALSE
# [8,] TRUE TRUE TRUE TRUE FALSE FALSE TRUE FALSE TRUE FALSE
# [9,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
#[10,] TRUE FALSE TRUE TRUE FALSE FALSE FALSE TRUE FALSE FALSE
#[11,] FALSE TRUE TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE
#[12,] TRUE TRUE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
#[13,] FALSE TRUE TRUE FALSE TRUE FALSE FALSE TRUE FALSE TRUE
#[14,] FALSE FALSE TRUE FALSE TRUE FALSE FALSE TRUE FALSE TRUE
#[15,] TRUE TRUE TRUE TRUE FALSE TRUE FALSE FALSE TRUE FALSE
#[16,] FALSE FALSE FALSE TRUE TRUE FALSE TRUE TRUE TRUE TRUE
#[17,] TRUE FALSE TRUE TRUE FALSE FALSE TRUE FALSE TRUE FALSE
#[18,] FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE TRUE TRUE
#[19,] FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE FALSE
#[20,] TRUE FALSE TRUE TRUE FALSE TRUE TRUE FALSE FALSE TRUE
@ alexis_laz：是的，现在可以了。不幸的是，我意识到 ans[cbind(cond1$i, cond1$j)] = cond1$isNA | cond2$isNA | cond3$isNA 不起作用。我收到错误消息：要替换的项目数不是替换长度的倍数。我认为它只需要 cond1 进行更换。 （我很抱歉我的示例数据集在所有情况下都为 cond2 和 cond3 提供 'FALSE'，但它仍然应该检查 code.Even 中的 'OR'，尽管结果看起来与 cond1 相同) 我想出了以下代码，它可以工作，但肯定不是很好或高效，因为我不太熟悉布尔值。也许你可以优化我的代码或编辑你的行（因为我的真实数据集很大，我会很高兴能进行任何优化）。在远端，我需要所有 True 条件（即 NA）为 0，所有 FALSE 条件为 1。这就是为什么我已经在我的代码中这样做了。
ans = matrix(NA, dim(na_array)[1], dim(na_array)[2])
cond1_bool <- ans
cond1_bool[cbind(cond1$i, cond1$j)] = cond1$isNA
cond2_bool <- ans
cond2_bool[cbind(cond2$i, cond2$j)] = cond2$isNA
cond3_bool <- ans
cond3_bool[cbind(cond3$i, cond3$j)] = cond3$isNA
ans_bool <- ans
ans_bool[which(cond1_bool == T|cond2_bool == T|cond3_bool == T)] <- 0
ans_bool[which(is.na(ans_bool))] <- 1