June 12, 2000

You can contact me at 303-497-1708 or by emailing thoar@ucar.edu

If you are at NCAR, my extension is simply 1708.

Matlab is an interpreted language that handles matrices with a
natural syntax. This generally means a **very** short spin-up
period for people unfamiliar with Matlab. One of Matlab's founding fathers
(who is still president of the MathWorks) is Cleve Moler of Linpack fame.

As with any interpreted language, it is possible to write some incredibly slow F77 code. We recommend against that ... Take advantage of the vector- or matrix- based syntax and avoid loops whenever possible.

- Preparing your environment
- Starting Matlab
- Help
- Syntax
- Public Tools
- Example 1 Import ascii, make a scatterplot.
- Example 2 Make a contour plot.
- Example 3 Meridional/Zonal averages
- Example 4 Linear Least Squares
- Example 5 Masking
- Example 6 Masked Contour Plot
- Example 7 Cheap Projections
- Scripts
- Functions
- Example 8 Profiling
- Example 9 Handle Graphics
- Printing
- Example 10 Importing binary files
- Example 11 netCDF
- Movies movies/animations

Nothing is really required, but there are some things that make life nicer, naturally.

mkdir ~/matlab mkdir /local/scr/data/ mkdir /local/scr/ export EDITOR=vi

Matlab searches for a **~/matlab** directory
upon startup, this is a nice place to have general purpose scripts/functions.
If there exists a **~/matlab/startup.m** file, it is automatically
executed upon startup. The file may contain any sequence of valid Matlab
commands. If there are directories that contain certain third-party
toolboxes, the obvious way to add them is with your **startup.m** file,
which could be as simple as this:

% Startup.m is executed by the matlabrc script automatically at startup. % path(path,'/local/scr/data'); % append these directories to path

Simply type **matlab** at the UNIX prompt. If your PATH is correct,
this will work. If it is not, you might as well give up now. You can fire up
Matlab from any directory and your MATLABPATH will be set to search
the normal Matlab installation, anything specified in your
**~/matlab/startup.m** file, and the current working directory.

The matlab prompt is:

>>if you didn't get the little splash screen, it means your

>> quit

There are lots of ways to get help in Matlab. Simply typing **help**
at the command line gets you started. If you know the command you want,
(for example "fft") simply type **help fft**.

If you do not know the name of the command, you can look for functions
by keywords -- **lookfor fft**, for example. This works for ALL functions
in the MATLABPATH, including the ones you create. More on this when we get
into the function section, ditto for "help".

Matlab also has an interactive WWW-based help system, called HelpDesk.
This is the 800# gorilla of help. Search all the manuals, connect to the
MathWorks FAQ site, list by topics, alphabetically, etc.
It will commandeer your browser, simply by typing **helpdesk**.
"If you don't have one; GET ONE!"

Matlab is short for *Mat*rix *Lab*oratory. If you can
think of the algebra in 2D, it will probably work exactly as you expect.
The thing to get used to is that any command without a ";" at the end
will print the result of the command to the screen. Not particularly
desirable for large results. I am in the habit of terminating **every**
command with a ";", whether it needs one or not.

General help for all the relational operators can be had with **help
relops**. The only tricky part is that the square brackets are used for
defining arrays, parentheses are for indexing, and the variables MUST be
conformable.

a = [0:127]; % "a" is a row vector b = cos(2*pi*a/length(a)); % so is "b" plot(b(1:64)); % only plot the first 64 elements title('My mental retention prowess.'); xlabel('(minutes)'); ylabel('fraction of capacity');

Multiplying a 10x3 matrix by a 10x3 matrix will generate an error message, unless you want the operator to work on an element-by-element basis. For example, if you want to multiply each element of an array by the corresponding element of another array:

c = a.*b; % should be reminiscent of the "dot product" d = sum(c); % wow

It is best to simply ignore the existence of vectors and just admit
they are really N x 1 or 1 x N matrices. All the rules of matrix algebra
apply. Some functions are tailored for 1D matrices and can be used on
higher-dimensional matrices, like **sum**.

c = b'*b; % premultiplying b by its transpose --> 128 x 128 d = sum(c) % what is "d"? (semicolon intentionally left off) imagesc(c); % color-code the result and plot colorbar; % plot a legend

The classic one is **fft**. Taking the fft of a 1D variable is pretty
clear. What happens if you feed **fft** a M x N matrix? Matlab grew from
Fortran beginnings, so it takes ffts of each **column**,
since these are stored sequentially in memory (and it can take unit strides).
If you want a real 2D fft, use **fft2**. This is a common theme in Matlab.
Operations that have a different connotation for 2D as opposed to 1D have an
additional "2" in the function name.

