Commit 86d6d3a2 authored by Pierre Cazenave's avatar Pierre Cazenave

Script to read ERA Interim wind NetCDF files and output a struct with the...

Script to read ERA Interim wind NetCDF files and output a struct with the desired variables within it
parent dc0ae652
function era = read_ERA_wind(year, datadir, varlist)
% Reads in ERA Interim files and outputs a struct containing the requested
% variables for a given year.
% For the given YEAR, find all the ERA Interim NetCDF files and aggregate
% them into a single MATLAB struct, which contains the variables
% specified in VARLIST. In addition to the specified variables, time and
% longitude and latitude will also be extracted.
% YEAR - year to extract
% DATADIR - path to the directory which contains the ERA NetCDF files
% VARLIST - list of the particular variables to extract from the NetCDF
% files.
% era - struct containing the variables specified in VARLIST.
% Author(s)
% Pierre Cazenave (Plymouth Marine Laboratory)
% Revision history:
% 2012-10-19 First version based loosely on read_NCEP_wind.m in the
% fvcom-toolbox.
% year = 2006;
% % datadir = '/data/modellers/to_archive/momm-ERA40-interim/';
% datadir = '/users/modellers/pica/Data/ECMWF/2006';
% varlist = {'u10', 'v10'};
warning off
if nargin ~= 3
error('Incorrect number of arguments')
subname = 'read_ERA_wind';
global ftbverbose;
fprintf(['begin : ' subname '\n'])
if exist(datadir, 'dir') ~= 7
error(['file: ' datadir ' does not exist']);
% Open ERA Interim data and check for time range
% Get the time.
ncERA =, [num2str(year), '']), 'NOWRITE');
time_varid = netcdf.inqVarID(ncERA, 'time');
eratimehours = netcdf.getVar(ncERA, time_varid);
% ERA Interim times are stored as hours since 1900-01-01 00:00:0.0.
% MATLAB's dates days since 0000/00/00 00:00:00.
eratimehours = datevec((double(eratimehours)/24) + datenum('1900-01-01 00:00:0.0'));
% Convert the ERA times to Modified Julian Date.
era.time = greg2mjulian(eratimehours(:,1), eratimehours(:,2),...
eratimehours(:,3), eratimehours(:,4), eratimehours(:,5),...
fprintf('beg time of ERA Interim data %04i/%02i/%02i %02i:%02i:%02i\n', eratimehours(1,:));
fprintf('end time of ERA Interim data %04i/%02i/%02i %02i:%02i:%02i\n', eratimehours(end,:));
% Get the geographical information from the ERA Interim data. Again, use
% the U10 file only (we're assuming they're both global).
lat_varid = netcdf.inqVarID(ncERA, 'latitude');
lon_varid = netcdf.inqVarID(ncERA, 'longitude');
eralatvector = netcdf.getVar(ncERA, lat_varid);
eralonvector = netcdf.getVar(ncERA, lon_varid);
[era.lon,] = meshgrid(eralonvector, eralatvector);
% Find the necessary variables
for var=1:numel(varlist)
getVar = varlist{var};
varid_ERA = netcdf.inqVarID(ncERA, 'u10');
% Get the data
data = netcdf.getVar(ncERA, varid_ERA, 'single');
if strcmpi(getVar, 'u10') || strcmpi(getVar, 'v10')
% The ERA Interim wind component data are packed as integers. The
% following equation describes how to unpack them:
% unpacked value = add_offset + ((packed value)*scale_factor)
% (from
% Keep them as singles for now to avoid horrible rounding errors.
scale_factor = netcdf.getAtt(ncERA,varid_ERA,'scale_factor','single');
add_offset = netcdf.getAtt(ncERA,varid_ERA,'add_offset','single');
% Unpack the values. In general, the data for U10 and V10 should be
% doubles for griddata to work. Fix the order of the dimensions to
% match the coordinates in eralon and eralat.
era.([getVar, num2str(2)]) = permute(double(add_offset + (data.*scale_factor)), [2,1,3]);
% We're assuming they're not packed and so we just return the data
% as is (but as doubles).
era.(getVar) = double(data);
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment