JISAO data

MATLAB scripts for reading and writing COARDS netCDF files

Simple | Sophisticated

Converting 3-dimensional arrays to 2-dimensions




Simple

Both the United States Geological Survey (USGS) and Australia's CSIRO WWW sites provide free MATLAB scripts for reading netCDF files. In the following I describe the USGS software. Please send me some email if you have tried the CSIRO scripts. My impression is that the USGS code is now a part of the standard MATLAB software. To see if this software is installed on your machine, type help netcdf in a MATLAB session and see if it provides you with information. This software is installed on the University of Washington Atmospheric Sciences / JISAO computers. The following are quick instructions to read files and two examples. With luck, you can actually be reading data in 15 minutes. Rene Garreaud provided this script. Contact Todd Mitchell if you have problems.

1) Let filename.nc be the netCDF file that you want to read. Type ncdump -h filename.nc in UNIX (not MATLAB) to see the header of the netCDF file. This will tell you the variable names and indicate whether or not the variables need to be unpacked when you read them into MATLAB. The outputs of this call will be the dimensions and variables.

  • Dimensions tells you the variable names and the dimensions of each variable. For example, lat (short for "latitude") could be a vector of latitude grid points, and there could be 180 of them for 2-degree latitude resolution global data.
  • Variables provides information on the variables. A variable will be introduced as
    variable_type variable_name, for example, float lat. The variable name is then lat. The scale_factor and add_offset for the variable, if they are defined, are used to unpack the data:
    unpacked_data = ( packed_data * scale_factor ) + add_offset;

    2) Type nc = netcdf('filename.nc','nowrite'); in MATLAB to open the netCDF file.

    3) The basic command you need in MATLAB is variable=nc{'variable_name'}(:,:,:,:);
    where variable_name is the variable that you want to read in. The (:,:,:,:) at the end specifies the dimensions (for example, time,vertical level,lat,lon) that you want to read in. (:,:,:,:) will read in all of a 4-dimensional data volume. The order of the arguments is given in the ncdump output. For NCEP-NCAR reanalysis data provided by CDC, the vertical level ranges from 1,2,3, ... (don't type in 500 hoping to get 500mb data). Often the data has been scaled to make for a physically smaller file. The best example is sea-level pressure (SLP), where a typical value of 1025.35mb will be written as 2525 in the file. The scaling information will be included as "add_offset" and "scale_factor" in the attributes for a variable. Matlab will rescale the data if you add a "1" to the above command: variable=nc{'variable_name',1}(:,:,:,:);

    Two examples:

    Example 1:

    ncdump -h sstoi8299.nc
    netcdf sstoi8299 {
    dimensions:
            lat = 180 ;
            lon = 360 ;
            time = UNLIMITED ; // (216 currently)
    variables:
            float lat(lat) ;
                    lat:title = "Latitude" ;
                    lat:units = "degrees_north" ;
                    lat:scale_factor = 1.f ;
                    lat:add_offset = 0.f ;
            float lon(lon) ;
                    lon:title = "Longitude" ;
                    lon:units = "degrees_east" ;
                    lon:scale_factor = 1.f ;
                    lon:add_offset = 0.f ;
            double time(time) ;
                    time:title = "Time" ;
                    time:units = "days    since 1982- 1- 1  0: 0: 0" ;
                    time:scale_factor = 1.f ;
                    time:add_offset = 0.f ;
            short data(time, lat, lon) ;
                    data:long_name = "deg.C" ;
                    data:add_offset = 20.f ;
                    data:scale_factor = 0.01f ;
                    data:missing_value = 32767s ;
                    data:units = "deg.C" ;
    
    nc = netcdf('sstoi8299.nc','nowrite');
    ygrid=nc{'lat'}(:);  
    xgrid=nc{'lon'}(:);  
    sst=nc{'data',1}(1,:,:);    
    
    Reads the latitude and longitude grid points, and the first month of SST data into MATLAB. The ,1 argument specifies that the SST data is unpacked using the add_offset and scale_factor information (see the ncdump output).

    Example 2:

    ncdump -h air.mon.mean.nc 
    netcdf air.mon.mean {
    dimensions:
            lon = 144 ;
            lat = 73 ;
            level = 17 ;
            time = UNLIMITED ; // (484 currently)
    variables:
     ...
            short air(time, level, lat, lon) ;
    
    
    That is, we have a global grid of air temperature, with 17 levels and 484 time samples (monthly from 1958-1997). The variable air have four dimensions: time, vertical level, lat and lon.
    We want to read the global fields at the 3rd vertical level at time=21 Then, try:
     nc = netcdf('air.mon.mean.nc','nowrite');
     T=nc{'air'}(21,3,:,:);
     whos
      Name      Size         Bytes  Class
    
      T        73x144        84096  double array
      nc        4-D           1690  netcdf object
    
    The result is in the matlab variable T
    
    
    
    Sophisticated

    The netCDF files written written to the COARDS standard include a time variable (year, month, day, hour, ...), which is nice to have and time consuming to set up. The rest of this page describes libraries that enable time and the other variables to be read more easily.

    Your MATLAB code will need both the MexCDF and MexEPS libraries to run these scripts. I apologize for this inconvenience, but I am a scientist and not a programmer, so I am forced to use everyone elses libraries. Use the following code at your own risk.

    The scripts comply with the netCDF standard of the Cooperative Ocean-Atmosphere Research Data Service (COARDS). The metadata information (header) for a simple file written to COARDS specifications is here.

    The COARDS standard is for time to be written as "so many units of time since some reference time." Writing and reading time in this way complies with the udunits standard, and keeping track of leap years and historical changes in calendars requires the MexEPS libraries, which are expertly managed by Willa Zhu of NOAA PMEL. The "nc_read.m" script described below will also read files where time is written as "yyyymmddhhmmss" (year, month, day, ...), which is an older in-house convention. The "nc_write.m" script writes the time variable in the COARDS fashion.

    Reading a map:
    [ data, i4, r4, nam, xgrid, ygrid ] = nc_read( filename, record);

    where:
    "data" is a vector of the data,
    "i4, r4" are vectors containing metadata: number of latitude and longitude points, first latitude point, etc.,
    "nam" (optional) variable name
    "xgrid, ygrid" (optional) vectors of longitude and latitude points, respectively
    "filename" name of file to be read in
    "record" number of the record to be read in

    The header of the MATLAB script explains these variables in more detail.

    nc_read.m -- the MATLAB script.

    Writing a map:
    nc_write( filename, record, data, nam, i4, r4 );

    where
    "filename" is the output filename,
    "record" is the record number to write,
    "data" is a vector containing the data,
    "nam" a variable name or units information,
    "i4, r4" are vectors containing metadata: number of latitude and longitude points, first latitude point, etc.,

    The header of the MATLAB script explains these variables in more detail.

    nc_write.m -- the MATLAB script.


    Matlab now handles arrays with greater than 2 dimensions, which is useful for calculating means along latitude or longitude circles or spatial differences. The great majority of calculations, however, treat each grid point as an independent timeseries, and it is simpler to manipulate the data as a two-dimensional matrix, A(t,x). The following 2 lines will reduce an input 3-dimensional array to 2-dimensions.

    Consider A = A(nt,ny,nx) where nt, ny, nx are the number of time, latitude, and longitude points, respectively.
    A = permute( A, [ 3 2 1 ] );
    A = squeeze( reshape( A, nx*ny,nt,1 ) )' ;
    "A" will be dimensioned A(nt,nx*ny), and the second index will first span the first latitude circle and then the second latitude circle, and so forth.

    threetotwo.m -- the MATLAB script.


    January 2010: I am at present figuring out how to get matlab and netcdf-input/output to work on my Apple computer. The Regional Ocean Modeling System (ROMS) WWW page seems very good so far. I have been able to get "snctools", a set of scripts written on top of low order input/output routines, to work. Yay!

    January 2009
    Todd Mitchell ( mitchell@atmos.washington.edu )
    Rene Garreaud
    JISAO data