How to extract seasonal means?

Hey MATLAB universe!

I have a file 360 x 40 x 444 (longitude x latitude x time). These are MONTHLY temperature values. The 444 are time values ("days since 1800"). I wanted to extract seasonal mean data for each year.

Example: Winter- Dec,Jan,Feb for each year and then take a mean for each year. So technically I should get a matrix of 360 x 40 x 37 (since this data set contains monthly values from January 1982 to December 2018). I have used the following code but I get wrong values if I map them:

lon=ncread('sstarc.nc','lon');
lat=ncread('sstarc.nc','lat');
sst=ncread('sstarc.nc','sst');
DJF = sort([[1:12:444],[2:12:444],[12:12:444]]);
Wacc = cumsum(sst(:,:,DJF));
Wavg = Wacc(:,:,3:3:end)/3;

I would like to know where I am going wrong. (Data file is attached and arranged as Jan-1982, Feb-1982,...Dec-2018.)

https://drive.google.com/open?id=12LAotPGE9WaGppXV8NYmr-fB5GqfXhW1

Thank you

Answer 1

[Cris LaPierre]

For those just coming here, you can see additional details and explanation here.

You seem to be grouping Jan, Feb and Dec from the same year together. I would be more inclined to group the winters as Dec yr1 with Jan yr2, Feb yr2. In doing so, I see the challenge being.

1. Grouping Dec from one year with Jan, Feb of the next
2. The first year will only have Jan, Feb while the last year will only have Dec.
3. If no data is removed, this means your will have data covereing 38 winters (Jan, Feb in last year as well as Dec in last year).

You asked what you were doing wrong. I don't think cumsum is doing what you want. The first value is itself. The second is the first plus the second. This pattern continues until the last value is the sum of all previous values plus itself. Instead, you want to sum each group separately.

Here's an approach that uses findgroups/splitapply

% Load data from nc file
lon=ncread('sstarc.nc','lon');
lat=ncread('sstarc.nc','lat');
sst=ncread('sstarc.nc','sst');

% Create vector of dates to use for extracting winter months
startDate = datetime(1982,1,1,"Format","MMM-uuuu");
endDate = datetime(2018,12,31,"Format","MMM-uuuu");
dates = startDate:calmonths(1):endDate;
winterMonth = month(dates) == 1 | month(dates) == 2 | month(dates) == 12;

% Extract winter data from sst
sstWinter = sst(:,:,winterMonth);

%% Create grouping variable for winter year 
% (treating consecutive Dec, Jan Feb as winter that year)
% Extract winter year as the year corresponding to January
winterYr = year(dates(winterMonth));
winterYr(month(dates(winterMonth))==12) = winterYr(month(dates(winterMonth))==12)+1;

% Find groups
G = findgroups(winterYr);
% Splitapply can't be used on 3D matrices. Reshaped to 2D to take the mean
sstTemp = reshape(sstWinter,size(sstWinter,1)*size(sstWinter,2),size(sstWinter,3));
f = @(x) mean(x,2,"omitnan");
tempYrMean = splitapply(f,sstTemp,G);
% reshape back to original number of rows, columns
sstYrMean = reshape(tempYrMean,size(sstWinter,1),size(sstWinter,2),[]);

%% Visualizing data
% Create variables for visualizing
[LAT, LON] = meshgrid(lat,lon);
winter1 = sstYrMean(:,:,1); % Change the index of the 3rd dimension to change the year

% Visualize the data for Winter 1982
geoscatter(LAT(:),LON(:),[],winter1(:))
colorbar(gca,"eastoutside")
title("Winter " + num2str(winterYr(1)))

[Keegan Carvalho]

Thank You Cris, I have understood most of this; the "sstYrMean" was the product that I wanted and happy that I can apply this to other seasons. I do have one last question: if I wanted to do an annual mean of the same dataset to get 37 values again (interannual data), should I use the following code:

winterMonth = month(dates) == 1 | month(dates) == 2 | month(dates) == 12;

But instead repeat the same from 1 to 12 like this:

annualMonth = month(dates) == 1 | month(dates) == 2 | month(dates) == 3 winterMonth = month(dates) == 4 | month(dates) == 5 | ... month(dates) == 12 |;
    

I think the above might give values in steps of 1,13,15 and so on. Or is there a simpler way?

[Cris LaPierre]

If you just wanted to keep it to within the same year, you would just have to remove this line. Everything else would stay the same.

winterYr(month(dates(winterMonth))==12) = winterYr(month(dates(winterMonth))==12)+1;

[Keegan Carvalho]

I'm sorry, but I didn't follow this part of the asnwer. I tried the code, got a map, but haven't understood it completely.

The following code:

winterMonth = month(dates) == 1 | month(dates) == 2 | month(dates) == 12;

Wouldn't I have to change the above code? Since it is coded for 1,2 and 12? (I thought the above code gives "winter months only"). So if I remove this:

winterYr(month(dates(winterMonth))==12) = winterYr(month(dates(winterMonth))==12)+1;

Would I get 37-year data from January-December each year?

[Cris LaPierre]

Winter months are Jan, Feb and Dec. Unless you want to change the season, don't change this line.

I grouped the data by year. The line I said to remove changes the year associated with December to the next year so that it groups with the contiguous Jan and Feb.

In pseudo code, it says "set the winter year of all months==12 (Decembers) to their current year + 1"

[Keegan Carvalho]

Okay I have folllowed this. Very grateful for your help, Cris. I can finally get back to my project with no obstructions. Thank you once again!!!

Answer 2

[Brian DeCicco]

I have a question on this topic. I'm doing something similar, but my dataset contains data at 6-hrly intervals daily from 1979-2018 (total of 58440 time steps). Just trying to extract Jan, Feb, Dec from the original time series. Here's how I have part of my code setup so far, but I am unsure how to point Matlab to know which month is which based on this setup. I appreciate any help you might be able to provide!

startDate = datetime(1979,1,1,0,0,0);

endDate = datetime(2018,12,31,18,0,0;

dates = startDate:hours(6):endDate;

winterMonth = month(dates) == 1 | month(dates) == 2 | month(dates) == 12;

sst_JFD = sst(:,:,winterMonth);

[Cris LaPierre]

The result of this code does not distinguish between the selected months. It's looking for any month.

The details about which month each value came from is stored in dates. You can see that with this code

dates(winterMonth)
% or
month(dates(winterMonth))

[Brian DeCicco]

Thanks Chris! So if I am trying to extract all times (6-hrly time steps per day) that are just assocaited with all Jan, Feb, and Dec each year, should my code look like this? I don't think that my "winterMonth" line is correct because I'm getting some mismatches for the month values (i.e. Feb is showing up as "1").

% Load data from nc file
lon=ncread('sst.nc','lon');
lat=ncread('sst.nc','lat');
sst=ncread('sst.nc','sst');

% Create vector of dates to use for extracting winter months
startDate = datetime(1979,1,1,0,0,0);
endDate = datetime(2018,12,31,18,0,0);
dates = startDate:hours(6):endDate;
winterMonth = month(dates) == 1 | month(dates) == 2 | month(dates) == 12;

sst_JFD = sst(:,:,winterMonth);