Live Script edama_02
This Live Script supports Chapter 2 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_02_01, print numeric values of a few data
clear all; % clears previously-defined variables
% performing a clear at the top of a LiveSript is a good idea
D=load('../Data/brf_temp.txt'); % read file into matrix D
t=D(:,1); % time vector t from column 1
d=D(:,2); % data (temperature) vector d from column 2
[N,M]=size(D); % D has N rows and M cols
N
M
% print the first few rows
D(1:5,1:2)
% edama_02_02: initial plot of the data
% Black Rock Forest temperature data used as an example
% read the data
D=load('../Data/brf_temp.txt');
t=D(:,1);00
d=D(:,2);
Ns=size(D);
N=Ns(1);
M=Ns(2);
% plot the data
figure(1); % open figure 1
clf; % clear figure
set(gca,'LineWidth',2); % thick lines
set(gca,'FontSize',14); % big font
hold on; % save changes
plot(t,d,'k-','LineWidth',2); % plot data
title('Black Rock Forest Weather Station'); % set title
xlabel('time (days)'); % set horizontal axis label
ylabel('temperature (deg C)'); % set vertical axis label
% edama_02_03: enlagred plot of part of the data
% Black Rock Forest temperature data used as an example
% make a plot of a section of data centered at a a time tc
% read the data
D=load('../Data/brf_temp.txt');
t=D(:,1);
d=D(:,2);
Ns=size(D);
N=Ns(1);
M=Ns(2);
tc = 1560.0; % center tine in seconds
Dt = t(2)-t(1); % sampling interval
wh = 25.0; % half width of window in days
iwh = floor(wh/Dt); % half width of window in samples
i=find((t>=tc),1); % find first index of t greater than tc
figure(1); % close-up plot, centered at tc, 2*w2 wide
clf;
set(gca,'LineWidth',2);
hold on;
% omit axis command, so plot auto-scales
plot( t(i-iwh:i+iwh), d(i-iwh:i+iwh), 'k-','LineWidth',2 );
plot( t(i-iwh:i+iwh), d(i-iwh:i+iwh), 'k*','LineWidth',2 );
title('Close-up: Black Rock Forest Temp');
xlabel('time in days after Jan 1, 1997');
ylabel('temperature');
% full plot in light gray with close-up highlighted in black
figure(2);
clf;
set(gca,'Linewidth',2);
hold on;
plot(t,d,'-','Linewidth',2,'Color',[0.8,0.8,0.8] );
plot( t(i-iwh:i+iwh), d(i-iwh:i+iwh), 'k-','LineWidth',3 );
title('Black Rock Forest Weather Station');
xlabel('time, days');
ylabel('temperature, C');
% edama_02_04: enlarged plot via mouse-click
% Black Rock Forest temperature data used as an example
% make an enlarged plot of a section of data
% centered at a time tc
% and with a half with of w2 samples
% read the data
D=load('../Data/brf_temp.txt');
t=D(:,1);
d=D(:,2);
Ns=size(D);
N=Ns(1);
M=Ns(2);
figure(1); % the interactive plot
clf;
set(gca,'Linewidth',2);
hold on;
plot(t,d,'k-','Linewidth',2);
title('Black Rock Forest Weather Station');
xlabel('time, days');
ylabel('temperature, C');
% plot section based on cursor click
[tc, dc] = ginput(1); % detect mouse click
i=find((t>=tc),1); % find first index of t greater than tc
close(1); % close the interative plot
figure(2); % close-up plot
clf;
set(gca,'LineWidth',2);
hold on;
Dt = t(2)-t(1); % sampling interval
wh = 25.0; % half width of window in days
iwh = floor(wh/Dt); % half width of window in samples
plot( t(i-iwh:i+iwh), d(i-iwh:i+iwh), 'k-','LineWidth',2 );
plot( t(i-iwh:i+iwh), d(i-iwh:i+iwh), 'k*','LineWidth',2 );
title('Black Rock Forest Temp');
xlabel('time in days after Jan 1, 1997');
ylabel('temperature');
% full plot in light gray with close-up highlighted in black
figure(3);
clf;
set(gca,'Linewidth',2);
hold on;
plot(t,d,'-','Linewidth',2,'Color',[0.8,0.8,0.8] );
plot( t(i-iwh:i+iwh), d(i-iwh:i+iwh), 'k-','LineWidth',2 );
plot( t(i-iwh:i+iwh), d(i-iwh:i+iwh), 'k.','LineWidth',2 );
title('Black Rock Forest Weather Station');
xlabel('time, days');
ylabel('temperature, C');
% edama_02_05: histogram of data
% Black Rock Forest temperature data used as an example
% read the data
D=load('../Data/brf_temp.txt');
t=D(:,1);
d=D(:,2);
Ns=size(D);
N=Ns(1);
M=Ns(2);
% compute histogram
Lh = 100; % number of bins in histogram
dmin = min(d); % lowest temperature bin
dmax = max(d); % highest temperature bin
centers = dmin + (dmax-dmin)*[1:Lh]'/Lh; % bin centers
counts = hist(d, centers)'; % counts in each bin
% plot histogram
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
plot(centers, counts, 'k', 'LineWidth',2);
title('Histogram of temperature');
xlabel('tempeature bin');
ylabel('counts');
% edama_02_06: plot histogram two ways
% Black Rock Forest temperature data used as an example
% computes and plot histogram two ways, as
% a graph and as a grey scale column-vector
% read the data
D=load('../Data/brf_temp.txt');
t=D(:,1);
d=D(:,2);
Ns=size(D);
N=Ns(1);
M=Ns(2);
% compute histogram
Lh = 100; % number of bins in histogram
dmin = min(d); % lowest temperature bin
dmax = max(d); % highest temperature bin
centers = dmin + (dmax-dmin)*[1:Lh]'/Lh; % bin centers
counts = hist(d, centers)'; % counts in each bin
figure(1);
clf;
subplot(1,2,1);
set(gca,'LineWidth',2);
hold on;
axis ij;
axis( [0, max(counts), -60, 40]);
plot(counts, centers, 'k','LineWidth',2);
title('Histogram of temperature');
xlabel('tempeature bin');
ylabel('counts');
subplot(1,2,2);
set(gca,'LineWidth',2);
axis([0, 1, -60, 40]);
hold on;
axis ij;
% axis off;
set(gca,'xtick',[]);
% create grey-scale color map
bw=0.9*(256-[0:255]')/256;
colormap([bw,bw,bw]);
% plot the vector dhist
imagesc( [0.4, 0.6], [centers(1), centers(end)], counts);
cb=colorbar('vert');
set(get(cb,'title'),'string','counts');
% edama_02_07, moving-window histogram of data
% Black Rock Forest temperature data used as an example
% read the data
D=load('../Data/brf_temp.txt');
t=D(:,1);
d=D(:,2);
Ns=size(D);
N=Ns(1);
M=Ns(2);
Dt = 1.0/24.0;
% compute histogram
Lh = 100; % number of bins in histogram
dmin = -60.0; % lowest temperature bin
dmax = 40.0; % highest temperature bin
centers = dmin + (dmax-dmin)*[1:Lh]'/Lh; % bin centers
% Lw moving window histograms, each toffset days long
toffset = 30.0; % offset of 30 days
offset = floor(toffset/Dt); % offset in samples
Lw=floor(N/offset)-1; % number of windows
tend = toffset*(Lw-1); % last time
Dw = zeros(Lh, Lw); % matrix of histograms
for i = [1:Lw] % loop over windows
j=1+(i-1)*offset; % index of first data in window
k=j+offset-1; % index of last data
Dw(:,i) = hist(d(j:k), centers)'; % histogram
end
% Set up graphics
figure(1);
clf;
set(gca,'LineWidth',2);
axis([0, tend, -60, 40]);
set(gca, 'XAxisLocation', 'top')
hold on;
axis ij;
% grey-scale color map
bw=0.9*(256-[0:255]')/256;
colormap([bw,bw,bw]);
% plot the matrix Dhist
imagesc( [0, tend], [dmin, dmax], Dw);
xlabel('time (days)');
ylabel('temperature (deg C)');
title('moving window histogram');
cb=colorbar('vert');
set(get(cb,'title'),'string','counts');
% edama_02_08: superimposed plots with different line types
% create sample data
N=51;
Dt = 1.0;
t = [0:N-1];
tmax=t(N);
d1 = sin(pi*t/tmax);
d2 = sin(2*pi*t/tmax);
% plot the sample data
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis xy;
axis([0, tmax, -1.1, 1.1]);
plot(t,d1,'k-','LineWidth',2);
plot(t,d2,'k:','LineWidth',2);
title('data consisting of sine waves');
xlabel('time');
ylabel('data');
% edama_02_09, sideways side-by-side plots
% create sample data
N=51;
Dt = 1.0;
t = [0:N-1];
tmax=t(N);
d1 = sin(pi*t/tmax);
d2 = sin(2*pi*t/tmax);
% plot the sample data
figure(1);
clf;
subplot(1,2,1);
set(gca, 'XAxisLocation','Top');
set(gca, 'YAxisLocation','Left');
set(gca,'LineWidth',2);
hold on;
axis([-1.1, 1.1, 0, tmax]);
axis ij;
plot(d1,t,'k-','LineWidth',2);
plot([0,0],[0,tmax],'k:','LineWidth',2);
title('d1');
ylabel('time');
xlabel('data');
subplot(1,2,2);
set(gca, 'XAxisLocation','Top');
set(gca, 'YAxisLocation','Left');
set(gca,'LineWidth',2);
hold on;
axis ij;
axis([-1.1, 1.1, 0, tmax]);
plot(d2,t,'k-','LineWidth',2);
plot([0,0],[0,tmax],'k:','LineWidth',2);
title('d2');
ylabel('time');
xlabel('data');
% eda02_10, colormaps
% random image
N = 11;
D = random('Uniform',0.0,1.0,N,N);
% MATLAB-provided grey-scale colormap
gray(256);
% plot
figure(1);
clf;
imagesc(D);
title('MATLAB''s grey colormap');
% alternaitive is to build a user-defined colormap
% note that this colormap is reversed; white is low
bw=0.9*(256-[0:255]')/256;
colormap([bw,bw,bw]);
figure(1);
clf;
imagesc(D);
title('image plotted with MATLAB''s grey colormap');
% edama_02_11: example of plotting an image
% based on edama_02_07, running histogram of temperature
% read the data
D=load('../Data/brf_temp.txt');
t=D(:,1);
d=D(:,2);
Ns=size(D);
N=Ns(1);
M=Ns(2);
Dt = 1.0/24.0;
% compute histogram
Lh = 100; % number of bins in histogram
dmin = -60.0; % lowest temperature bin
dmax = 40.0; % highest temperature bin
centers = dmin + (dmax-dmin)*[1:Lh]'/Lh; % bin centers
% Lw moving window histograms, each toffset days long
toffset = 30.0; % offset of 30 days
offset = floor(toffset/Dt); % offset in samples
Lw=floor(N/offset)-1; % number of windows
tend = toffset*(Lw-1); % last time
Dw = zeros(Lh, Lw); % matrix of histograms
for i = [1:Lw] % loop over windows
j=1+(i-1)*offset; % index of first data in window
k=j+offset-1; % index of last data
Dw(:,i) = hist(d(j:k), centers)'; % histogram
end
% grey-scale color map
bw=0.9*(256-[0:255]')/256;
colormap([bw,bw,bw]);
% Set up graphics
figure(1); % opem figure
clf; % clear previous contents
set(gca,'LineWidth',2); % thicker line in axe
axis([0, tend, dmin, dmax]); % plot axes
set(gca, 'XAxisLocation', 'top'); % axis on top of plot
hold on; % don't override previous setting
axis ij; % 'matrix-sytpe orientaion
vlo = 0.0;
vhi = max(max(Dw));
imagesc( [0, tend], [dmin, dmax], Dw, [vlo, vhi]); % plot image
xlabel('time (days)'); % horizontal axis label
ylabel('temperature (deg C)'); % vertical axis label
title('moving window histogram'); % title
cb=colorbar('vert'); % plot colorbar
set(get(cb,'title'),'string','counts'); % set colorbar title
% edama_02_12: plot and histgram of rates of change
% using the Neuse River Hydrograph data as an example
% load data
D=load('../Data/neuse.txt');
t=D(:,1);
d=D(:,2);
N=length(d);
% plot only these bounds (in samples)
N1=500;
N2=1000;
figure(1);
clf;
% plot discharge
subplot(1,3,1);
set(gca,'LineWidth',2);
axis ij;
hold on;
plot(d(N1:N2),t(N1:N2),'k-','LineWidth',2);
title('A)');
xlabel('discharge');
ylabel('time');
% compute rate
dddt = (d(2:N)-d(1:N-1)) ./ (t(2:N)-t(1:N-1));
% plot rate
subplot(1,3,2);
set(gca,'LineWidth',2);
axis ij;
hold on;
plot(dddt(N1:N2),t(N1:N2),'k-','LineWidth',2);
title('B');
xlabel('d/dt discharge');
ylabel('time');
% compute histogram
Nh=100;
dddtmin=-500;
dddtmax=500;
bins = dddtmin + (dddtmax-dddtmin)*[0:Nh-1]'/(Nh-1);
dddth = hist(dddt,bins);
% delete first and last bins
dddth = dddth(2:Nh-1);
bins = bins(2:Nh-1);
% plot histogram
subplot(1,3,3);
set(gca,'LineWidth',2);
axis ij;
hold on;
plot(dddth,bins,'k-','LineWidth',2);
% plot dotted line at counts=0
a = axis;
plot( [0,a(2)]', [0,0]', 'k:','LineWidth',2);
title('C)');
xlabel('counts');
ylabel('d/dt discharge');
% edama_02_13: segregating data based on the sign of rate
% using the Neuse River Hydrograph data as an example
% load data
D=load('../Data/neuse.txt');
t=D(:,1);
d=D(:,2);
Ls=size(d);
L=Ls(1);
figure(1);
clf;
% compute rate
dddt=(d(2:L)-d(1:L-1)) ./ (t(2:L)-t(1:L-1));
% divide (d, dddt) pairs into two groups depending
% upon whether dddt is positive or negative
pos = find(dddt>0);
neg = find(dddt<0);
dpos = d(pos);
dneg = d(neg);
dddtpos = dddt(pos);
dddtneg = dddt(neg);
% scatter plot of data vs rates for positive rates
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis xy;
axis( [0 10000 0 5000] );
plot(dpos,dddtpos,'k.','LineWidth',2);
title('A) positive rates');
xlabel('discharge');
ylabel('d/dt discharge');
% scatter plot of data vs rates for negative rates
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis xy;
axis( [0 10000 -5000 0] );
plot(dneg,dddtneg,'k.','LineWidth',2);
title('B) negative rates');
xlabel('discharge');
ylabel('d/dt discharge');
% edama_02_14: scatter plots of Atlantic Rock chemical data
% load 'rocks.txt' and assign columns to mnemonic variable names
% file contaisn: SiO2 TiO2 Al2O3 FeOT MgO CaO Na2O K2O
D = load('../Data/rocks.txt');
sio2 = D(:,1); % SiO2
tio2 = D(:,2); % TiO2
als03 = D(:,3); % Al203
feot = D(:,4); % FeO-total
mgo = D(:,5); % MgO
cao = D(:,6); % CaO
na20 = D(:,7); % Na2O
k20 = D(:,8); % K2O
Ns = size(D);
N = Ns(1);
M = Ns(2);
% note the use of a cellstr to hold chemical names
names = { 'SiO2', 'TiO2', 'Al203', 'FeO-total', ...
'MgO', 'CaO', 'Na2O', 'K2O' };
[temp, Ns] = size(names); % number of chemical names
% Np: number of permulations of Ns things taken two at a time
Np = Ns*Ns-(Ns*(Ns+1)/2);
M = 4; % plots in a row
Nr = ceil(Np/M); % number of rows
% scatter plots
figure(1);
clf;
sp=1; % count plots
for i = [1:Ns-1]
for j = [i+1:Ns]
subplot(Nr,M,sp);
hold on;
set(gca,'LineWidth',2);
set(gca,'FontSize',12);
plot( D(:,i), D(:,j), 'k.', 'LineWidth',2 );
xlabel( char(names(i)) );
ylabel( char(names(j)) );
sp=sp+1; % increment coumter
end
end
% edama_02_15: character strings
% Part 1. String variables
my_name = 'Bill';
my_path = 'C:/bill';
disp(my_name);
disp(my_path);
disp(" "); % blank line
% Part 2. Using cellstrs to hold strings
names = { 'lithium', 'calcium', 'sodium', 'potassium' };
[temp, Ns] = size(names); % number of names
disp('number of names:');
disp(Ns);
disp('name 1:');
disp(char(names(1)));
disp('name 2:');
disp(char(names(2)));
disp('name 3:');
disp(char(names(3)));
disp('name 4:');
disp(char(names(4)));
disp(" "); % blank line
% Part 3. Formatting strings using the % operator
% the next two examples use sprintf() to create a string
% and then print the string using disp()
% prints "element 2"
i=2;
mystring1 = sprintf('element %d', i);
disp(mystring1);
% prints "a=4.72"
a=4.7213;
mystring2 = sprintf('a=%.2f', a);
disp(mystring2);
% the next two examples use fprintf() to create and print
% the string, without creating a string variable. Note that
% a newline character (the character '\n') must be added
% to the end of the format to advance to the next line
% prints "row 2 column 4"
i=2;
j=4;
fprintf('row %d column %d\n', i, j);
% prints "(2,4): 4.721"
i=2;
j=4;
a=4.7213;
fprintf('(%d,%d): %.3f\n', i, j, a);
disp(" "); % blank line
% Part 4. Combining lists and formats
for i = [1:Ns]
fprintf('element %d is %s\n', i, char(names(i)) );
end
disp(" ");