Matlab/Octave tutorial

Alexander Barth, Aida Alvera-Azcárate

http://gher-diva.phys.ulg.ac.be/MatlabOctaveTutorial/

This is a basic introduction to Matlab/Octave.

Numbers

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

Strings

Vectors and matrices

Indexing

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

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

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

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

NetCDF format

Scripts

Functions

Dates

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: Example plot figure 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,...

color maps

colorbar
show color bar relating values and colors

Annotating your graphs

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

if-statement

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
    mask
           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.
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: Example plot figure

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: Example plot figure

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: Example plot figure

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: Example plot figure

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: Example plot figure

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: Example plot figure

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: Example plot figure

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: Example plot figure

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