How to get decade means of monthly 3D satellite data?

Hello Matlab World!

I have been trying to convert the monthly netcdf data (file attached "bob") to decade data and later plot a map of the same. I have visited a similar question asked before and tried the following code, but seem to get the wrong results:

lon=ncread('bob.nc','lon');                             % longitude
lat=ncread('bob.nc','lat');                             % latitude
time=ncread('bob.nc','time');                           % time
temp=ncread('bob.nc','sst');
t=size(time,3);
[groups, groupID] = findgroups(floor(t/12)*10);

The temperature variable, temp= 17 x 15 x 444 (where 444 is the monthly data from January 1986 to December, 2019). So I basically need to get 4 groups (1986-95, 96-2005, 2006-2015 and 2016-2019). The above code gives 1 group which is obviously not correct.

Would be grateful to receive some help in this.

Answer 1

[Cris LaPierre]

Do you want a single temperature for each decade, or you you want an average for each lat/lon coordinate pair?

I wasn't able to get splitapply or groupsummary to work with a 3D dataset, but I was able to achieve group means with a for-loop.

The challenge here is perhaps identifying the groups. I want to convert your time variable to a datetime so I can group by year. The following code allows we to inspect the attributes of time.

info = ncinfo('keeganCarvalho_bob.nc')

From here, I can see that time is the number of days since 1800-1-1 00:00:00, and that is 1 day. I performed the convertion with the following code.

dTime = datetime(1800,1,1)+days(time)

I then extracted the year and used the discretize function to group the years into your specified bins.

yr = year(dTime);
G = discretize(yr,[1986 1996 2006 2016 2019]);

All that is left is to get the mean.

for d = 1:max(G)
    decadeMean(d) = mean(temp(:,:,G==d),"all");
end
plot(decadeMean)

If instead you want the mean for each lat/lon location, you could try something like this. You'd have to change the plotting code to be something appropriate for matrix data.

for d = 1:max(G)
    decadeMean(:,:,d) = mean(temp(:,:,G==d),3);
end
% rearrange data so decade means are in rows, and columns are lat/lon pairs
plot(reshape(permute(decadeMeanLL,[3,1,2]),max(G),[]))

[Keegan Carvalho]

Thanks Cris,

Yes, I wanted to get a single temperature for each decade at each coordinate (lon-lat). This pretty much gives me what I need. But the end product (map) is a problem. I got an error when I reach this:

for d = 1:max(G)
    decadeMean(:,:,d) = mean(temp(:,:,G==d),3);
end

And thats becasue the "decadeMean" is a 1 by 4 matrix while the right side matrix is 17 by 15. Also I tried the last code but get a 2d plot.

I was wondering if there was a way to loop the group means in the "temp" variable directly?

Another part I didn't quite understand. In the above codes, shouldn't it be

d = 1:mean(G)....
    

Since, I want to find the mean for each decade?

[Cris LaPierre]

That error means the variable already exists. You likely ran the code I shared previously. If you want to keep both, change this variable name. Otherwise, initialize the variable before running the for-loop.

decadeMean = zeros(size(temp,1),size(temp,2),max(G));

As for your second question, no, it should be max. This allows it to loop through all the groups and calculate the mean for each one. You can learn more about for loops by completing Chapter 13 of MATLAB Onramp.

[Cris LaPierre]

I don't understand your question for "a way to loop the group means in the "temp" variable directly?"

What do you envision? What would the final size of temp be? It is currently 17x15x408. You mention wanting a single temperature for each decade at each coordinate. That is what decadeMean contains, with a size of 17x15x4.

[Keegan Carvalho]

Apologies for my bad phrasing...

What I meant was to get a matrix such as 17 x 15 x 4 (which I managed to get). I did try the above code (revised) but the 3rd dimension of decadeMean variable is 0 for all "4" decades (17 x 15 x "4"). Here's my attempts so far:

dTime = datetime(1800,1,1)+days(time);
yr = year(dTime);
G = discretize(yr,[1986 1996 2006 2016 2019]);
for d = 1:max(G)
    decadeMean(d) = mean(temp(:,:,G==d),"all");
end
    for d = 1:max(G)
    decadeMean = zeros(size(temp,1),size(temp,2),max(G));
end

[Cris LaPierre]

Looks like I'm the one confusing you now. I had proposed two different solutions. It looks like you are trying to combine pieces of both, which will not work. Now that I know what it is you want, let me propose a single solution.

dTime = datetime(1800,1,1)+days(time);
yr = year(dTime);
G = discretize(yr,[1986 1996 2006 2016 2019]);

decadeMean = zeros(size(temp,1),size(temp,2),max(G));
for d = 1:max(G)
    decadeMean(:,:,d) = mean(temp(:,:,G==d),3);
end

[Keegan Carvalho]

Thanks Cris!

It was my bad but I understood the code now. Thank you once again for your continued help. Cheers!!