How to plot an average of two netcdf variables from different files?

Hello, thank you for taking time to read my question. I have a dataset for a uni coursework that I have downloaded from Nasa Copernicus website. All of the downloaded files are .nc format and I am trying to plot sea surface temp (sst) in the Baltic sea, I currently have monthly averages from 1993-2019. My problem arises when I try to combine or add or average two sst's from different months/years. The outcome I want is an map average of sst over the 20 or so odd years for specific month each year. How would I go about averaging 20 or so monthly averaged .nc files for each year? I want to use a month from every year since 1992 and have an average throughout the years? I have tried loading in multiple files however I am unsure how to go about this with a high volume of files and data like this? I created a working code to geoshow the monthly average which I will paste below, I have tried troubleshooting and tried various methods such as ncwriteschema and accumarray but I don't think I fully understand how the functions need to be used. If anyone could push me in th rright direction or explain a way I could do this it would be much appreciated. I have also trialled averaging per lat long degree bin however, I failed and couldnt not get it to work. I will percentage my trial shots at working it out in the code. Thank you so much for any help it is much appreciated.

clc; close all ;


ncfile = 'CMEMS_BAL_PHY_reanalysis_monthlymeans_201509.nc' ; %loads temperature file
ncdisp(ncfile) %info about ncfile

 depth = ncread(ncfile,'depth');
    lon = ncread(ncfile,'longitude'); %lat long
    lat = ncread(ncfile,'latitude');
    time = ncread(ncfile,'time'); nt = length(time);
    sst = ncread(ncfile, 'thetao');
    sst0 =sst(:, :, 1);
   
    sal = ncread(ncfile,'so'); %defines salinity variable

    figure('Color','white'); %creates figure


[Lat, Lon] = meshgrid(lat, lon);  % transform longitudes and latitudes into a grid of lon long x lat dimension
%depth = ncread(ncfile,'depth');
%lon = ncread(ncfile,'longitude');
%lat = ncread(ncfile,'latitude');
axes=worldmap([45 90],[0 45]); % loads world map with the limits for BALTICS
hold on

%lon_binID =  1 + floor( Lon ./ 1 ) ;
%lat_binID = 37 + floor( Lat ./ 1 ) ;
%binID     = [lon_binID, lat_binID] ;
%nBin      = 72 ;         

%sum_bin = accumarray( binID, sst0, [nBin, nBin] ) ;
%mean_bin = accumarray( binID, sst0, [nBin, nBin], @mean ) ;    this is where i tried dividing into bins and averaging but i wasn't able to get accumarray to work

%onesVector = ones( size(sst0,1), 1 ) ;
%count_bin  = accumarray( binID, onesVector, [nBin, nBin] ) ;
 
%sum_bin1  = accumarray( binID, data(:,2), [nBin, nBin] ) ;
%mean_bin1 = sum_bin1 ./ count_bin ;



%setm(axes,'mapprojection','mercator','Origin', [-180 0 180]) %changes projection and changes origin reference for coordinates
%setm(gca,'mapprojection','ortho')
levels = [ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 ]; % creates levels for contours
hold on

%nsst= reshape(sst, [523 383 56]);

cb = contourcbar(axes, 'Location', 'southoutside'); %creates colorbar
caxis([1 22]) %defines limits for colorbar
colormap(jet) %sets color scheme
hold on
  set(get(cb,'XLabel'),'String','SST (^oC)') %title for color bar
  set(cb,'Position',[0.25 0.05 0.55 0.03]) %adjust location and size of color bar
  hold on
  
%levels=[28]; %sets level for contour defining specific
%pcolorm(Lat,Lon,mean(sst0,3,'omitnan')) %average plot

geoshow(double(Lat),double(Lon),sst0,'DisplayType', 'contour','LevelList',levels,'LineColor','black') %plots temperature contour defining the western Pacific Warm Pool
hold on
geoshow(double(Lat),double(Lon), sst0,'DisplayType', 'texturemap') %creates surface map of temperature
hold on 
land = shaperead('landareas.shp', 'UseGeoCoords', true); %define land areas
geoshow(land, 'FaceColor', [0 0 0]) % plots land areas in black
hold on


%pcolorm(Lat,Lon,sst(:,:,3))
%time = ncread(ncfile,'time') ;  ; %reads the time variable and obtains its size
%% in case salinty wants to be used
%ncfileS = 'woa18_decav_s00_04.nc' ; %loads salinity file


%This is the file details 
Source:
           E:\2015\CMEMS_BAL_PHY_reanalysis_monthlymeans_201511.nc
Format:
           netcdf4_classic
