Loading netcdf in Matlab, era-interim and ncep-ncar

Purpose

This code is a good example of how to load a netcdf file into matlab. I give an example here for reanalysis data. Some reanalysis already come in the netcdf format, others may come in grib. I usually convert the grib files into netcdf before loading them with Matlab. Every reanalysis dataset uses different convention so the code needs to be adapted to the way your data is organized but also to the logic behind the file-names. With some modification, it can be used to load any netcdf of any reanalysis dataset.

Code for loading ERA-INTERIM in matlab

function [ value pressure latitude longitude] = get_ERA_INTERIM_all_levels(year,month,field)
 
folder=['/where_your_netcdf_on_pressure_levels_are/',num2str(year,'%04i'),'/'];
list=dir([folder,'ERAI.pl.',num2str(year,'%04i'),'-',num2str(month,'%02i'),'*']);
 
            for i=1:length(list)
                    filepath=[folder,'/',list(i).name];
                    ncid=netcdf.open(filepath,'NC_NOWRITE');
                    varidp=netcdf.inqVarID(ncid,'lv_ISBL1');
                    varidlat=netcdf.inqVarID(ncid,'g0_lat_2');
                    varidlon=netcdf.inqVarID(ncid,'g0_lon_3');
 
 
		if strcmp(field,'geopotential_height')
		textid='Z_GDS0_ISBL';
		end
		if strcmp(field,'temperature')
		textid='T_GDS0_ISBL';
		end
		if strcmp(field,'u')
		textid='U_GDS0_ISBL';
		end
		if strcmp(field,'v')
		textid='V_GDS0_ISBL';
		end
		if strcmp(field,'w')
		textid='W_GDS0_ISBL';
		end
		if strcmp(field,'vorticity')
		textid='VO_GDS0_ISBL';
		end
		if strcmp(field,'pv')
		textid='PV_GDS0_ISBL';
		end
 
		varid=netcdf.inqVarID(ncid,textid);
                temp=netcdf.getVar(ncid,varid);
                temp=permute(temp,[4 3 2 1]);
		temp=single(temp);
		if(strcmp(field,'geopotential_height'))
                	temp=temp/9.81;
            	end
                pressure=single(netcdf.getVar(ncid,varidp));
                latitude=single(netcdf.getVar(ncid,varidlat));
                longitude=single(netcdf.getVar(ncid,varidlon));
                value{i}=temp;			
                netcdf.close(ncid)
                avalueisset=1;                     
            end    
 
            value=cat(1,value{:});
            value=flipdim(value,2);
            pressure=flipdim(pressure,1);            
            pressure=single(pressure);
            latitude=single(latitude);
            longitude=single(longitude);
end

Code for loading NCEP-NCAR in matlab

function [ value pressure latitude longitude] = get_NCEP_NCAR(year,month,field)
 
%Field is the type of variable you want : geopotential_height,temperature,u,v,w
%Dimension of output value is (time,pressure_level,latitude,longitude)        
 
folder='/where_the_netcdf_are/';
 
%List the files in selected folder
air_list=dir([folder,'air.',num2str(year,'%04i'),'*']);
hgt_list=dir([folder,'hgt.',num2str(year,'%04i'),'*']);
omega_list=dir([folder,'omega.',num2str(year,'%04i'),'*']);
uwnd_list=dir([folder,'uwnd.',num2str(year,'%04i'),'*']);
vwnd_list=dir([folder,'vwnd.',num2str(year,'%04i'),'*']);
 
if(strcmp(field,'temperature'))
filename=air_list.name;
filepath=[folder,filename];
ncid=netcdf.open(filepath,'NC_NOWRITE');
varid=netcdf.inqVarID(ncid,'air');                
end  
 
if(strcmp(field,'geopotential_height'))
filename=hgt_list.name;
filepath=[folder,filename];
ncid=netcdf.open(filepath,'NC_NOWRITE');
varid=netcdf.inqVarID(ncid,'hgt');                                               
end
 
if(strcmp(field,'w'))
filename=omega_list.name;
filepath=[folder,filename];
ncid=netcdf.open(filepath,'NC_NOWRITE');
varid=netcdf.inqVarID(ncid,'omega');
end
 
