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

24 views (last 30 days)
Daryna Butash on 10 Mar 2021
Edited: ANKUR KUMAR on 10 Mar 2021
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
##### 0 CommentsShowHide -1 older comments

Sign in to comment.

### Accepted Answer

ANKUR KUMAR on 10 Mar 2021
Edited: ANKUR KUMAR on 10 Mar 2021
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)')
##### 0 CommentsShowHide -1 older comments

Sign in to comment.

R2020a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by