%*************************************************************************************************
% Matlab Script to read and plot NCEP/NCAR Reanalysis data in NetCDF for 4x daily data, either 3D or 4D. Made by Paula Doubrawa Moreira @ University of Alaska Fairbanks.
%
% Need to have MexCDF SNCTools (very easy to download and install! Don't forget to set them on your matlab path with the command 'addpath' on your matlab window).
%
% You have to modify somethings according to the times you want to plot. Look through for directions and examples, you'll find them. Don't forget that in order to run this script, you need to type the full path to it on the matlab window or run matlab from the same folder where you saved it, so matlab can find it! (or try saving it to your matlab path!)
%*************************************************************************************************
%
% Clear workspace. This command "reinitializes" matlab.
%
clear
%
% Ask for file to be opened. Path typed will be stored in variable "file". Write the whole path! Ex: /home/john/files/air.2009.nc. This next command will store the input on a variable named "file", so you don't have to input the long path again.
%
file=input('Enter path to the file you wish to open: ','s');
%
% See variables/dimensions/attributes
%
nc_dump(file)
%
% Open file using Mathwork's NetCDF package. This will generate a NetCDF ID (ncid).
%
ncid=netcdf.open(file,'NC_NOWRITE');
%
% Inquire about the number of dimensions (Will be 3 for surface data files and 4 for files that have data on various pressure levels). Dimension 0 is time, 1 and 2 are lon and lat. Dont get confused... 3 dimensions with ids 0,1 and 2. Variable 3 is not a dimension, but the variable! (Variable 4 if you have more than 1 z levels)
%
[numdims,numvars]=netcdf.inq(ncid);
%
% Get name of the variable that will be plotted. Variable ID will be 3 if it's a 3D file. Itll be 4 if its a 4D file. Name of variable will be returned under variable name var (in square brackets below).
%
if numvars==4
[var]=netcdf.inqvar(ncid,3);
elseif numvars==5
[var]=netcdf.inqvar(ncid,4);
end
%
% Define dimensions. Only define level dimension for 4D files.
%
if numdims==4
lvl=nc_varget(file,'level');
% Define lon and lat for all files.
lon=nc_varget(file,'lon');
lat=nc_varget(file,'lat');
disp('**************You chose a 4D file**************')
elseif numdims==3
% Define lon and lat for all files.
lon=nc_varget(file,'lon');
lat=nc_varget(file,'lat');
disp('**************You chose a 3D file**************')
end
%
% In case of a 4D file, ask for pressure level.
%
if numdims==4
nivel=input('Enter pressure level you wish to plot (in hPa): ');
%
%---------------------------------------------------------------------------------------------------
% LOOP IN CASE YOU HAVE A 4D FILE
%---------------------------------------------------------------------------------------------------
%
switch nivel
case 1000
z=1
case 925
z=2
case 850
z=3
case 700
z=4
case 600
z=5
case 500
z=6
case 400
z=7
case 300
z=8
case 250
z=9
case 200
z=10
end
%**************************************************************************************************
% This is the TIME LOOP. Modify it for the times you want to plot! It depends on your file. If you have yearly files, t=1 will be January 1st at 00Z, t=2 will be the same day but at 06Z and so on. Here it's just plotting the first 4 times of the file - one day worth of data.
for i=1:4
data=nc_varget(file,var,[i z 0 0], [1 1 73 144]);
%
% Create map and name it. Set latitudes on first brackets. Longitudes on second. Naming the map is useful for manipulating its settings later on. You can change the values of lat / lon for whatever suits you!
%
mapa = worldmap([-90 90],[0 360]);
%
% Plot data for the current time.
%
contourfm(lat,lon,data,205:5:300);
%
% Insert boundaries.
%
geoshow('landareas.shp','Facecolor','none');
%
% Define color and insert colorbar. The first number is the interval between each value. Give a name to this as you'll need to edit the labels on the colorbar to make it look nicer.
%
barra=contourcmap(5,'jet','colorbar','on','location','horizontal','coloralignment','center');
%
% Set map grid lines to black color (k is short for black)
%
gridm('gcolor','k');
%
% This will create black axis with a white background and a box. Set them off. GCA means you're dealing with the graphic's current axis. It'd be gcf to deal with the current figure's properties.
%
set(gca,'color','none','box','off','visible','off');
%
% The colorbar will probably have way too many marks. Make it less! This command will make it only show the first, third, fifth, etc.
set(barra,'xtick',[1 3 5 7 9 11 13 15 17 19],'fontname','helvetica narrow');
%
% Get rid of the 90degree latitude label, or it'll be over the title and make things look ugly. This command sets parallel labels only to latitudes 45, 60 and 75. Now the name of the map comes in handy!
%
setm(mapa,'plabellocation',[45 60 75]);
set(mapa,'fontname','helvetica narrow');
%
% Define date and time of plot
%
% Jan 1st (merely an example!)
switch i
case 1
date='WRITE THE DATE CORRESPONDING TO t=1'
fname=(sprintf('%s_WRITE DATE TO BE STORED AS FILE NAME WHEN ITS SAVED_%d',var,nivel))
case 2
date='Jan 01 at 06UTC'
fname=(sprintf('%s_01_01_06_%d',var,nivel))
case 3
date='WRITE THE DATE CORRESPONDING TO t=3'
fname=(sprintf('%s_01_01_12_%d',var,nivel))
case 4
date='WRITE THE DATE CORRESPONDING TO t=4'
fname=(sprintf('%s_10_01_18_%d',var,nivel))
end
%
% Write title. Modify to name of variable that you're plotting.
%
title(sprintf('\\fontname{Times}Name of variable - %s',date),'fontsize',16)
%
% Save figure. The -depsc tells it to save in .eps format.
%
print('-depsc',fname)
%
% End the loop
%
i=i+1
end
disp('**************SUCCESSFUL COMPLETION OF PLOTS**************')
end
%
if numdims==3
%
%---------------------------------------------------------------------------------------------------
% LOOP IN case YOU HAVE A 3D FILE
%---------------------------------------------------------------------------------------------------
%
% Create one matrix per time and plot it (in this case from Oct 8th at 00Z through Oct 12th at 18Z).
%
for i=1:4
data=nc_varget(file,var,[i 0 0], [1 73 144]);
%
% Create map and name it. Set latitudes on first brackets. Longitudes on second. Naming the map is useful for manipulating its settings later on.
%
mapa = worldmap([-90 90],[0 360]);
%
% Plot data for the current time.
%
contourfm(lat,lon,data,205:5:300);
%
% Insert boundaries.
%
geoshow('landareas.shp','Facecolor','none');
%
% Define color and insert colorbar. The first number is the interval between each value. Give a name to this as you'll need to edit the labels on the colorbar to make it look nicer.
%
barra=contourcmap(5,'jet','colorbar','on','location','horizontal','coloralignment','center');
%
% Set map grid lines to black color (k is short for black)
%
gridm('gcolor','k');
%
% This will create black axis with a white background and a box. Set them off. GCA means you're dealing with the graphic's current axis. It'd be gcf to deal with the current figure's properties.
%
set(gca,'color','none','box','off','visible','off');
%
% The colorbar will probably have way too many marks. Make it less!
set(barra,'xtick',[1 3 5 7 9 11 13 15 17 19],'fontname','helvetica narrow');
%
% Get rid of the 90degree latitude label, or it'll be over the title and make things look ugly. This command sets parallel labels only to latitudes 45, 60 and 75. Now the name of the map comes in handy!
%
setm(mapa,'plabellocation',[45 60 75]);
set(mapa,'fontname','helvetica narrow');
%
% Define date and time of plot. %s puts the variable name, defined in the beginning of the script, to the filename fname!
%
switch i
case 1
date='MODIFY AS SUITS YOU'
fname=(sprintf('%s_DATE_sfc',var))
case 2
date='MODIFY AS SUITS YOU'
fname=(sprintf('%s_DATE_sfc',var))
case 3
date='MODIFY AS SUITS YOU'
fname=(sprintf('%s_DATE_sfc',var))
case 4
date='MODIFY AS SUITS YOU'
fname=(sprintf('%s_DATE_sfc',var))
end
%
% Write title
%
title(sprintf('\\fontname{Times}Name of variable - %s',date),'fontsize',16)
%
% Save figure. The -depsc tells it to save in .eps format.
%
print('-depsc',fname)
%
% End the loop
%
i=i+1
end
disp('**************SUCCESSFUL COMPLETION OF PLOTS**************')
end