if(strcmp(field,'u'))
filename=uwnd_list.name;
filepath=[folder,filename];
ncid=netcdf.open(filepath,'NC_NOWRITE');
varid=netcdf.inqVarID(ncid,'uwnd');
end  
 
if(strcmp(field,'v'))
filename=vwnd_list.name;
filepath=[folder,filename];
ncid=netcdf.open(filepath,'NC_NOWRITE');
varid=netcdf.inqVarID(ncid,'vwnd');
end
 
varidp=netcdf.inqVarID(ncid,'level');
varidlat=netcdf.inqVarID(ncid,'lat');
varidlon=netcdf.inqVarID(ncid,'lon');
varidt=netcdf.inqVarID(ncid,'time');
 
value=netcdf.getVar(ncid,varid,'short');
scale = netcdf.getAtt(ncid,varid,'scale_factor');
offset= netcdf.getAtt(ncid,varid,'add_offset');
value=double(value)*scale+offset;
value=single(permute(value,[4 3 2 1]));
 
pressure=single(netcdf.getVar(ncid,varidp));
latitude=single(netcdf.getVar(ncid,varidlat));
longitude=single(netcdf.getVar(ncid,varidlon));
%to keep only the time asked for
time=netcdf.getVar(ncid,varidt);     
dayssince111=time/24;
datevalue=dayssince111+datenum(1,1,1)-2;	
mm=datevec(datevalue); mm=mm(:,2);
value=value(month==mm,:,:,:);                
netcdf.close(ncid);
 
value(value>10^14)=nan;
 
%pads the data for the month with NAN if not all month is available yet
theoreticaltime=datenum(year,month,1):(datenum(year,month+1,1)-1);
ts=length(theoreticaltime)*4;
dummy=zeros([ts,length(pressure),length(latitude),length(longitude)]);
dummy=dummy*0/0;
dummy(1:size(value,1),:,:,:)=value;
value=dummy;
 
end

