calculate mean of 4D array
Show older comments
I have 41 years monthly data of temperature (lat,lon,depth,time,temperature). I wanted to calculate annual mean and seasonal mean [(summer: Dec, Jan, Feb), (autumn: Mar, Apr, May), (winter = Jun, Jul, Aug), (spring = Sep, Oct, Nov)] of each grid, eache depth level for the period of 1982-2022.
2 Comments
Paul
on 11 May 2024
Hi @Gurumoorthi K.,
You're more likely to get help if you upllad your data (or something like it) in .mat file using the Paperclip icon on the Insert menu so that people can see exactly how the data is organized.
Gurumoorthi K
on 11 May 2024
Moved: Voss
on 11 May 2024
Accepted Answer
Nivedita
on 23 Jul 2024
Hi Gurumoorthi,
As I do not have access to the data,I am providing the answering based on my undertanding from your query. You can make necessary modifications according to the data available with you. We can achieve your requirement by following the processes which involves loading the data, reshaping it, and then calculating the means for each grid point and depth level.
- We will first initialize the annual_mean and seasonal_mean variables to accumulate temperature data.
- Then we will loop through each year from 1982 to 2022 and each month from January to December.
- For each month, we will construct the filename and load the temperature data using ncread.
- We then accumulate the temperature data for annual and seasonal means.
- After accumulating the data, we will calculate the annual mean by dividing by the total number of months (41 years * 12 months). Similarly, we will also calculate the seasonal means by dividing by the number of months in each season (41 years * 3 months per season).
- Finally, we will save the annual and seasonal means to .mat files.
% Define the path to your data files
dataDir = 'E:\GURU_2024\metoffice_data_1982_2022_croped_SO\';
% Initialize variables to accumulate the data
annual_mean = [];
seasonal_mean = struct('summer', [], 'autumn', [], 'winter', [], 'spring', []);
% Loop through each year and month
for year = 1982:2022
for month = 1:12
% Construct the filename
filename = sprintf('%s%d%02d.nc', dataDir, year, month);
% Load the temperature data
temperature = ncread(filename, 'temperature');
% Accumulate the temperature data for annual mean
if isempty(annual_mean)
annual_mean = temperature;
else
annual_mean = annual_mean + temperature;
end
% Determine the season and accumulate the temperature data for seasonal mean
if ismember(month, [12, 1, 2])
% Summer
if isempty(seasonal_mean.summer)
seasonal_mean.summer = temperature;
else
seasonal_mean.summer = seasonal_mean.summer + temperature;
end
elseif ismember(month, [3, 4, 5])
% Autumn
if isempty(seasonal_mean.autumn)
seasonal_mean.autumn = temperature;
else
seasonal_mean.autumn = seasonal_mean.autumn + temperature;
end
elseif ismember(month, [6, 7, 8])
% Winter
if isempty(seasonal_mean.winter)
seasonal_mean.winter = temperature;
else
seasonal_mean.winter = seasonal_mean.winter + temperature;
end
elseif ismember(month, [9, 10, 11])
% Spring
if isempty(seasonal_mean.spring)
seasonal_mean.spring = temperature;
else
seasonal_mean.spring = seasonal_mean.spring + temperature;
end
end
end
end
% Calculate the annual mean by dividing the accumulated sum by the number of months
annual_mean = annual_mean / (41 * 12);
% Calculate the seasonal means by dividing the accumulated sums by the number of months in each season
seasonal_mean.summer = seasonal_mean.summer / (41 * 3);
seasonal_mean.autumn = seasonal_mean.autumn / (41 * 3);
seasonal_mean.winter = seasonal_mean.winter / (41 * 3);
seasonal_mean.spring = seasonal_mean.spring / (41 * 3);
% Save the results to .mat files
save('annual_mean.mat', 'annual_mean');
save('seasonal_mean.mat', 'seasonal_mean');
Make sure to adjust the dataDir variable to point to the directory containing your NetCDF files. This script assumes that the files are named in the format YYYYMM.nc. If your file naming convention is different, you may need to adjust the filename construction accordingly.
For more information on handling NetCDF files and computing means over different dimensions, please refer to these following community queries:
https://www.mathworks.com/matlabcentral/answers/203585-how-to-calculate-average-of-netcdf-daily-data
For more information on the "ncread" function, type the following command in the command window:
web(fullfile(docroot, 'matlab/ref/ncread.html'))
Hope this helps!
More Answers (1)
Kris Fedorenko
on 27 Jan 2026
Edited: Kris Fedorenko
on 27 Jan 2026
I know this is an old post, but I was curious and taking a look at it anyway, so I thought I'd post my stab at it just in case.
Would have been great to have the data attached. I wasn't able to find exactly the same data, but I found similar one at https://www.metoffice.gov.uk/hadobs/en4/download-en4-2-2.html . I only downloaded data for years 2021 and 2022 (for all the months).
In my approach I first concatenate all the temperature data into longitude-by-latitude-by-depth-by-date matrix, so this is more memory-intensive than the other proposed solution. But this would be useful if there is need to calculate any other statistics over the same data.
ncFiles = dir('*.nc');
% this will be dimension vector for the 4th dimension (dates)
dates = NaT(length(ncFiles),1);
temperatures = [];
for i = 1:length(ncFiles)
filename = ncFiles(i).name;
% extract date from filename and add it to the dimension vector
[~, dateText, ~] = fileparts(filename);
dates(i) = datetime(dateText, InputFormat="yyyyMM");
% alternatively extract date from the file
% daysSince1800_01_01 = ncread(filename, "time");
% dates(i) = datetime(1800, 01, 01) + days(daysSince1800_01_01);
% extract longitudes from the file (dimension vector for 1st dimension)
curLons = ncread(filename, "lon");
% check if longitudes are the same in each file;
% might be unnecessary if we can already assume that.
if i==1
longitudes = curLons;
else
assert(isequal(longitudes, curLons), "Longitudes are not the same in all files")
end
% extract latitudes from the file (dimension vector for 2nd dimension)
curLats = ncread(filename, "lat");
% check if latitudes are the same in each file;
% might be unnecessary if we can already assume that.
if i==1
latitudes = curLats;
else
assert(isequal(latitudes, curLats), "Latitudes are not the same in all files")
end
% extract depths from the file (dimension vector for 3rd dimension)
curDeps = ncread(filename, "depth");
% check if depths are the same in each file;
% might be unnecessary if we can already assume that.
if i==1
depths = curDeps;
else
assert(isequal(depths, curDeps), "Depths are not the same in all files")
end
% extract temperature and concatenate it all together
% (along 4th demension)
curTemp = ncread(filename, 'temperature'); % lon x lat x depth x time
temperatures = cat(4, temperatures, curTemp);
end
Now that we have all temperature data in one matrix and also the dimension vectors for Longitude, Latitude, Depth, and Date, we can calculate any needed statistic from it (like the mean for any given year). We can also extract data for any specific Lat-Long-Depth grid point. For example:
% yearly mean for 2021 (for each lat-lon-depth grid)
mean2021 = mean(temperatures(:,:,:, dates.Year==2021), 4); % mean over dates (4th dim), the result is long x lat x depth
% 2021 yearly mean for Latitude = 71, Longitude = 42, at Depth closest to 15 meters
[~, closestInd] = min(abs(depths-15));
mean2021(longitudes==42, latitudes==71, closestInd) % ~277.4722 Kelvin
% summer mean over all years
summer = [12, 1, 2]; % Dec, Jan, Feb - summer months in southern hemisphere
meanSummers = mean(temperatures(:,:,:, ismember(dates.Month, summer)), 4);
% summer mean for Latitude = 71, Longitude = 42, at Depth(s) close to 110 meters
meanSummers(longitudes==42, latitudes==71, isapprox(depths, 110, "veryloose")) % ~275.1186 Kelvin
% plot summer mean for Latitude = 71, Longitude = 42, at all depths
meanSummersAtCoords = squeeze(meanSummers(longitudes==42, latitudes==71, :)); % depth x 1
plot(depths, meanSummersAtCoords, 'x-')
title('Mean temperature over Jan, Feb, Dec months 2021-2022 at 71 N, 42 E')
xlabel("Depth (meters)")
ylabel("Temperature (Kelvin)")
Alternatively, you could also flatten the data into a timetable and calculate any needed statistics using groupsummary or rowfun (for an example of a timetable approach see https://uk.mathworks.com/matlabcentral/answers/2176414-era5-average-hourly-data-and-yearly-data?s_tid=srchtitle ).
Categories
Find more on Visualization and Analytics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
