重塑 3D 数组并删除缺失值

Question

我有一个 NxMxT 数组，其中数组的每个元素都是地球网格。如果网格在海洋上空，则值为 999。如果网格在陆地上空，则它包含观测值。 N是经度，M是纬度，T是月。 特别是，我有一个名为 tmp60 的数组，用于从 1960 年到 1969 年的十年，因此每个网格有 120 个月。

为了测试 1960 年 1 月的全球平均值，我写道：

tmpJan60=tmp60(:,:,1);
tmpJan60(tmpJan60(:,:)>200)=NaN;
nanmean(nanmean(tmpJan60))

这给了我 5.855。

我对 reshape 函数感到困惑。我认为以下代码应该产生相同的平均值，即 5.855，但事实并非如此：

load tmp60

N1=size(tmp60,1)
N2=size(tmp60,2)
N3=size(tmp60,3)

reshtmp60 = reshape(tmp60, N1*N2,N3);
reshtmp60( reshtmp60(:,1)>200,: )=[];
mean(reshtmp60(:,1))

这给了我-1.6265，这是不正确的。

我已经在 Excel (!) 中检查了结果，并且 5.855 是正确的，所以我认为我在 reshape 函数中犯了一个错误。

理想情况下，我想要一个包含每个网格的矩阵，首先沿着 N 维向下移动，然后生成 720 行和 120 列（每列是一个月）。这些前 720 行将代表地球周围相同纬度的一个经度带。接下来，我想将纬度增加 1，从而增加 720 行和 120 列。最终，我想为所有 360 纬度做到这一点。 如果输入经度和纬度，例如第 1 列和第 2 列，则矩阵应如下所示：

temp = [-179.75 -89.75 -1  2 ...
        -179.25 -89.75 2  4 ...
           ...
         179.75 -89.75 5 9 ...
        -179.75 -89.25 2  5 ...
        -179.25 -89.25 3  4 ... 
          ...
        -179.75  89.75 2  3 ...
          ...
         179.75  89.75 6  9 ...]

所以temp(:,3) 应该是 1960 年 1 月的所有观测值。

一种方法是：

grid1 = tmp60(1,1,:);
g1 = reshape(grid1, [1,120]);
grid2 = tmp60(2,1,:);
g2 = reshape(grid2,[1,120]);
g = [g1;g2];

但显然很麻烦。

我无法为 N*M 元素自动执行此过程，因此感谢 cmets！

A link to the file tmp60.mat

Answer 1

您的代码中的主要问题是处理nans。请注意以下示例：

a = randi(10,6);
a(a>7)=nan
m = [mean(a(:),'omitnan') mean(mean(a,'omitnan'),'omitnan')]

m =
       3.8421       3.6806

m 中的两个元素只是a 中所有元素的平均值。但他们不一样！原因是取所有值的平均值，mean(a(:),'omitnan') 就像对所有非 nan 值求和，然后除以我们求和的值的数量：

sum(a(:),'omitnan')/sum(~isnan(a(:)))==mean(a(:),'omitnan') % this is true

但取第一个维度的平均值，我们得到 6 个平均值：

sum(a,'omitnan')./sum(~isnan(a))==mean(a,'omitnan') % this is also true

当我们取它们的平均值时，我们除以一个更大的数字，因为所有nans 都已经被省略了：

mean(sum(a,'omitnan')./sum(~isnan(a)))==mean(a(:),'omitnan') % this is false

---

这是我认为你想要的代码：

% this is exactly as your first test:
tmpJan60=tmn60(:,:,1);
tmpJan60(tmpJan60>200) = nan;
m1 = mean(mean(tmpJan60,'omitnan'),'omitnan')

% this creates the matrix as you want it:
result = reshape(permute(tmn60,[3 1 2]),120,[]).';
result(result>200) = nan;
r = reshape(result(:,1),720,360);
m2 = mean(mean(r,'omitnan'),'omitnan')

isequal(m1,m2)

要创建矩阵，您首先要置换维度，以便您想要保持原样（时间）的那个将是第一个。然后将数组重塑为 Tx(lon*lat)，因此所有时间步长为 120 行，所有坐标组合为 259200 列。剩下的就是转置它。

m1 是您的第一个计算，m2 是您尝试在第二个计算中执行的操作。它们在这里是相等的，但它们的值不是 5.855，即使我使用你的代码。

但是，我认为正确的解决方案是将所有值的平均值放在一起：

mean(result(:,1),'omitnan')