23 thoughts on “Loading netcdf in Matlab, era-interim and ncep-ncar

  1. Namaoui

    Hello
    I just wanted to know if this program is to interpolate the pressure values ​​for the model of ERA-interim.
    Otherwise thank you.

    Reply
    1. admin Post author

      This code does not interpolate the data. It simply loads in Matlab the data from Era-Interim already formatted on pressure levels.

      Does this answer your question?

      Reply
  2. smd

    Hi,

    whether the structure of this code will be the same for loading NCEP Reanalysis data on pressure levels for air temperature/ geopotential height or where I have to change?

    Reply
    1. admin Post author

      Thanks for asking, NCEP-NCAR data is structured in such a way that there is one file for each year and there are different files for different variables. I added on this page my code to load NCEP-NCAR. You can modify it to suit your own needs.
      I hope it helps,
      Patrick Martineau

      Reply
  3. smd

    How can I save one field variable in netcdf format in one degree latitude x one degree longitude(in NCEP-NCAR data spatial coverage is 2.5 degree) and for time ranging form 1948 to 2012 monthly in matlab? I am unable to create in netcdf format in matlab.

    Reply
    1. Patrick Martineau Post author

      The best way to read multiple netcdf files in matlab is probably to loop through every file and read them using a method similar to what I posted here. It is very important to close the NetCdf files after reading them using (netcdf.close(ncid)) as in my above examples. Opening too many netCdf in Matlab at the same time will crash the program. After reading the NetCdf in each step of the loop, I usually store variables in cells (value{i}=temperature). After the loop is done executing, cells can usually be easily concatenated into one single array.

      If you inspect further my sample “Code for loading ERA-INTERIM in matlab”, it is already looping through many files (one per day) and concatenates them for a full month.

      Reply
  4. Jason

    Hello Patrick

    I am thinking about using to Matlab to deal with netcdfs, currently I am using CDAT which is a python thing. I am now sure in Matlab how to reference or slice time periods, for instance slicing out the data in between 1979-1-1 and 1979-5-1. Because it could be tedious to count the time points to find out the corresponding indices for the time axis. And similarly if I need data within the Nino-3.4 region, could I use something like “nino=var(latitude=(-5,5),longitude=((360-170),(360-120))”?

    Reply
    1. Patrick Martineau Post author

      Hi Jason,
      In the sample code I gave here, the output named “value” is a 4d array with dimensions (time,pressure,latitude,longitude)
      The “time” “pressure” “latitude” and “longitude” are vectors with the values of the grid. (I loaded the time in the NCEP-NCAR example only, and you can see in the code how I picked one month only.
      To obtain only a subset of the large array you could do something like:

      new_value=value( time>=starttime & time < =endtime , pressure==500 , latitude>=-5 & latitude < =5 , longitude>=170 & longitude < =360)
      
      Once you learn about array manipulation in Matlab, it becomes a very powerful and customizable tool (Matlab reference on arrays)
      I also recommend using the Matlab functions to manipulate time
      Hope it helps,
      Patrick
      Reply
  5. NArun

    I had downloaded CMAP enhanced (with NCEP reanalysis) data for precipitation from the ESRL website which was done year wise. I read it using the basic function code as mentioned in help menu of Matlab. It gives two values for precipitation for each month. Now I am unable to understand why two values and what do each stand for? It would be grateful if you could clear my doubt as I am a beginner in this area. Thanks in advance.

    Reply
    1. Patrick Martineau Post author

      I am unable to help you since I’ve never used the dataset you mention. I advise to use ncdump -h filename.nc to see the description of the precipitation variables and sort out what they are.

      Reply
  6. smd

    Hi,
    Can you suggest me how to combine world ocean atlas(e.g., WOA2009) monthly and annual temperature/salinity data in one netcdf file in matlab? Actually monthly WOA data has 24 vertical level where as annual WOA data has 33 vertical level. I want to add the remaining vertical levels from annual WOA data to monthly WOA data so that monthly WOA data has 33 levels.

    Thanking you.

    Reply
  7. NArun

    I recently downloaded monthly data (NCEP/NCAR reanalysis 1) for temperature corresponding to the geographic extent of my interest . On opening the files I could see 32 point data each point with 73 values each. I am unable to interpret them, how is the data stored in netcdf files? It would great help if you could explain me what does it actually represent.

    Reply
    1. Patrick Martineau Post author

      To understand the structure of the NetCdf file, you should use ncdump -h filename, all necessary information to understand the structure of the data is well documented there, the dimensions and the variables available. From the netcdf file you should load the variable (like the temperature) and the arrays describing the grid, like latitude longitude pressure and time. By comparing the size of those coordinates arrays to the size of the variable array, you’ll figure out the dimensions. Ounce understood, Matlab is a very powerful tool to manipulate arrays.

      Reply
  8. NArun

    Can you help me provide a code for the latest version of Matlab for regridding of climate data. Like for example to reduce to the grid size of any GCmodel to that given by NCEP reanalysis.

    Reply
    1. Patrick Martineau Post author

      Hi Narun,
      To reduce grid size or interpolate to finer grid, I would look for the function interp2 in Matlab. The function is documented at that link http://www.mathworks.com/help/matlab/ref/interp2.html. You can use different interpolation methods like linear, nearest, spline and cubic. I think this would work well for reducing the grid size. It would require you to loop through pressure levels and interpolate horizontally at each level.

      Hope it helps,
      Patrick

      Reply
  9. Pascal Appelghem

    Hello Patrick
    Can you suggest me a function to convert degrees (included in netcdf files) in vector.
    I know converting U and V data, but I do not know the function to transform it into a vector value (0 to 360 ° in native netcdf files).
    Thank you in advance.
    Best regards.
    Pascal

    Reply
  10. Zohaib

    Dear Patrick
    I downloaded NCEP/NCAR Reanalysis monthly mean data separately for Geo-potential height, temperature and U-wind as .nc file. I want a time series of a month e.g. April, for whole time period i.e 1948 to present. Please suggest me any help

    Reply

Leave a Reply

Your email address will not be published. Required fields are marked *