# Matlab/Octave tutorial

This is a basic introduction to Matlab/Octave.

## Numbers

• Use a dot (.) as decimal separator (e.g. 3.14 and not 3,14)
• You can use the scientific notation $$a \times 10^b$$ using the e-notation. (e.g. $$3 \times 10^{-7}$$ becomes 3e-7)

### Useful constants

Various constants are also pre-defined: pi, e (Euler's number), i and j (the imaginary number), inf (Infinity, result from e.g. 1/0) and NaN (Not a Number - result from e.g. 0/0)

NaN is very special: NaN+ any number is NaN

## Variables

• Numbers (and any other data type) can be put into variables
• The value of a variable is referenced by its name
• A variable name can be composed by letters (a-z and A-Z, no accents), numbers (0-9) and underscore (_). The first character must be a letter.
• Example:
         temp = 21;

The variable temp has now the value 21. The value of the variable can be changed later on.
• An assignment without a final semicolon echos its value to the screen
• Careful: the value of constants can be overwritten. The following is allowed but not encouraged:
         pi = 3;

• The command whos lists the currently defined variables and their size.

## Strings

• Delimited by single quotes
       'Hello world'

• How to use a single quote in a string? -> two consecutive single quotes become one single quote. For example

'That''s great'


## Vectors and matrices

• Vectors are a list of numbers
• Matlab/Octave distinguishes between :
• a row vector: [1 4 5 3]
• a column vector: [1; 4; 5; 3]
• Matrices are tables of numbers. Rows are separated by a semicolon

[1 2; 3 4]

• Consecutive elements can be written as
         [first:step:last]

or simply, if the step is 1,
         [first:last]

• For example [1:5] is the same as [1 2 3 4 5]
• For example [1:2:6] is the same as [1 3 5]

### Indexing

• Individual elements of a vector or matrix can be addressed by their index
• The second element of a vector a is for example a(2)
• The element at the second row and the first column of a matrix A is for example A(2,1)
• The special word end refers to the last index.
• One can also use a list of indexes to extract a part of the vector or matrix
• The symbol colon : is a short-hand for 1:end
• If a matrix is indexed with only one subscript, the matrix is treaded as a vector where all columns are concatenated
• Matlab/Octave supports also higher-dimensional arrays

#### Example

Assume the following matrix A

         A = [ 1  2  3; ...
4  5  6; ...
7  8  9; ...
10 11 12];


What is A(1,1)

1

A(1,1)

What is A(end,1)

10

A(end,1)

What is A(1:2,1:2)

[1 2; 4 5]

A(1:2,1:2)

What is A(2,:)

is it a row or column vector?

[4 5 6]

A(1:2,1:2)

What is A(:,2)

is it a row or column vector?

[2; 5; 8; 11]

A(1:2,1:2)

### Operators

• Scalar and matrix operations: + sum, - difference, * multiplication, / division
• Element-wise matrix operations: .* multiply element-wise, ./ divide element-wise

Try the matrix multiplication and the element-wise multiplication of the matrices [1 0; 0 1] and [1 2; 3 4]. Enter the difference (multiplication minus the element-wise multiplication) below.

use * and .*

[0 2; 3 0]

A = [1 0; 0 1];
B =  [1 2; 3 4];
A*B - A.*B


### Comparison operators

• element-wise equal (==), element-wise different (~=)
• logical "and" (&&) logical "or" (||) (with short-circuit evaluation)
• logical element-wise "and" (&) logical element-wise "or" (|)
• The results of such operators can also be used to index an array
• For example return all elements in the variable T which are greater than 10 but less than 20.
T(10 < T & T < 20)

• The function find(condition) returns the indexes of all elements where the condition is true.
• false is zero and true is 1. For instance to count the number of elements in the vector T which are larger than 20 one can use sum(T > 20).

Using the matrix A, compute the sum of all elements larger than 5

63

sum(A(A(:) > 5))
% or
sum(sum(A(A > 5)))


### Useful functions

sin, cos, tan
trigonometric functions
asin, acos, atan
inverse trigonometric functions
log, log2, log10
natural, base 2 and base 10 logarithms:
exp
exponentiation
abs
absolute value
sqrt
square root
mean
mean
median
median
std
standard deviation
var
variance
isnan
Check if variable is NaN. Note that NaN == NaN is false!
inv
inverse of a matrix
sum
sum of all elements
prod
product of all elements
max,min
maximum,minimum value

These function can also operate of a given dimension: sum(array,dimension)

Find out more of these function by typing help function name.

What is the product of all integer numbers between 1 and 10 (included)?

use prod

3628800

prod(1:10)

What is the product of all even integer numbers between 1 and 10 (included)?

use prod

3840

prod(2:2:10)

What is the maximum in row 2 of matrix A given above?

use max

6

max(A(2,:))

What is the minimum in column 3 of matrix A given above?

use min

3

min(A(:,3))

What is the general maximum of matrix A given above?

use max

12

max(A(:))

What is the sum of A over the second dimension?

use sum

[6; 15; 24; 33]

sum(A,2)

What is the length of the hypotenuse of a triangle with sides of length 4 and 7?

use this formula from the Pythagorean theorem

8.0623

sqrt(4^2 + 7^2)

To find new commands try lookfor (or your favorite search engine).

What function would you use to fit a polynomial to your data?

A may try to Google "matlab fit a polynomial"

polyfit

lookfor fit

## The file system

• On every current operating system, files are organized in a tree of directories starting from a root directory
• The absolute path of a directory or file defines which directories to follow starting from the root directory to the given directory or file
• In Linux/UNIX/Max OS X, files and directory names are separated by a slash (/), on windows by a backslash (\)
• In order to avoid to deal with long path names, every program has a current working directory
• The current working directory from Matlab/Octave can be queried with the command pwd.
• The relative path of a directory or file defines which directories to follow starting from the current directory to the given directory or file
• In relative path, two dots (..) represent the parent directory.
• To change the current directory, you can use the command cd('new/path')

Matlab/Octave is currently in the directory /home/user/work. You need to access a file in /home/user/data/file1.txt. What would be the relative path to this file on a Linux system? (shortest possible form)

You may need to use ..

../data/file1.txt

## Importing/Exporting data

### ASCII format

• To read ASCII data in Matlab/Octave, tables should be saved as an ASCII text file using space or tab as element separator. Each line corresponds to one row. Make sure that a dot is used as a decimal separator.
• Load a file named data.txt, creates a variable called data  load data.txt 
• or to load the file in a variable called var:  var = load('data.txt'); 
• Saving the variable data in the file data.txt using the ASCII format  save -ascii data.txt data 

### NetCDF format

• Reading a variable called var data from a NetCDF file  data = ncread('file.nc','var'); 
• Writing a variable called var data to a NetCDF file  ncwrite('file.nc','var',data); 

## Scripts

• A series of commands can be collected in a script file
• A script file has the extension .m
• How can Matlab/Octave find your script file?
• it must be either in your current work directory
• the directory containing the script file must be added to the search path using addpath. For example
  addpath('/path/to/directory');


## Functions

• Functions are similar to scripts
• Unlike scripts, functions can have input/output parameters
• For example a function calculating the speed of ocean current based on the zonal and meridional component
function speed = current_speed(u,v)
speed = sqrt(u^2 + v^2);

• File name and function name should match

## Dates

• Matlab/Octave unit for time is days since the January 0, 0000 (or rather December 31, -1) at 00:00 UTC.
• You obtain this serial date number with the function datenum(year,month,day,hours,minutes,sec)
• To convert back to year, month, day, hours, minutes and seconds you can use datevec.
• A nicely formatted string can be obtained by datestr
• On figures, you can use datetick for nice labels

How many days are between 01-Jan-2013 and the 04-Apr-2013

use max

93

datenum(2013,4,4)- datenum(2013,1,1)


## Plotting

### One-dimensional data (e.g. time series, vertical profiles,...)

plot(x,y,format)
Draws a line with the values in x and in y as x- (horizontal) and y-axis (vertical) respectively. With format one can specify the color (blue ('b'), red ('r'), green ('g'), ...) and style of the line (solid (-), dots (.), dotted (:), ...)

Download these data (sea level time series in the West Florida Shelf, in the 6th column) and make a plot, with a solid line in green. The date can be derived from the 5 first columns, using the command datenum. Include labels with the variable units (meters) and the date, and add a legend. Your image should look like this: Did it work?

you may answer yes or no

yes

data = load('8762075.sealevel.txt');
date = datenum(data(:,1),data(:,2),data(:,3),data(:,4),data(:,5),0);
plot(date,data(:,6),'g')
datetick('x',19)
ylabel('m','rotation',0)
legend('sea level')


At which date did the sea level reach its maximum? Enter the day in the format yyyy-mm-dd HH:MM? You may want to use datestr.

Remember that you can you can index an array by a condition or you may use the function find.

2004-10-10 13:00

data = load('8762075.sealevel.txt');
date = datenum(data(:,1),data(:,2),data(:,3),data(:,4),data(:,5),0);
max_sealevel = max(data(:,6));
max_data = date(max_sealevel == data(:,6));
datestr(data_max,'yyyy-mm-dd HH:MM')


### Two-dimensional data (e.g. horizontal sections, ...)

pcolor(x,y,v)
The value within a rectangle defined by x and y is drawn by color depending on v and on the color map. Per default pcolor plots also the grid line which can result in a completely black image if the number of rows and columns is very high. This can be avoided with the command shading flat
colormap(map)

Change the color map. Possible values for map are jet, hsv, hot,...

colorbar
show color bar relating values and colors

title('my figure')
give a title to the current figure
xlabel('my label'), ylabel('my label')
give a name to the x- and y-axis
print -dpng file.png
Save the figure as a PNG file. For a EPS file use -depsc. Do not save images in JPEG as it degrades the quality

## Loops

• Loops have a counter which takes successively all elements of a row vector
for i = [1 2 10 20]
i
end

• Loops are often used with a range of values
for i = 1:100
i
end

• Loops are slow and can often be avoided.
• For example, sum all integer from 1 to 10
total = 0;
for i = 1:10
total = total + i;
end
total

• Can simply be computed as sum(1:10).

## if-statement

• Sometimes your code needs to behave differently depending on some conditions.
• if-statement has the following structure.
if some_conditions
% do something
else
% do something else
end

• The else section can be omitted.
• For example.
if x < 0
x = -x;
end


Which Matlab/Octave function implements the last code example?

It was seen in useful function

abs

## Large data set

 >> ncdisp('WesternMedSST.nc')
Source:
WesternMedSST.nc
Format:
classic
Global Attributes:
Conventions = "CF-1.5"
title       = "DINEOF analysis of Western Mediterranean sea surface temperature"
author      = "Igor Tomazic"
institution = "GHER, University of Liege"
source      = "http://gher-dineof01.phys.ulg.ac.be:8080/opendap/DW/dineof/SVRI_L3C_SST1H/medwest2_005_327x217/fc_cv5_k2c/20130801_20130816/df_al0_0008_mc0_nk56_nm40/IO/data.nc.html"
history     = "DINEOF analysis"
references  = "DINEOF http://modb.oce.ulg.ac.be/mediawiki/index.php/DINEOF"
Dimensions:
lon  = 327
lat  = 217
time = 384
Variables:
lon
Size:       327
Dimensions: lon
Datatype:   double
Attributes:
standard_name = "longitude"
units         = "degree_east"
lat
Size:       217
Dimensions: lat
Datatype:   double
Attributes:
standard_name = "latitude"
units         = "degree_north"
time
Size:       384
Dimensions: time
Datatype:   double
Attributes:
standard_name = "latitude"
units         = "days since 1900-00-00 00:00:00"
seviri_sst
Size:       327x217x384
Dimensions: lon,lat,time
Datatype:   single
Attributes:
standard_name = "sea_water_temperature"
long name     = "sea surface temperature"
units         = "degree_Celsius"
_FillValue    = -9999
seviri_sst_filled
Size:       327x217x384
Dimensions: lon,lat,time
Datatype:   single
Attributes:
standard_name = "sea_water_temperature"
long name     = "sea surface temperature"
units         = "degree_Celsius"
_FillValue    = -9999
Size:       327x217
Dimensions: lon,lat
Datatype:   int8
Attributes:
long name = "land sea mask"
comment   = "1-sea; 0-land"



### Useful functions

ncdisp
display the content of a NetCDF file.
Read a variable from a NetCDF file.
ncwrite
Write a variable from a NetCDF file.
nccreate
Create a variable in a NetCDF file.
ncwriteatt
Write an attribute.

Load all variables in Octave/MATLAB. Use the variable names lon, lat, time, SST, SSTf and mask. Did it work?

you may answer yes or no

yes

lon = ncread('WesternMedSST.nc','lon');


What is the minimum value of longitude?

0.025

min(lon)


Check also the maximum and minimum value of longitude and latitude.

What is the spatial resolution of the data set?

0.05

lon(2)-lon(1) and lat(2)-lat(1)

What is the date and time of the first time instance in the data set? (in yyyy-mm-dd HH:MM:SS)

Use datestr and datenum.

2013-08-02 00:00:00

datestr(time(1) + datenum(1900,1,1),31)

Check also the 2nd and last time instance.

Plot the first time instance of the data set with pcolor. Your image should look like this:

you may answer yes or no

yes

pcolor(lon,lat,SST(:,:,1)'); shading flat, colorbar
print -dpng SST_first.png


What are the areas without color?

Plot the percentage of valid data grid point over time. Your image should look like this:

you may answer yes or no

yes

count = 100 * sum(~isnan(SST),3) / size(SST,3);
print -dpng ../html/images/SST_count.png


For all time instances, what is the percentage of sea grid points not covered by clouds? Look for the command squeeze. Your image should look like this:

you may answer yes or no

yes

count = squeeze(sum(sum(~isnan(SST),1),2));
percentage = 100 * count / sum(mask(:));
plot(time + datenum(1900,1,1),percentage)
datetick('x')
print -dpng ../html/images/SST_count_time.png


Plot the time average of SST. Your image should look like this:

Be aware than NaN + any number is NaN.

yes

SST2 = SST;
SST2(isnan(SST2)) = 0;
meanSST = sum(SST2,3) ./ sum(~isnan(SST),3);
print -dpng ../html/images/SST_time_mean.png

You could have also used nanmean which computes a mean ignoring NaNs:
SST2 = SST;
meanSST = nanmean(SST2,3);
print -dpng ../html/images/SST_time_mean.png


Plot the space average of SST (assuming that all pixels have the same area). Your image should look like this:

Be aware than NaN + any number is NaN.

yes

meanSSTt = sum(sum(SST2,1),2) ./ sum(sum(~isnan(SST),1),2);
meanSSTt = squeeze(meanSSTt);
plot(time + datenum(1900,1,1),meanSSTt)
datetick('x')
print -dpng ../html/images/SST_space_mean.png


On the figure of the previous exercise plot the result of nanmean(nanmean(SST,1),2). Your image should look like this:

yes

meanSSTt2 = nanmean(nanmean(SST,1),2);
meanSSTt2 = squeeze(meanSSTt2);
hold on
plot(time + datenum(1900,1,1),meanSSTt,'b')
plot(time + datenum(1900,1,1),meanSSTt2,'r')
hold off
legend('result with sum','result with nanmean');
datetick('x')
print -dpng ../html/images/SST_space_mean2.png


Why is there a difference?

Make a time serie with the number of pixels with the temperature larger than 25 degree Celsius. Your image should look like this:

yes

count = squeeze(sum(sum(SST > 25,1),2));
plot(time + datenum(1900,1,1),count)
datetick('x')
print -dpng ../html/images/SST_count25.png


Make a time serie of the area (in km2) with the temperature larger than 25 degree Celsius. You might want to use reshape and repmat. Your image should look like this:

The radius of the Earth is 6371 km.

yes

% Earth Radius (in km)
R = 6371;
% surface of each cell
dx = pi * 0.05 * R/180;
dy = pi * 0.05 * R/180 * cos(pi*lat/180);

dS = dx*dy;
dS2 = reshape(dS,[1 217 1]);
dS3 = repmat(dS2,[327 1 384]);
dS3(SST <= 25) = 0;
area = squeeze(sum(sum(dS3,1),2));
plot(time + datenum(1900,1,1),area)
datetick('x')
print -dpng ../html/images/SST_area25.png