Live Script edama_01
This Live Script supports Chapter 1 of Environmental Data Analysis with MATLAB or PYTHON, Third Edition, by WIlliam Menke
The script naming system has the form edama_CC_SS, where CC is the chapter number and SS is the script number.
% edama_01_01: print the date
date
% edama_01_02: directory manipulation
clear all; % clear all previously-defined variables
% routinely doing a clear at the top a LiveScript is a good idea
% print working directory
pwd
% the presumption is that the current directory
% is edama/LiveScripts, so change directory to one up
% and then down into TesFolder
cd ../TestFolder
pwd
% list files and sub-folders in the current directory
dir
% change directory back to ch01
cd ../LiveScripts
pwd
% edama_01_03: simple algebra
% example of simple algebra, c=a+b with a=3.5 and b=4.1
a=3.5;
b=4.1;
c=a+b;
c
% edama_01_04: algebra with square root and power
% example of simple algebra, c=sqrt(a^2+b^2); with a=3 and b=4
a=3;
b=4;
c=sqrt(a^2+b^2);
c
% edama_01_05: algebra with sine
% example of simple algebra, c=sin(n*pi*(x-x0)/L) with n=2, x=3, x0=1; L=5
n=2; x=3; x0=1; L=5;
c=sin(n*pi*(x-x0)/L);
c
% edama_01_06: cell array
% this section mainly is provided to match the Python
% section, which is on list-like data objects.
% In MATLAB, the functionality of a list is provided
% by a data type called a "cell array". A "cell" can hold
% any kind of data; furthermore, cells can be grouped
% together into "cell arrays"
% Her ewe create a cellarray containing four items, in this case just
% numeric data
mycellarray = {1, 2, 4, 8};
mycellarray
% determine the length of the cell array
N=length(mycellarray);
N
% set the variable "a" to the contents of cell 3
a=mycellarray{3};
a
% cell arrays are not especially useful for storing numeric
% data but are much more useful in storing more complicated
% data items, and expecially, data items of variable size
% edama_01_07: vectors of zeros and ones
% a length-4 column-vector of zeros
c0 = zeros(4,1);
c0
% a length-4 column-vector of ones
c1 = ones(4,1);
c1
% a length-6 row-vector of zeros
r0 = zeros(1,6);
r0
% a length-6 row-vector of ones
r1 = ones(1,6);
r1
% 3x3 square matrix A0 of zeros, and
% 3x3 square matrix A1 of ones
A0 = zeros(3,3);
A0
A1 = ones(3,3);
A1
% length of column vector c1
N = length(c1);
N
% size of matrix A1
[N,M] = size(A0);
N
M
% edama_01_08: defining vectors and matrices
% define row-vector, column-vector and a matrix
r = [2, 4, 6];
c = [1, 2, 3]';
M = [ [1, 4, 7]', [2, 5, 8]', [3, 6, 9]'];
% print
r
c
M
% this works as well
M = [1, 2, 3; 4, 5, 6; 7, 8, 9];
M
% edama_01_09: vector and matrix addition
% column vector c
a = [1.,3.,5.]';
a
% column vector b
b = [2.,4.,6.]';
b
c = a+b;
c
% square matrix A
A = [1., 2., 3.; 4., 5., 6.; 7., 8., 9. ];
A
% square matrix B
B = ones(3,3);
B
C = A+B;
C
% edama_01_10: vector and matrix multiplication
% define matrices and vectors
a=[1, 3, 5]';
b=[2, 4, 6]';
A=[ [1, 0, 2]', [0, 1, 0]', [2, 0, 1]'];
B=[ [1, 0, -1]', [0, 2, 0]', [-1, 0, 3]'];
% dot product
s = a'*b;
s
% outer product
T = a*b';
T
% matrix times vector
c = A*a;
c
% matrix times matrix
P = A*B;
P
% element-wise multiplication of two vectors
d = a.*b;
d
% edama_01_11: accessing elements of a matrix
% define a column-vector a and a matrix B
a = [1, 2, 3]';
a
B = [ [1, 4, 7]', [2, 5, 8]', [3, 6, 9]'];
B
% 2nd element of vector a
s = a(2);
s
% row=2 column=3 element of matrix B
t = B(2,3);
t
% 2nd column of matrix B
b = B(:,2);
b
% 2nd row of matrix B, converted to column vector
c = B(2,:)';
c
% 2x2 sub-matrix in the bottom right of M
T = B(2:3,2:3);
T
% examples of using the colon notation, alone
x = [1:4];
x
y = [1:2:9];
y
z = [4:-1:1];
z
% edama_01_12: derivatives and integrals
clear all;
% independent variable t
N = 21; % length of t
tmin = 0; % minimum time
tmax = 1; % maximum time
Dt = (tmax-tmin)/(N-1); % time increment
t = tmin + Dt*[0:N-1]'; % time vector
% exemplary funtion
x = sin(pi*t);
% finite difference approximation of derivative
% note that its length is N-1
s = (x(2:N)-x(1:N-1))/Dt;
s'
% Riemann sum appoximation for integral
a = Dt*cumsum(x);
a'
% Riemann sum approximation for last value only
atotal = Dt*sum(x);
atotal
% the rest of this section makes plots used in the text,
% but since plotting has not yet been discussed, readers
% may prefer to skip it
% true derivative and integral
strue = pi*cos(pi*t); % derivative
atrue = 1/pi-cos(pi*t)/pi; % integral, with integration constant
% chosen so integral is zero at t=0
% time series plot
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis( [0, 1, 0, 1.1 ] );
xlabel('t');
ylabel('x(t)');
% plot function as continuous curve
plot( t, x, 'g-', 'LineWidth', 4, 'Color', [0.8,0.8,0.8] );
% plot function as circles
plot( t, x, 'ko', 'LineWidth', 2 );
% plot vertical tics
for i=[1:N-1]
plot([t(i),t(i)], [0,0.1], 'k-', 'LineWidth', 2 );
end
% calculus plot
figure(2);
clf;
subplot(3,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, 1, 0, 1.1 ] );
xlabel('t');
ylabel('x(t)');
% plot exemplary rectangles, shaded
K=5;
for k=[1:5]
if(x(k)>0)
r=rectangle( 'Position', [t(k),0,Dt,x(k)], 'FaceColor', [0.9,0.9,0.9] );
end
end
% plot function as continuous curve
plot( t, x, 'g-', 'LineWidth', 4, 'Color', [0.9,0.9,0.9] );
% plot function as circles
plot( t, x, 'ko', 'LineWidth', 2 );
% embolden one slope
K=5;
plot( [t(K),t(K+1)], [x(K),x(K+1)], 'k-', 'LineWidth', 6, 'Color', [0,0,0]);
% plot polygons
for i=[1:N-1]
plot([t(i),t(i),t(i+1),t(i+1)], [0,x(i),x(i+1),0], 'k-', 'LineWidth', 2 );
end
% plot derivative
subplot(3,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, 1, -1.1*pi, 1.1*pi ] );
xlabel('t');
ylabel('s=dx/dx');
plot( t, strue, 'k-', 'LineWidth', 4, 'Color', [0.8,0.8,0.8] );
plot( t(1:N-1), s, 'k-', 'LineWidth', 2 );
% plot integral
subplot(3,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, 1, 0, 1.1*2/pi ] );
xlabel('t');
ylabel('integral of x');
plot( t, atrue, 'k-', 'LineWidth', 4, 'Color', [0.8,0.8,0.8] );
plot( t, a, 'k-', 'LineWidth', 2 );
plot( t(N), atotal, 'ko', 'LineWidth', 3 ); % last value
% edama_01_13: for-loops
clear all;
% define a matrix A and a vector a
A = [ [1, 4, 7]', [2, 5, 8]', [3, 6, 9]'];
a = [0, 0, 0]';
% for loop that extracts the diagonal elements of a matrix
for i = [1:3]
a(i) = A(i,i);
end
A
a
% two nested for loops used to rearrange (flip left-right)
% the elements of a matrix
B = zeros(3,3);
for i = [1:3] % loop over all rows
for j = [1:3] % loop over all cols
B(i,j) = A(i,4-j); % left-right-flip
end
end
% print
A
B
% column vector of length 12
a = [ 1, 2, 1, 4, 3, 2, 6, 4, 9, 2, 1, 4 ]';
N = length(a);
% copy the vector a to the vector b, except that when
% an element of vector a exceeds 6, set the corresponding
% element of b to 6
b = zeros(N,1);
for i = [1:N]
% if greater than 6
if ( a(i) >= 6 )
% set to 6
b(i) = 6;
else
% else set to a
b(i) = a(i);
end
end
% side-by-size printing of a and b
c=[a, b];
conv2
% edama_01_14: alternatives to for-loops
% define a matrix M
B = [ [1, 4, 7]', [2, 5, 8]', [3, 6, 9]'];
% extract diagonal
a = diag(B);
A = [ [1, 4, 7]', [2, 5, 8]', [3, 6, 9]'];
% reverse order of rows
B = fliplr(A);
A
B
% another way of reversing order order of rows
B = A(:,end:-1:1);
A
B
% define vector a
a = [ 1, 2, 1, 4, 3, 2, 6, 4, 9, 2, 1, 4 ]';
% clip to 6 using logical arithmetic
b = a.*(a<6)+6.*(a>=6);
a'
b'
% alternative way to clip to 6, using find
b=a;
k = find(a>6);
b(k)=6;
a'
b'
% yet another alternative way to clip to 6, using logical addressing
b=a;
b(a>6)=6;
a'
b'
% edama_01_15: matrix inverse
% define a matrix A
A = [ [1, 5, 13]', [2, 7, 17]', [3, 11, 19]'];
% compute inverse
B = inv(A);
% print
A
B
% clek that B really is the inverse of A, by checking
% that multiplication yields the identity matrix
I1=A*B;
I2=B*A;
% print
I1
I2
% edama_01_16: solving linear system
% example of the matrix inverse via slash operator
% define matrices and vector
A = [ [1, 5, 13]', [2, 7, 17]', [3, 11, 19]'];
B = [ [1, 2, 3]', [4, 5, 6]', [7, 8, 9]'];
b = [2, 4, 6]';
% solve b=Ac for c via c=Ainv b
c = A\b;
% print
A
b
c
% check result is correct
error1=b-A*c;
error1
% solve D A = B for D via D = B Ainv
D = B/A;
% print
A
B
D
error2 = B-D*A;
error2
% edama_01_16
% loading the neuse.txt dataset
D=load('../Data/neuse.txt');
[N,K] = size(D); % determine length of data
t=D(:,1);
d=D(:,2);
N
K
% edama_01_17: load text file and make very-quick plot
% this example uses the neuse.txt stream flow dataset
D=load('../Data/neuse.txt');
t=D(:,1); % time in column 1
d=D(:,2); % discharge data in column 2
% quick plot, to visually inspect that data loosk reasonable
figure(1);
clf;
plot(t,d,'k-');
xlabel('time (days)');
ylabel('discharge (cfs)');
% edama_01_18: load text file and mmke nicer plot
% this example uses the neuse.txt stream flow dataset
D=load('../Data/neuse.txt');
t=D(:,1);
d=D(:,2);
figure(1);
clf;
hold on;
set(gca,'LineWidth',2);
set(gca,'FontSize',16);
axis( [0, 4500, 0, 1.1*max(d)] );
plot(t,d,'k-','LineWidth',2);
title('Neuse River Hydrograph');
xlabel('time (days)');
ylabel('discharge (cfs)');
% edama_01_19: log-log plot
% load the earthquake_magnitudes.txt data set
% convert from magniutude to enegy
% and make a log-log plot
clear all;
% get 10-year list of earthquake magnitudes from file
% note: since the file contains only one column of data, the
% np.genfromtxt puts the data into an N by 0 array. So we
% call eda_cvec to put it into a standard columnn vector
d = load('../Data/earthquake_magnitudes.txt');
time = 10; % 10 years of data
% use the histogram function to count up the
% numberof earthquake in a serie of magnitude bins
m = [5:0.1:10]'; % bins from 5 to 10 in increments of 0.1
count = hist(d,m)'; % count up earthquakes in each bin
r = count/time; % rate r is number of earthquakes divided by time
% Gutenberg-Richter magnitude-energy relation
log10Ee = 1.5*m + 11.8; % convert magnide to log-ergs
log10Ej = log10Ee - 7; % convert log-ergs to log-joules
E = 1e-7 * 10.^log10Ej; % convert log-joules to joules
% quick and dirty log log plot
figure(1);
clf;
hold on;
set(gca,'Xscale','log','Yscale','log');
plot( E, r, 'k-', 'LineWidth', 2);
plot( E, r, 'ko', 'LineWidth', 2);
xlabel('energy (joules)');
ylabel('earthquakes rate (per year)');
% nice looking log-log plot
figure(2);
clf;
hold on;
set(gca,'LineWidth',2);
set(gca,'FontSize',12);
set(gca,'Xscale','log','Yscale','log');
set(gca,'TickLength', [0.04 0.035]);
hold on;
axis([1e5,1e10,0.1,1000]);
plot( E, r, 'ko', 'LineWidth', 2 );
plot( E, r, 'k-', 'LineWidth', 2 );
xlabel('energy (joules)');
ylabel('earthquakes rate (per year)');
title('Earthquake rate vs energy');
% edama_01_20: write matrix to file
% load the neuse.txt data set
% convert discharge from cfs to m3/s
% and write to a file neuse_metric.txt
% get matrix from file
D=load('../Data/neuse.txt');
t=D(:,1); % time
d=D(:,2); % discharge
% conversion from cfs to m3/s
f=35.3146;
% convert data
dm = d/f;
% make simple plot to veryify result "by eye"
figure(1);
clf;
plot(t,dm,'k');
title('Neuse River Hydrograph');
xlabel('time in days');
ylabel('discharge in m3/s');
% concatenate time, discharge into a matrix C
C = [t, dm];
% save matrix to file of tab-separated values
dlmwrite('../Data/neuse_metric.txt',Dm,'\t');
% edama_01_21: big-picture commenting (recommended)
% This script evaluates the 2D Normal p.d.f. p(d1,d2) with
% mean (d1bar,d2bar) and covariance C on a L-by-L grid
% define da and d2 axes
L=101;
d1=[0:L-1]'/(L-1);
d2=[0:L-1]'/(L-1);
% mean
d1bar=0.5;
d2bar=0.5;
% covariance
C = [0.25, 0.1; 0.1, 0.5];
% evaluate the 2D normal probability density
% function p(d1,d2)
CI=inv(C);
N=1/(2*pi*sqrt(det(C))); % normalization
Pp=zeros(L,L);
for i = [1:L] % d2
for j = [1:L] % d1
dd = [d1(i)-d1bar,d2(j)-d2bar]';
% p(d1,d2)
p(i,j)=N*exp(-0.5*dd'*CI*dd);
end
end
figure(1);
clf;
imagesc(p);
title("normal p.d.f.")
% edama_01_22: minutia commenting (not recommended)
% create a matrix version of p(d1,d2)
% define length L vectors d1 and d2 on interbval (0,1)
L=101;
d1=[0:L-1]'/(L-1);
d2=[0:L-1]'/(L-1);
% two constants
d1bar=0.5;
d2bar=0.5;
% C is a 2x2 matrix
C = [0.25, 0.1; 0.1, 0.5];
% evaluate p(i,j)
CI=inv(C); % take inverse
N=1/(2*pi*sqrt(det(C))); % compute N
p=zeros(L,L); % set p to zero
for i = [1:L] % loop over i
for j = [1:L] % loop over j
% dd is a 2-vector
dd = [d1(i)-d1bar,d2(j)-d2bar]';
% cpmute d(i,j)
p(i,j)=N*exp(-0.5*dd'*CI*dd);
end
end
% plot p(i,j)
figure(1);
clf;
imagesc(p);
title("normal p.d.f.");