calculate mean of 4D array

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.

[Paul]

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]

>> ncdisp('198201.nc'); % display .nc file

Source:

E:\GURU_2024\metoffice_data_1982_2022_croped_SO\198201.nc

Format:

netcdf4_classic

Global Attributes:

CDI = 'Climate Data Interface version 2.0.4 (https://mpimet.mpg.de/cdi)'

Conventions = 'CF-1.0'

_NCProperties = 'version=2,netcdf=4.8.1,hdf5=1.10.7'

institution = 'UK Met Office'

title = 'Temperature and salinity analysis'

DSD_entry_id = 'UKMO-L4UHFnd-GLOB-v01'

references = 'Website and paper: https://www.metoffice.gov.uk/hadobs/en4/; Good, S. A., M. J. Martin and N. A. Rayner, 2013. EN4: quality controlled ocean temperature and salinity profiles and monthly objective analyses with uncertainty estimates, Journal of Geophysical Research: Oceans, 118, 6704-6716, doi:10.1002/2013JC009067'

contact = 'Rachel Killick - [e-mail]'

GDS_version_id = 'v1.7'

netcdf_version_id = '3.5'

creation_date = '2021-03-20 20:33:41.770 -00:00'

product_version = '1.0'

grid_resolution = ' 1.00000 degree'

start_date = '2001-01-01 UTC'

start_time = '00:00:00 UTC'

stop_date = '2001-01-01 UTC'

stop_time = '00:00:00 UTC'

southernmost_latitude = -90.5

northernmost_latitude = 89.5

westernmost_longitude = 0.5

easternmost_longitude = 362.5

file_quality_index = '0'

licence = 'EN4 is distributed under the Non Commercial Government Licence: http://www.nationalarchives.gov.uk/doc/non-commercial-government-licence/version/2/. The data are available for non-commercial use with attribution to the data providers, please see the product website (see global attribute: references) for terms and conditions.'

history = 'Fri May 10 12:04:04 2024: cdo sellonlatbox,-180,180,-75,0 EN.4.2.2.f.analysis.g10.198201.nc EN.4.2.2.f.analysis.g10.198201_SO.nc'

CDO = 'Climate Data Operators version 2.0.4 (https://mpimet.mpg.de/cdo)'

Dimensions:

time = 1 (UNLIMITED)

bnds = 2

lon = 360

lat = 76

depth = 42

Variables:

time

Size: 1x1

Dimensions: time

Datatype: single

Attributes:

standard_name = 'time'

bounds = 'time_bnds'

units = 'days since 1800-01-01 00:00:00'

calendar = 'standard'

axis = 'T'

time_bnds

Size: 2x1

Dimensions: bnds,time

Datatype: double

lon

Size: 360x1

Dimensions: lon

Datatype: single

Attributes:

standard_name = 'longitude'

long_name = 'longitude'

units = 'degrees_east'

axis = 'X'

lat

Size: 76x1

Dimensions: lat

Datatype: single

Attributes:

standard_name = 'latitude'

long_name = 'latitude'

units = 'degrees_north'

axis = 'Y'

depth

Size: 42x1

Dimensions: depth

Datatype: single

Attributes:

standard_name = 'depth'

long_name = 'depth'

units = 'metres'

positive = 'down'

axis = 'Z'

bounds = 'depth_bnds'

depth_bnds

Size: 2x42

Dimensions: bnds,depth

Datatype: single

temperature

Size: 360x76x42x1

Dimensions: lon,lat,depth,time

Datatype: single

Attributes:

standard_name = 'sea_water_potential_temperature'

long_name = 'temperature'

units = 'kelvin'

add_offset = 273.15

scale_factor = 1

_FillValue = -32768

missing_value = -32768

salinity

Size: 360x76x42x1

Dimensions: lon,lat,depth,time

Datatype: single

Attributes:

standard_name = 'sea_water_salinity'

long_name = 'salinity'

units = '1'

_FillValue = -32768

missing_value = -32768

temperature_uncertainty

Size: 360x76x42x1

Dimensions: lon,lat,depth,time

Datatype: single

Attributes:

long_name = 'temperature error standard deviation'

units = 'kelvin'

_FillValue = -32768

missing_value = -32768

salinity_uncertainty

Size: 360x76x42x1

Dimensions: lon,lat,depth,time

Datatype: single

Attributes:

long_name = 'salinity error standard deviation'

units = '1'

_FillValue = -32768

missing_value = -32768

temperature_observation_weights

Size: 360x76x42x1

Dimensions: lon,lat,depth,time

Datatype: single

Attributes:

long_name = 'temperature observation weights'

_FillValue = -32768

missing_value = -32768

comment = 'The total weighting given to the observation increments when forming the analysis'

salinity_observation_weights

Size: 360x76x42x1

Dimensions: lon,lat,depth,time

Datatype: single

Attributes:

long_name = 'salinity observation weights'

_FillValue = -32768

missing_value = -32768

comment = 'The total weighting given to the observation increments when forming the analysis'

Answer 1

[Nivedita]

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.

1. We will first initialize the annual_mean and seasonal_mean variables to accumulate temperature data.
2. Then we will loop through each year from 1982 to 2022 and each month from January to December.
3. For each month, we will construct the filename and load the temperature data using ncread.
4. We then accumulate the temperature data for annual and seasonal means.
5. 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).
6. 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/512283-how-to-extract-seasonal-means

https://www.mathworks.com/matlabcentral/answers/829998-mean-max-min-map-from-multiple-netcdf-file

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!

Answer 2

[Kris Fedorenko]

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