Reshape a 3D array and remove missing values

Question

I have an NxMxT array where each element of the array is a grid of Earth. If the grid is over the ocean, then the value is 999. If the grid is over land, it contains an observed value. N is longitude, M is latitude, and T is months. In particular, I have an array called tmp60 for the ten years 1960 through 1969, so 120 months for each grid.

To test what the global mean in January 1960 was, I write:

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

which gives me 5.855.

I am confused about the reshape function. I thought the following code should yield the same average, namely 5.855, but it does not:

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))

this gives me -1.6265, which is not correct.

I have checked the result in Excel (!) and 5.855 is correct, so I assume I make a mistake in the reshape function.

Ideally, I want a matrix that takes each grid, going first down the N-dimension, and make the 720 rows with 120 columns (each column is a month). These first 720 rows will represent one longitude band around Earth for the same latitude. Next, I want to increase the latitude by 1, thus another 720 rows with 120 columns. Ultimately I want to do this for all 360 latitudes. If longitude and latitude were inputs, say column 1 and 2, then the matrix should look like this:

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 ...]

So temp(:,3) should be all January 1960 observations.

One way to do this is:

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

But obviously very cumbersome.

I am not able to automate this procedure for the N*M elements, so comments are appreciated!

A link to the file tmp60.mat

You show reshtemp1 and reshtmp60, is this a typo? ... what does reshtemp1 equal? — Flynn
Sure, and I understand that N is longitude, but what is N1, N2, N3? — Flynn
Also, shouldn't mean(reshtmp60(:,1)) be mean(reshtmp60(:,3))... the same column that you just removed the Nan's (999) from? — Flynn

EBH EBH · Accepted Answer · 2017-07-20T17:50:01

The main problem in your code is treating the nans. Observe the following example:

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

m =
       3.8421       3.6806

Both elements in m are simply the mean on all elements in a. But they are different! The reason is the taking the mean of all values together, with mean(a(:),'omitnan') is like summing all not-nan values, and divide by the number of values we summed:

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

but taking the mean of the first dimension, we get 6 mean values:

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

and when we take the mean of them we divide by a larger number, because all nans were omitted already:

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

Here is what I think you want in your code:

% 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)

To create the matrix you first permute the dimensions so the one you want to keep as is (time) will be the first. Then reshape the array to Tx(lon*lat), so you get 120 rows for all time steps and 259200 columns for all combinations of the coordinates. All that's left is to transpose it.

m1 is your first calculation, and m2 is what you try to do in the second one. They are equal here, but their value is not 5.855, even if I use your code.

However, I think the right solution will be to take the mean of all values together:

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

Reshape a 3D array and remove missing values

1 Answers