load topo; wow = topo - mean(topo(:)); C = fft(wow); % take the fft of each column of "c" whos % "who and size" (semicolon intentionally left off) D = fft2(wow); subplot(2,2,1) imagesc(real(C)) title('Real part of FFT of each column') colorbar subplot(2,2,2) imagesc(fftshift(real(D))) title('Real part of 2DFFT') colorbar subplot(2,2,3) imagesc(imag(C)) title('Imaginary part of FFT of each column') colorbar subplot(2,2,4) imagesc(fftshift(imag(D))) title('Imaginary part of 2DFFT') colorbar

There are simply hundreds of third-party books and toolboxes. Check out the MathWorks home page and poke around.

first, we copy the ascii file to
/local/scr/data/soi.ascii and get on with it. The matlab variable will be
the filename without the extension, so we should wind up with
a variable called **soi**.

load /local/scr/data/soi.ascii whos; % the command you will use most years = soi(:,1); % save years as Nx1 matrix a = soi(:,2:13); % strip out years, other columns are monthly values. inds = find(a < -900); % locate the "missing" values a(inds) = NaN; % replace with Matlab's "missing" flag. [nrows,ncols] = size(a); b = reshape(a',nrows*ncols,1); % reformat into Nx1 matrix. plot(b) % plot, but pretty ugly % % Lets zoom in on just the first 10 years (120 months) % axis([0 10*12 -Inf Inf]) % [xmin xmax ymin ymax] (Inf == unknown) grid % % It would be convenient to have a better time axis. % I will make a matrix the same size as the input with each entry % corresponding to the month midpoint == year.fraction % t = ones(nrows,1) * [0.5/12 : 1/12 : 11.5/12]; % matrix of fractions tmat = years * ones(1,ncols) + t; % matrix of times T = reshape(tmat',nrows*ncols,1); % array of times plot(T,b)

The only thing tricky about contour plots in matlab is the resulting orientation. (1,1) is at the bottom left.

clear % clears all variables in your workspace!!!!! load topo % a default dataset whos contour(topo) colorbar

Changing the number of contour lines is easy.

contour(topo,5) colorbar

Specifying which contour levels is also easy.

contour(topo,[0:500:4500]) colorbar

Personally, I like the imageplots a lot better, they give clues as to the original resolution. However, the plot has (1,1) in the upper left.

imagesc(topo) colorbarThis is the introduction to "Handle Graphics". It is generally trivial to tweak almost any aspect of a figure.

imagesc(topo) set(gca,'YDir','normal') % GetCurrentAxis attribute "YDir"; set to ... colorbarTo see _what_ you can tweak, do a:

get(gca)There are other neat plots, like

[x,y,z] = peaks(360); % CREATE SOME FAKE DATA (z is 360x360) imagesc(z); % TAKE A QUICK LOOK AT THE DATA

newz = z(1:2:360,:); % GRAB EVERY OTHER ROW OF DATA imagesc(newz); % IMAGE THE SUBSETTED MATRIX set(gca,'YDir','normal'); % PLOT 1,1 IN THE LOWER LEFTTaking the mean is simple. In general, Matlab applies operators to each column, there is no need to "loop" over them. Since we conceptualize our matrix

clf; % CLEAR OUT THE FIGURE WINDOW (GOOD STYLE) subplot(2,2,1); % CHOP PLOT WINDOW INTO A 2X2 MATRIX imagesc(newz); % CHEAP PLOT OF DATA set(gca,'YDir','normal'); % PLOT 1,1 IN THE LOWER LEFT title('Peak data'); subplot(2,2,4); % use the 4th window of the 2x2 contour(newz); % duhhhhhhhh title('Contoured peak data'); subplot(2,2,3); % USE THE THIRD WINDOW OF THE 2x2 set mmeans = mean(newz); % CREATE MEANS OF EACH COLUMN (ALL LATITUDES) plot(mmeans); % PLOT THEM (versus their index [1,360]) title('Meridional means'); subplot(2,2,2); zmeans = mean(newz'); % TRANSPOSE THE MATRIX AND FIND THE MEANS plot(zmeans); % PLOT THEM (versus their index [1,180] title('Zonal means');

[cs,h] = contour(peaks); % this is kinda sneaky clabel(cs,h,'fontsize',15,'color','r','rotation',0)

clf; % CLEAR OUT THE FIGURE WINDOW (GOOD STYLE) a = hamming(length(zmeans)); % CREATE WEIGHTING FUNCTION (column vector) b = zmeans'; % CHANGE ROW VECTOR TO COLUMN VECTOR c = a.*b; % CREATE THE WEIGHTED AVERAGE lats = [-89:1:90]; % CREATE A "REALISTIC" LATITUDE AXIS plot(lats,zmeans,'-.',lats,c,'-'); grid; % PUTS A GRID ON THE AXIS legend('unweighted zonal mean','weighted zonal mean')Here's where you click/drag the legend (multiple times if you wish). Instead of

The solution to **A*X = B** is simply **X = A\B;**
Pretty boring, but true.

help \or

help relop

Just for grins, put worldmap.m in your ~/matlab directory.

clear; clf; % clear workspace, clear figure load topo; % get some elevation data [180x360] lon = [0.5:359.5]; lat = [-89.5:89.5]; a = peaks(360); % make some bogus data (too big) cph = a(1:2:360,:); % cph is same size as elevation data mask = topo < 0; % make a land/sea mask (sea is 1.0) clear a topolegend topomap1 topomap2; % for illustrative purposes figure(1); clf; subplot(3,1,1); imagesc(lon,lat,topo); set(gca,'YDir','normal'); colorbar; subplot(3,1,2); imagesc(lon,lat,mask); set(gca,'YDir','normal'); colorbar; colormap(jet); subplot(3,1,3); imagesc(lon,lat,cph); set(gca,'YDir','normal'); colorbar; figure; % GET ANOTHER GRAPHICS WINDOW newdat = cph.*mask; % Use only "cph" data over land imagesc(lon,lat,newdat,[-7 7]); % quick look set(gca,'YDir','normal'); % correct orientation title('"Peaks" Masked by the Earth''s topography.'); colorbar; % because I said so tim = jet; % snag a colormap tim(32:33,:) = 0.9; % make colors around 0 == gray colormap(tim); % apply new colormap worldmap

Just for grins, put sinusoid.m in your ~/matlab directory.

figure(1); clf a = sinusoid(360,1,1,90); % make some phony data b = a*a'; c = b(1:2:360,:); subplot(2,1,1) imagesc(lon,lat,c); % look at the phony data set(gca,'YDir','normal'); % correct orientation % MASK OUT THE "LAND" i = find(mask > 0); % Find locations of land c(i) = NaN; % set land elements to NAN subplot(2,1,2) contour(lon,lat,c,[-10:2:10]); % "gapped" contours hold on; % add to current plot [d,h] = contour(lon,lat,topo,[0 0]); % continents

clf; load topo [x,y,z] = sphere(35); % MAKE A 3D SURFACE h = surface(x,y,z,'FaceColor','texture','Cdata',topo) % DO TEXTURE MAP, % KEEP GRAPHICS HANDLE set(h,'EdgeColor','none') % USE HANDLE TO CHANGE % AN ATTRIBUTE axis square; % MAKE AXES EQUAL axis off; % TURN OFF LABELLING colormap(jet); % A DIFFERENT PALETTE

Also try the free mapping routines in the SEA-MAT distribution.

By now, you should be sick of cutting and pasting. Enter the script,
generically called a ".m" file. By inserting your matlab commands
into a file called (for example) **foo.m**, you can simply type
"foo" at the Matlab prompt and the entire file is checked for syntax
and executed. Put the following lines in a file called
**~/matlab/foo.m**:

% This is a generic matlab script. % load topo; % GET [180x360] ELEVATION DATASET lats = [-89.5:89.5]; % CREATE LAT ARRAY FOR TOPO MATRIX lons = [0.5:359.5]; % CREATE LON ARRAY FOR TOPO MATRIX imagesc(lons,lats,topo); % CREATE SOME PLOT w/ true x,y limits set(gca,'YDir','normal'); % CORRECT ORIENTATION worldmap;

Then, (one at a time) from Matlab's command line:

which foo help foo type foo fooNotice that the initial block of comment lines forms the "help" entry for the script! Pretty cool. A nice way to keep documentation up-to-date with the program.

A Matlab Function is a special .m file. The first line of the file is the function declaration.

You can peruse most functions with **type**.

The FIRST line of the help file should contain some keywords.
This is the only line searched by the **lookfor** function.

profile your favorite script/function --- VERY impressive. follow the example from "help profile"

Perhaps the best part of Matlab. Infinitely customizable graphics.

You will definitely want wysiwyg.m in your ~/matlab directory.

% A pretty ugly figure h1 = subplot('position',[0.2 0.4 0.4 0.2]); % l, b, w, h plot(1:10) xlabel('wow') ylabel('unbelievable!') h2 = subplot('position',[0.1 0.1 0.8 0.2]); plot(11:20) xlabel('m/s^2') ylabel('t_i') h3 = subplot('position',[0.2 0.6 0.6 0.3]); imagesc(peaks) title('trumpets blare, twice') h = colorbar; orient tall % fill a page -- portrait style wysiwyg % What You See Is What You Get ... wysiwyg % Clean it up. get(h) set(h,'Position',[0.743 0.6 0.025 0.3]); % thinner colorbar child = get(h,'Ylabel'); % get handle to the Ylabel. set(child,'String','meters/furlong') get(h3) set(h3,'YTick',[1 13 22 38]); % change tick positions set(h3,'YTickLabel',[' 1';' 13';'kahuna';' 38']); % change labels child = get(h3,'Title') set(child,'String',{'Trumpets Blare','Two times'},'FontSize',18) get(h1) set(h1,'Position',[0.2 0.4 0.6 0.1], ... 'YAxisLocation','right'); % line continuation, multiple % attributes in one "set"

Easy, as long as you know how the data is written.
For this exercise, save
some binary data as
**/local/scr/data/chi_jan96.ieee**

With a clever combination of **fread**, I feel supremely confident in
the ability to read any valid file. Read the help file for **fread** and
**fopen**

The binary file below was written in the following manner:

real datmat(129,64) integer nt do i = 1,nt ... write(iunit)datmat enddoso each record is (1+129*64+1)*4 bytes, the file is this_same_number x nt. Logically the format is:

[4bytes][129*64*4 bytes of data][4bytes] [4bytes][129*64*4 bytes of data][4bytes] ... [4bytes][129*64*4 bytes of data][4bytes] [4bytes][129*64*4 bytes of data][4bytes]The following segment just reads the data and makes a crude contour plot. For more information on contour plots -- try

nlon = 129; nlat = 64; nt = 5; fid = fopen('/local/scr/data/chi_jan96.ieee','r') frewind(fid); for islice = 1:nt, expr = sprintf('slice %.0f',islice); dum1 = fread(fid, 1,'float32'); udat = fread(fid,[nlon nlat],'float32'); dum1 = fread(fid, 1,'float32'); contour(udat); title(expr); disp('Hit a key to continue ...'); pause; end

Basically exactly the same as the unformatted binary without the little 4 byte record headers/enders. Every record is constrained to be the same, however. The logical format is:

[xxx bytes of data] [xxx bytes of data] ... [xxx bytes of data] [xxx bytes of data]Your Matlab script to read it could be:

fid = fopen('/local/scr/data/DirectAccess.ieee','r') frewind(fid); for islice = 1:nt, udat = fread(fid,[xxx],'float32'); end

The most flexible yet. "Let the buyer beware". C binary generally comes with no header bytes, everything is simply data pasted together. If you get off by a byte, you get garbage. There is no error-checking. You can mix reads of different lengths and data types in any fashion. The following is snytactically correct, but may not match the data at all.

fid = fopen('/local/scr/data/wildcard.ieee','r') udat = fread(fid,[10],'float32'); a = fread(fid,[10 30],'int32'); b = fread(fid,[1],'float64');

- Color Postscript:

print -dpsc -P[printer]

- Black & White Postscript:

print -dps -P[printer]

- Color Postscript:

print -dpsc [filename]

- Black & White Postscript:

print -dps [filename]

Lots of postscript files are simply ascii files of postscript commands. You can view and change the source with your favorite editor. This shifts the burden of rendering the figure from the computer to the printer. There can be some relatively small postscript files that take an extremely long time to print because of this. By making the computer render the figures to some finite raster density, the computational burden stays on the computer at the expense of making potentially large (but simple and relatively fast to print) postscript files for the printer. This is particularly true for maps.

If you have figures that are taking too long to print, you may consider having Matlab do the rendering. There is a huge tradeoff between size of the file and the raster level (dots per inch).

print -dpsc2 -zbuffer -r200 [filename]

NetCDF [overview] [history] [documentation] [reading] [creating] [low-level examples] [KeithLindsays examples]

path(path,'/contrib/matlab'); ncstartupOr you can beat the rush by simply putting those two lines in your

- CSIRO matlab/netCDF
interface (pay particular attention to the
`getnc`function) useful for reading netCDF files (in my opinion), - the netCDF Toolbox users guide the first thing you should read,
- the netCDF Toolbox Home Page (a bit broader than the users guide),
- and the "help" facility in Matlab -- particularly
`help netcdf`.

fname = '/contrib/matlab/ncexample.nc'; % saves typing depth = getnc(fname,'depth'); % depth is a *gasp* matrix

cdf = '/data/csm/b003.01/atm/01-0011.to.0060.nc'; [cdfid,rcode ]=mexcdf('open',cdf,'NOWRITE'); [varid]=mexcdf('varid',cdfid,var) [var_name,var_type,nvdims,var_dim,natts]=mexcdf('varinq',cdfid,varid) mexcdf('close',cdfid);

I am still building this ... as of May 1, 2002

NetCDF [overview] [history] [documentation] [reading] [creating] [examples]

ftp://ftp.mathworks.com/pub/contrib/v5/graphics/mpgwrite

ftp://ftp.mathworks.com/pub/contrib/v5/graphics/mpgread

Tim Hoar