Global Attributes:
           references            = 'http://www.smhi.se'
           institution           = 'Swedish Meterological and Hydrological Institute'
           history               = 'See source and creation_date attributees'
           Conventions           = 'CF-1.5'
           comment               = 'Provided by SMHI as a Copernicud Marine Environment Monitoring Service production unit'
           bullentin_type        = 'reanalysis'
           cmems_product_id      = 'BALTICSEA_REANALYSIS_PHY_003_011'
           title                 = 'CMEMS V4 Reanalysis: NEMO model 3D fields (monthly means)'
           deepest_depth         = 711.0592
           shallowest_depth      = 1.5014
           source                = 'SMHI reanalysis run NORDIC-NS2_1d_20151101_20151103'
           file_quality_index    = 1
           creation_date         = '2019-05-09 UTC'
           bullentin_date        = '20151101'
           easternmost_longitude = 30.2358
           northernmost_latitude = 65.8914
           westernmost_longitude = 9.0138
           southernmost_latitude = 48.4917
           stop_date             = '2015-11-01 UTC'
           start_time            = '00:00 UTC'
           start_date            = '2015-11-01 UTC'
           stop_time             = '00:00 UTC'
Dimensions:
           depth     = 56
           latitude  = 523
           longitude = 383
           time      = 1     (UNLIMITED)
Variables:
    depth    
           Size:       56x1
           Dimensions: depth
           Datatype:   single
           Attributes:
                       units         = 'm'
                       positive      = 'down'
                       standard_name = 'depth'
                       axis          = 'Z'
                       long_name     = 'depth'
    latitude 
           Size:       523x1
           Dimensions: latitude
           Datatype:   single
           Attributes:
                       units         = 'degrees_north'
                       long_name     = 'latitude'
                       standard_name = 'latitude'
                       axis          = 'Y'
    longitude
           Size:       383x1
           Dimensions: longitude
           Datatype:   single
           Attributes:
                       units         = 'degrees_east'
                       long_name     = 'longitude'
                       standard_name = 'longitude'
                       axis          = 'X'
    time     
           Size:       1x1
           Dimensions: time
           Datatype:   double
           Attributes:
                       units         = 'days since 1950-01-01 00:00:00'
                       long_name     = 'Validity time'
                       standard_name = 'time'
                       calendar      = 'gregorian'
                       axis          = 'T'
    sob      
           Size:       383x523x1
           Dimensions: longitude,latitude,time
           Datatype:   single
           Attributes:
                       _FillValue    = NaN
                       units         = '0.001'
                       long_name     = 'Sea water salinity at sea floor'
                       standard_name = 'sea_water_salinity'
                       missing_value = NaN
    bottomT  
           Size:       383x523x1
           Dimensions: longitude,latitude,time
           Datatype:   single
           Attributes:
                       _FillValue    = NaN
                       units         = 'degrees_C'
                       long_name     = 'Sea floor potential temperature'
                       standard_name = 'sea_water_potential_temperature_at_sea_floor'
                       missing_value = NaN
    mlotst   
           Size:       383x523x1
           Dimensions: longitude,latitude,time
           Datatype:   single
           Attributes:
                       _FillValue    = NaN
                       units         = 'm'
                       long_name     = 'Ocean mixed layer thickness defined by density (as in de Boyer Montegut, 2004)'
                       standard_name = 'ocean_mixed_layer_thickness_defined_by_sigma_theta'
                       missing_value = NaN
    so       
           Size:       383x523x56x1
           Dimensions: longitude,latitude,depth,time
           Datatype:   single
           Attributes:
                       _FillValue    = NaN
                       units         = '0.001'
                       long_name     = 'salinity'
                       standard_name = 'sea_water_salinity'
                       missing_value = NaN
    thetao   
           Size:       383x523x56x1
           Dimensions: longitude,latitude,depth,time
           Datatype:   single
           Attributes:
                       _FillValue    = NaN
                       units         = 'degrees_C'
                       long_name     = 'potential temperature'
                       standard_name = 'sea_water_potential_temperature'
                       missing_value = NaN
    uo       
           Size:       383x523x56x1
           Dimensions: longitude,latitude,depth,time
           Datatype:   single
           Attributes:
                       _FillValue    = NaN
                       units         = 'm s-1'
                       long_name     = 'eastward current'
                       standard_name = 'eastward_sea_water_velocity'
                       missing_value = NaN
    vo       
           Size:       383x523x56x1
           Dimensions: longitude,latitude,depth,time
           Datatype:   single
           Attributes:
                       _FillValue    = NaN
                       units         = 'm s-1'
                       long_name     = 'northward current'
                       standard_name = 'northward_sea_water_velocity'
                       missing_value = NaN

Answer 1

[ANKUR KUMAR]

If you have multiple netcdf files in a folder, and wish to take mean of SST, you can have a for loop to read the values, and subsequently take the mean.

F=dir('*.nc')
lon=ncread(F(1).name,'lon');
lat=ncread(F(1).name,'lat');
for i = 1:4
   sst=ncread(F(i).name,'thetao') ; 
   sst_values(:,:,i)=sst;
  % uncomment the below line if you wish to save the values in cell
   % sst_values_cell{i}=sst;    
end
contourf(lon,lat,nanmean(sst_values,3)')
% uncomment the below line if you have saved the values in cell
% contourf(lon,lat,nanmean(cat(3,sst_values_cell{:}),3)')