ESS 116 Introduction to Data Analysis in Earth Science
ESS 116
Introduction to Data Analysis in Earth Science
Image Credit: NASA
Instructor: Mathieu Morlighem
E-mail: mmorligh@uci.edu (include ESS116 in subject line)
Office Hours: 3218 Croul Hall, Friday 2:00 pm – 3:00 pm
This content is protected and may not be shared uploaded or distributed
Lecture 8 quick review
Lecture 9 – Time Series & Bivariate statistics
What is a time series?
Trend
Cycles (diurnal, seasonal)
Regressions
Bivariate statistics: Correlation between two variables x and y
Final thoughts
Today’s lecture
Two types of images:
Vector graphics (made of vector objects)
Raster graphics (made of pixels, each pixel has a uniform color)
RGB color model:
Each color is represented by its level of Red, Green and Blue
uint8: from 0 to 255
double: from 0 to 1
Example (for the class double):
Red: [1 0 0], White: [1 1 1], Yellow: [1 1 0], Gray: [0.5 0.5 0.5]
Reading an image with MATLAB: A=imread(‘filepath’)
Display an image with MATLAB: image(A)
Review – Lecture 8
Image indexed to a colormap
Represented by a 2D matrix of numbers referring to a row in the colormap
Review – Lecture 8
Image Matrix
Colormap Matrix
4 Pixel Image
Colormap size should be consistent with the image matrix
If not, use imagesc(A) instead (scales the colorbar)
Review – Lecture 8
“True color” Image (RGB)
Represented by a 3D matrix (3 slices, for Red Green and Blue)
(:, :, 1) R-values; (:, :, 2) G-values; (:, :, 3) B-values
N columns
M rows
i>Clicker question
To create a Russian flag, what should we use?
iMat(:,:,1) = [1;0;1];
iMat(:,:,2) = [1;0;0];
iMat(:,:,3) = [1;1;0];
iMat = [1;2;3]; colormap([1 1 1;0 0 1;1 0 0])
Both A and B
None of the above
Red = [1 0 0]
Blue = [0 0 1]
White = [1 1 1]
Lecture 9 – Time Series
& Bivariate statistics
Time series
Time Series: A collection of observations xt, each one being recorded at time t
Time Series
Objective of Time Series Analysis:
Provide compact mathematical description of the data
Describe trends and cycles (seasonality)
Signal processing: extracting a signal in the presence of noise
Prediction: use model to predict future values of the time series
Ozone concentration
Keeling curve (CO2)
Arctic sea ice
Decomposing a time series:
= + +
t: time
m(t): trend component (not periodic, slowly changing in time)
s(t) : cycle component (systematic, calendar related movements)
e(t) : random error component (unsystematic, short term fluctuations)
Components of a time series
Hands on Example
%Load two variables in file: time and %seaice
load(‘SeaiceDailyExtent.mat’);
%Plot time series
plot(time,seaice,’-b’);
axis([1979 2020 4 17])
title(‘Artic sea ice surface area’);
xlabel(‘time (year)’);
ylabel(‘millions km^2′);
Trend m(t): long term direction
Option 1: filters such as doing a moving average can also be used to get a more smooth inter-annual trend:
MATLAB command: m = movmean(y,12)
Option 2: Polynomial fitting (generally up to a degree 2)
MATLAB command: p = polyfit(t,y,2);
m = polyval(p,t);
Estimation of the components
Hands on Example
%Compute trend m (long term direction)
%———————————–
%OPTION 1: use MATLAB’s movmean (window averaging)
%we have 365 points per year, so we need to average with a window of 1
%year)
y = seaice;
m = movmean(y,365);
%OPTION 2: use polynomial fit. Here, it looks like a linear decrease, so
%use polynomial of degree 1 only
p = polyfit(time,y,1);
m = polyval(p,time);
%Plot time series (blue) and its trend (red)
plot(time,seaice,’-b’);
hold on
plot(time,m,’-r’);
axis([1979 2020 4 17])
title(‘Artic sea ice surface area’);
xlabel(‘time (year)’);
ylabel(‘millions km^2’);
legend(‘Artic sea ice surface area’,’Trend’);
Seasonal Cycles (i.e. sines and cosines):
Cycles have an amplitude and frequency (1/period)
Sine/cosine naturally model cycles
T is the period (use the same unit as the unit of the time variable!)
Example: for a cycle of 1 year (seasons)
If time is in year, T = 1
If time is in month, T = 12
If time is in days, T = 365
Estimation of the components
The concentration of O3 in the atmosphere varies with sunlight (i.e. it has a daily cycle). What would be the value of the period T to best fit a cycle for which the time is expressed in days?
T = 1
T = 12
T = 24
T = 3600*24
i>Clicker question
Seasonal Cycles FYI only: Why 2 terms?
A seasonal cycle “should” look like
Where are 2 unknowns are:
amplitude (⍺)
phase shift (β)
Using the formula cos(a+b) = cos(a) cos(b) – sin(a) sin(b)
Estimation of the components
If y is a time series, and y’ the detrended time series:
y’(t) = y(t) – m(t)
Detrended time series
Original data y(t)
Detrended data y’(t)
(centered around y=0)
We want to find A and B such that
is “as close as possible” to the detrended series
To do the regression, we are going to use a “Least square fit”!
Regression
We express the “modeled” values as:
Regression
To find our unknown parameters (A and B in vector β):
Construct the design matrix X
We want to be as close as possible to y’.
Solve (you don’t need to know “the magic” behind least square fit)
Hands on Example
%Detrend data
y_prime = y – m;
%Construct design matrix
X = [cos(2*pi/1*time) sin(2*pi/1*time)];
%Find coefficients (A, B)
beta = inv(X’*X)*X’*y_prime;
%Seasonal cycle
A = beta(1); B = beta(2);
s = A*cos(2*pi/1*time) + B*sin(2*pi/1*time);
%Plot detrended data and seasonal cycle
plot(time,y_prime,’-k’);
hold on
plot(time,s,’–r’);
xlim([1979 2020])
xlabel(‘time (year)’);
ylabel(‘millions km^2’);
legend(‘Detrended data’,’Seasonal cycle’)
Noise (i.e. what’s left…):
We found the trend and seasonal cycle, the remaining term is just noise that we could not fit (due to natural variability)
Estimation of the components
%Find noise (what’s left)
e = y – m – s;
%Plot all components!
subplot(2,1,1)
plot(time,seaice,’-b’);
axis([1979 2020 4 17])
title(‘Artic sea ice surface area’);
xlabel(‘time (year)’);
subplot(2,3,4)
plot(time,m,’-g’);
title(‘Trend’); xlim([1979 2020]);
subplot(2,3,5)
plot(time,s,’-k’);
title(‘Seasonal cycle’); xlim([1979 2020]);
subplot(2,3,6)
plot(time,e,’-r’);
title(‘Noise’); xlim([1979 2020]);
Estimation via regression
We can actually use the same methodology to find both the trend and the seasonal cycle
Write the expression of a model that would best describe the time series
Trend can be constant, linear, quadratic
Seasonal cycle has a given period
Write the corresponding design matrix X
Solve the following equation
beta = inv(X’*X)*X’*y;
Using regression for both
the trend and seasonal cycle
Which of the following decompositions best matches the hypothetical time-series to the right?
Don’t know
i>Clicker Question
For a quadratic trend and a seasonal cycle of period 12 months:
Example of regression
Hands on Example
%Design matrix: y = A + B*t + C*sin(omega t) + D*cos(omega t)
n = length(y);
X = [ones(n,1) time cos(2*pi/1*time) sin(2*pi/1*time)];
%Find coefficients (A, B, C and D)
beta = inv(X’*X)*X’*y;
%Find trend, seasonal cycle and noise
m = beta(1)+beta(2)*time;
s = beta(3)*cos(2*pi/1*time) + beta(4)*sin(2*pi/1*time);
yhat = m + s;
e = y – yhat;
%plot data and model
subplot(2,1,1)
plot(time,y,’-b’);
hold on
plot(time,yhat,’–r’);
legend(‘Data’,’Model (best fit)’);
xlim([1979 2020]);
%plot individual components
subplot(2,3,4)
plot(time,m,’-g’); title(‘Trend’); xlim([1979 2020]);
subplot(2,3,5)
plot(time,s,’-k’); title(‘Seasonal cycle’); xlim([1979 2020]);
subplot(2,3,6)
plot(time,e,’-r’); title(‘Noise’); xlim([1979 2020]);
If we want to fit a time series with a seasonal cycle (12 months), and no trend (m=0), what would the design matrix look like?
C.
D. don’t know
i>Clicker Question
Bivariate statistics (Correlation)
Correlation Coefficient (r)
Correlation Coefficient: quantifies (or measures) the extent and nature of the relationship
Pearson’s correlation coefficient:
Exactly –1 indicates a perfect negative linear relationship.
Close to –1 indicates a strong negative linear relationship.
Close to 0 means no linear relationship exists.
Close to +1 indicates a strong positive linear relationship.
Exactly +1 indicates a perfect positive linear relationship.
No unit
Sometimes referred to as r2 or r
-1 ≤ r ≤1
Change the order of X and Y won’t change the value of r.
MATLAB: R=corrcoef(A,B) returns a 2×2 matrix
Make sure to use R(1,2) or R(2,1)
Correlation Coefficient (r)
Some scatterplots and correlations
Correlation assesses the LINEAR relationship!
Which pair of data in the graph do you think has highest correlation coefficient?
A.
B.
i>Clicker Question
C. Both of them are low
Correlation
Looking at the previous section on time series, we can quantify the “goodness of fit” between the model ) and the data (y)
>> corrcoef(y,yhat)
ans =
1.0000 0.9771
0.9771 1.0000
Or other quantities!
Ice cream consumption vs air temperature
Temperature vs concentration of CO2
Etc…
Example
Correlation vs Causation…
Correlation Doesn’t Necessarily Mean Cause-and-Effect
Cause and Effect
Correlation Doesn’t Necessarily Mean Cause-and-Effect
Cause-and-effect relationship is when a change in X causes a change in Y.
increase of CO2 global air temperature increase
global air temperature increase drop of ice sheet mass
Correlation between ice cream consumption vs murder rates: strong linear relationship.
more ice cream consumption causes more murders ???
One way to address this problem is to use a Lead-Lag analysis
Cause and Effect
Lead-Lag Correlation
Rainfall in March
Correlation Coefficient between soil moisture and March precipitation
http://www.mvd.usace.army.mil/About/MississippiRiverCommission(MRC).aspx
Cause and Effect
Conclusion
MATLAB:
Use command line
Create functions and scripts
Process big matrices
Read almost any input file (using load or textscan)
Plot curves, Histograms and images
Curve fitting, interpolation & Extrapolation
Image processing (True color & indexed to a color map)
Programming:
Selection statements (if/elseif/else)
For loops
Nested loops
Type casting
fprintf
What you have learned
Statistics
Descriptive statistics & Sample visualization
Properties of a distribution:
Central tendency
Dispersion
General shape
Theoretical distributions
Central Limit Theorem
Hypothesis testing (Student’s t-test)
Time series decomposition
Regression
Correlation
What you have learned
We have covered many different disciplines:
Statistics
Probability
Data visualization
Hypothesis testing
Interpolation, Extrapolation, Prediction
Computer Science
You have acquired a new skill:
Programming !
What you have learned
MATLAB is a powerful tool for processing quantitative data
MATLAB is not the only tool for data analysis
Is not ideal for all analyses, but is very good for most
Not all Earth scientists know how to code… but they should!
Knowing how to write your own code gives you the freedom to create customized tools
Tools that save time
Reduce errors from repetitive tasks
Avoid re-doing data analyses multiple times
Coding allows you to develop new research methods
Don’t need to wait for a program to have a button. You can be the one that makes the buttons
Final Thoughts
Final Thoughts
When the quarter’s dust settles….
Update your resume! You’ve earned it!
Definitely use this in your job & grad school applications
Thanks to our TA (Chia-Chun) for all her hard work this quarter
Thanks to all of you for your hard work, I hope you leave with a better understanding of how we analyze data, process images and test a hypothesis. I hope you all enjoyed the power of programming !
Lab 9: Time series
Lecture 10: review + sample exam
Come with questions!
What’s next?
/private/var/folders/cy/_n4ycj5x7f9699n7qpqg5c200000gn/T/tp041e57a9_7887_434d_90df_2d587e41071f.eps
1980 1985 1990 1995 2000 2005 2010 2015 2020
time (year)
4
6
8
10
12
14
16
m
il
li
o
n
s
k
m
2
2000200220042006200820102012201420162018
year
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
o
z
o
n
e
m
i
x
i
n
g
r
a
t
i
o
[
p
p
b
]
1234567
Day of Week [week of July 1 2016]
-0.04
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
0.04
o
z
o
n
e
m
i
x
i
n
g
r
a
t
i
o
[
p
p
b
]
1234567
Day of Week [week of July 1 2016]
-0.01
-0.005
0
0.005
0.01
0.015
0.02
0.025
0.03
o
z
o
n
e
m
i
x
i
n
g
r
a
t
i
o
[
p
p
b
]
1234567
Day of Week [week of July 1 2016]
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
o
z
o
n
e
m
i
x
i
n
g
r
a
t
i
o
[
p
p
b
]
y (t) = m (t)
| {z }
Trend
+ s (t)
| {z }
seasonal cycle
+ e (t)
| {z }
noise
1234567
Day of Week [week of July 1 2016]
0.02
0.025
0.03
0.035
0.04
0.045
0.05
0.055
o
z
o
n
e
m
i
x
i
n
g
r
a
t
i
o
[
p
p
b
]
/private/var/folders/cy/_n4ycj5x7f9699n7qpqg5c200000gn/T/tp6f15cccb_8d68_4542_b665_c3df46ff923e.eps
1980 1985 1990 1995 2000 2005 2010 2015 2020
time (year)
4
6
8
10
12
14
16
m
il
li
o
n
s
k
m
2
Artic sea ice surface area
m (t) =
0.5⇥ xt�6 + xt�5 + …+ xt+5 + 0.5⇥ xt+6
12
m (t) = a0 + a1t+ …+ akt
k
/private/var/folders/cy/_n4ycj5x7f9699n7qpqg5c200000gn/T/tpe4d4c945_6329_4f5f_86f9_4e192f3a6d1d.eps
1980 1985 1990 1995 2000 2005 2010 2015 2020
time (year)
4
6
8
10
12
14
16
m
il
li
o
n
s
k
m
2
Artic sea ice surface area
/private/var/folders/cy/_n4ycj5x7f9699n7qpqg5c200000gn/T/tp5efb61ca_9a32_4302_be92_6958d4db53ba.eps
1980 1985 1990 1995 2000 2005 2010 2015 2020
time (year)
4
6
8
10
12
14
16
m
il
li
o
n
s
k
m
2
Artic sea ice surface area
Artic sea ice surface area
Trend
s (t) = A cos
✓
2⇡
T
t
◆
+B sin
✓
2⇡
T
t
◆
s (t) = ↵ cos (�)
| {z }
A
cos
✓
2⇡
T
t
◆
+ (�↵ sin (�))
| {z }
B
sin
✓
2⇡
T
t
◆
s (t) = ↵ cos
✓
2⇡
T
t+ �
◆
/private/var/folders/cy/_n4ycj5x7f9699n7qpqg5c200000gn/T/tp9826bfe7_09c8_412e_8690_710541fc6cea.eps
1980 1985 1990 1995 2000 2005 2010 2015 2020
time (year)
-8
-6
-4
-2
0
2
4
6
m
il
li
o
n
s
k
m
2
y0 (t) = y (t)�m (t) = s (t)
| {z }
seasonal cycle
+ e (t)
| {z }
noise
� =
�
XTX
��1
XT y0
2
6666666666666
4
ŷ0 (t1)
ŷ0 (t2)
ŷ0 (t3)
…
ŷ0 (tn)
3
7777777777777
5
=
2
6666666666666
4
A cos
✓
2⇡
T
t1
◆
+B sin
✓
2⇡
T
t1
◆
A cos
✓
2⇡
T
t2
◆
+B sin
✓
2⇡
T
t2
◆
A cos
✓
2⇡
T
t3
◆
+B sin
✓
2⇡
T
t3
◆
…
A cos
✓
2⇡
T
tn
◆
+B sin
✓
2⇡
T
tn
◆
3
7777777777777
5
=
2
6666666666666
4
cos
✓
2⇡
T
t1
◆
sin
✓
2⇡
T
t1
◆
cos
✓
2⇡
T
t2
◆
sin
✓
2⇡
T
t2
◆
cos
✓
2⇡
T
t3
◆
sin
✓
2⇡
T
t3
◆
…
cos
✓
2⇡
T
tn
◆
sin
✓
2⇡
T
tn
◆
3
7777777777777
5
| {z }
Design matrix, X
2
4
A
B
3
5
| {z }
Unkowns, �
= X�
/private/var/folders/cy/_n4ycj5x7f9699n7qpqg5c200000gn/T/tp98f47f36_dbc0_4cfb_9d51_7b08a913229a.eps
1980 1985 1990 1995 2000 2005 2010 2015 2020
time (year)
-8
-6
-4
-2
0
2
4
6
m
il
li
o
n
s
k
m
2
Detrended data
Seasonnal cycle
/private/var/folders/cy/_n4ycj5x7f9699n7qpqg5c200000gn/T/tp2a9d59e5_f23c_484a_a139_20dbc77123f4.eps
1980 1985 1990 1995 2000 2005 2010 2015 2020
time (year)
4
6
8
10
12
14
16
m
il
li
o
n
s
k
m
2
Artic sea ice surface area
1980 1990 2000 2010 2020
10
10.5
11
11.5
12
12.5
13
Trend
1980 1990 2000 2010 2020
-5
-4
-3
-2
-1
0
1
2
3
4
5
Seasonal cycle
1980 1990 2000 2010 2020
-3
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Noise
ŷ (t) = A cos
✓
2⇡
12
t
◆
+B sin
✓
2⇡
12
t
◆
ŷ (t) = A+B t+ C cos
✓
2⇡
12
t
◆
+D sin
✓
2⇡
12
t
◆
ŷ (t) = A+B t+ C t2 +D cos
✓
2⇡
12
t
◆
+ E sin
✓
2⇡
12
t
◆
2
666666666
4
ŷ (t1)
ŷ (t2)
…
ŷ (tn)
3
777777777
5
=
2
666666666
4
A+Bt1 + Ct
2
1 +D cos
✓
2⇡
T
t1
◆
+ E sin
✓
2⇡
T
t1
◆
A+Bt2 + Ct
2
2 +D cos
✓
2⇡
T
t2
◆
+ E sin
✓
2⇡
T
t2
◆
…
A+Btn + Ct
2
n +D cos
✓
2⇡
T
tn
◆
+ E sin
✓
2⇡
T
tn
◆
3
777777777
5
=
2
666666666
4
1 t1 t
2
1 cos
✓
2⇡
T
t1
◆
sin
✓
2⇡
T
t1
◆
1 t2 t
2
2 cos
✓
2⇡
T
t2
◆
sin
✓
2⇡
T
t2
◆
…
…
…
…
…
1 tn t
2
n cos
✓
2⇡
T
tn
◆
sin
✓
2⇡
T
tn
◆
3
777777777
5
| {z }
Design matrix, X
2
666666666666
4
A
B
C
D
E
3
777777777777
5
| {z }
Unkowns, �
= X�
< l a t e x i t s h a 1 _ b a s e 6 4 = " Q Z h d 3 b b Q K V d W l Z X O U j D q H N 4 2 6 8 8 = " >
A
A
A
G
J
3
i
c
j
V
R
d
b
9
M
w
F
M
1
G
C
y
N
8
d
f
D
I
i
0
V
H
N
b
S
p
a
q
J
K
8
D
J
p
t
J
v
E
4
5
D
W
r
V
J
d
K
s
d
x
W
q
u
J
E
9
n
u
W
B
X
5
3
/
D
C
X
+
E
F
C
R
C
C
R
/
4
J
T
l
y
V
d
i
1
r
r
M
T
3
6
t
7
r
e
0
7
O
b
e
0
l
I
R
W
y
0
f
i
9
t
X
2
n
V
L
5
7
b
+
e
+
/
e
D
h
o
8
d
P
K
r
t
P
L
0
Q
8
4
Z
h
0
c
B
z
G
v
O
s
h
Q
U
L
K
S
E
d
S
G
Z
J
u
w
g
m
K
v
J
B
c
e
u
N
2
l
r
+
8
I
l
z
Q
m
J
3
L
a
U
L
6
E
R
o
y
G
l
C
M
p
A
4
N
d
k
t
H
0
C
N
D
y
l
L
E
O
Z
q
q
l
O
N
Q
2
T
A
k
g
e
w
t
J
7
A
O
j
5
B
M
p
y
r
P
7
s
u
B
A
z
k
d
j
u
Q
r
C
O
3
s
W
U
6
6
y
8
k
r
P
5
Z
i
b
S
G
b
F
d
q
Q
M
H
8
G
Z
k
J
9
u
3
Z
U
+
x
+
X
t
w
e
g
p
S
m
A
A
9
D
W
5
o
O
r
n
R
O
I
Y
2
H
a
Q
j
/
g
C
K
c
u
T
K
h
K
z
9
U
/
r
r
r
s
F
E
B
B
2
a
Y
6
T
d
R
A
u
A
b
C
3
Q
j
h
F
o
R
w
F
y
D
m
u
h
g
s
Z
r
D
Y
R
i
x
W
E
I
s
t
Y
q
3
o
a
+
a
R
i
z
x
h
P
u
G
e
P
k
z
S
9
Y
p
n
S
9
k
O
q
I
F
M
9
X
z
X
L
G
u
g
i
O
a
1
o
o
q
b
x
q
7
Z
N
7
V
3
C
7
V
f
V
T
u
r
3
+
T
M
u
T
C
z
b
+
L
C
C
n
G
5
f
R
q
2
G
q
R
Q
k
m
u
Z
n
h
B
B
h
w
x
E
S
H
J
6
f
Q
j
2
u
n
t
K
F
R
i
R
/
l
O
Y
i
b
a
M
a
R
t
z
Y
s
z
p
r
Z
A
d
N
o
4
/
M
n
E
I
l
G
4
q
k
Z
r
/
N
F
L
o
B
a
C
r
c
m
N
S
i
3
3
s
Q
a
X
a
q
D
f
y
B
V
Y
d
Z
+
Z
U
r
d
k
6
G
1
S
+
Q
T
/
G
k
4
g
w
i
U
M
k
R
M
9
p
J
L
K
v
2
0
m
K
Q
5
J
9
p
i
A
J
w
m
M
0
J
D
3
t
M
h
Q
R
0
U
/
z
e
0
6
B
l
z
r
i
g
y
D
m
+
m
U
S
5
N
H
F
E
y
m
K
h
J
h
G
n
q
7
U
4
o
3
E
z
V
w
W
X
J
f
r
T
W
T
w
p
p
9
S
l
k
w
k
Y
d
g
A
B
Z
M
Q
y
B
h
k
l
y
b
w
K
S
d
Y
h
l
P
t
I
M
y
p
5
g
r
w
C
O
l
Z
S
H
2
1
Z
i
I
4
N
z
9
5
1
b
l
w
6
0
6
z
3
n
z
f
r
B
6
3
Z
n
L
s
W
M
+
t
F
9
a
+
5
V
i
v
r
W
P
r
n
X
V
m
d
S
x
c
+
l
T
6
U
v
p
e
+
l
H
+
X
P
5
a
/
l
n
+
Z
U
q
3
t
2
Z
n
n
l
l
L
q
/
z
n
L
2
v
M
D
g
E
=
< / l a t e x i t >
/private/var/folders/cy/_n4ycj5x7f9699n7qpqg5c200000gn/T/tpfe886685_8534_4093_a365_a611a8126a93.eps
1980 1985 1990 1995 2000 2005 2010 2015 2020
2
4
6
8
10
12
14
16
18
Data
Model (best fit)
1980 2000 2020
10
10.5
11
11.5
12
12.5
13
Trend
1980 2000 2020
-5
-4
-3
-2
-1
0
1
2
3
4
5
Seasonal cycle
1980 2000 2020
-3
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Noise
⇢ (x, y) =
1
n� 1
nX
i=1
(xi � x̄) (yi � ȳ)
sx sy
R (A,B) =
0
@
⇢ (A,A) ⇢ (A,B)
⇢ (B,A) ⇢ (B,B)
1
A =
0
@
1 ⇢ (A,B)
⇢ (A,B) 1
1
A