重塑 3D 数组并删除缺失值
Reshape a 3D array and remove missing values
我有一个 NxMxT 数组,其中数组的每个元素都是一个地球网格。如果网格在海洋上方,则该值为 999。如果网格在陆地上方,则它包含一个观测值。 N 是经度,M 是纬度,T 是月份。
特别是,我有一个名为 tmp60
的数组,表示从 1960 年到 1969 年的十年,因此每个网格有 120 个月。
为了测试 1960 年 1 月的全球平均水平,我写下:
tmpJan60=tmp60(:,:,1);
tmpJan60(tmpJan60(:,:)>200)=NaN;
nanmean(nanmean(tmpJan60))
这给了我 5.855。
我对重塑函数感到困惑。我认为下面的代码应该产生相同的平均值,即 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 个元素自动执行此过程,因此不胜感激!
您代码中的主要问题是处理 nan
。观察以下示例:
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
当我们取它们的平均值时,我们除以一个更大的数,因为所有 nan
s 都已经被省略了:
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')
我有一个 NxMxT 数组,其中数组的每个元素都是一个地球网格。如果网格在海洋上方,则该值为 999。如果网格在陆地上方,则它包含一个观测值。 N 是经度,M 是纬度,T 是月份。
特别是,我有一个名为 tmp60
的数组,表示从 1960 年到 1969 年的十年,因此每个网格有 120 个月。
为了测试 1960 年 1 月的全球平均水平,我写下:
tmpJan60=tmp60(:,:,1);
tmpJan60(tmpJan60(:,:)>200)=NaN;
nanmean(nanmean(tmpJan60))
这给了我 5.855。
我对重塑函数感到困惑。我认为下面的代码应该产生相同的平均值,即 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 个元素自动执行此过程,因此不胜感激!
您代码中的主要问题是处理 nan
。观察以下示例:
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
当我们取它们的平均值时,我们除以一个更大的数,因为所有 nan
s 都已经被省略了:
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')