This is a basic introduction to Matlab/Octave.

- 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)

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`

- 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.

- 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 are a list of numbers
- Matlab/Octave distinguishes between :
- a row vector:
`[1 4 5 3]`

- a column vector:
`[1; 4; 5; 3]`

- a row vector:
- 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]`

- 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

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)

- 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

- 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)))

- 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

- 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

- 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`

- 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);`

- 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 exampleaddpath('/path/to/directory');

- 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

- 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)

`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')

`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

- Let your computer do repetitive tasks!
- 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

- DINEOF analysis of Western Mediterranean sea surface temperature. Data is also available in mat format.
- Data in NetCDF
- Download file and view content of NetCDF file with
`ncdisp('filename.nc')`

>> 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 mask Size: 327x217 Dimensions: lon,lat Datatype: int8 Attributes: long name = "land sea mask" comment = "1-sea; 0-land"

- ncdisp
- display the content of a NetCDF file.
- ncread
- Read a variable from a NetCDF file.
- ncwrite
- Write a variable from a NetCDF file.
- nccreate
- Create a variable in a NetCDF file.
- ncreadatt
- Read an attribute.
- 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'); lat = ncread('WesternMedSST.nc','lat'); time = ncread('WesternMedSST.nc','time'); SST = ncread('WesternMedSST.nc','seviri_sst'); SSTf = ncread('WesternMedSST.nc','seviri_sst_filled'); mask = ncread('WesternMedSST.nc','mask');

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); pcolor(lon,lat,count'), shading flat, colorbar 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); pcolor(lon,lat,meanSST'), shading flat, colorbar 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); pcolor(lon,lat,meanSST'), shading flat, colorbar 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 km^{2}) 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