Loma Prieta Earthquake Analysis

Open Live Script

This example shows how to store timestamped earthquake data in a timetable and how to use timetable functions to analyze and visualize the data.

Load Earthquake Data

The example file quake.mat contains 200 Hz data from the October 17, 1989, Loma Prieta earthquake in the Santa Cruz Mountains. The data are courtesy of Joel Yellin at the Charles F. Richter Seismological Laboratory, University of California, Santa Cruz.

Start by loading the data.

load quake 
whos e n v

  Name          Size            Bytes  Class     Attributes

  e         10001x1             80008  double              
  n         10001x1             80008  double              
  v         10001x1             80008  double

In the workspace there are three variables, containing time traces from an accelerometer located in the Natural Sciences building at UC Santa Cruz. The accelerometer recorded the main shock amplitude of the earthquake wave. The variables n, e, v refer to the three directional components measured by the instrument, which was aligned parallel to the fault, with its N direction pointing in the direction of Sacramento. The data are uncorrected for the response of the instrument.

Create a variable, Time, containing the timestamps sampled at 200Hz with the same length as the other vectors. Represent the correct units with the seconds function and multiplication to achieve the Hz ( $s^{- 1}$ ) sampling rate. This results in a duration variable which is useful for representing elapsed time.

Time = (1/200)*seconds(1:length(e))';
whos Time

  Name          Size            Bytes  Class       Attributes

  Time      10001x1             80010  duration

Organize Data in Timetable

Separate variables can be organized in a table or timetable for more convenience. A timetable provides flexibility and functionality for working with time-stamped data. Create a timetable with the time and three acceleration variables and supply more meaningful variable names. Display the first eight rows using the head function.

varNames = ["EastWest","NorthSouth","Vertical"];
quakeData = timetable(Time,e,n,v,VariableNames=varNames)

quakeData=10001×3 timetable
      Time       EastWest    NorthSouth    Vertical
    _________    ________    __________    ________

    0.005 sec       5            3            0    
    0.01 sec        5            3            0    
    0.015 sec       5            2            0    
    0.02 sec        5            2            0    
    0.025 sec       5            2            0    
    0.03 sec        5            2            0    
    0.035 sec       5            1            0    
    0.04 sec        5            1            0    
    0.045 sec       5            1            0    
    0.05 sec        5            0            0    
    0.055 sec       5            0            0    
    0.06 sec        5            0            0    
    0.065 sec       5            0            0    
    0.07 sec        5            0            0    
    0.075 sec       5            0            0    
    0.08 sec        5            0            0    
      ⋮

Explore the data by accessing the variables in the timetable with dot notation. (For more information on dot notation, see Access Data in Tables.) Choose the East-West amplitude and plot it as function of the duration.

plot(quakeData.Time,quakeData.EastWest)
title("East-West Acceleration")

Figure contains an axes object. The axes object with title East-West Acceleration contains an object of type line.

Scale Data

Scale the data by the gravitational acceleration, or multiply each variable in the table by the constant. Since the variables are all of the same type (double), you can access all variables using the dimension name, Variables. Note that quakeData.Variables provides a direct way to modify the numerical values within the timetable.

quakeData.Variables = 0.098*quakeData.Variables;

Select Subset of Data for Exploration

Examine the time region where the amplitude of the shockwave starts to increase from near zero to maximum levels. Visual inspection of the above plot shows that the time interval from 8 to 15 seconds is of interest. For better visualization draw black lines at the selected time spots to draw attention to that interval. All subsequent calculations involve this interval.

t1 = seconds(8);
t2 = seconds(15);
xline(t1,LineWidth=2)
xline(t2,LineWidth=2)

Figure contains an axes object. The axes object with title East-West Acceleration contains 3 objects of type line, constantline.

Store Data of Interest

Create another timetable with data in this interval. Use timerange to select the rows of interest.

tr = timerange(t1,t2);
quakeData8to15 = quakeData(tr,:)

quakeData8to15=1400×3 timetable
      Time       EastWest    NorthSouth    Vertical
    _________    ________    __________    ________

    8 sec         -0.098       2.254          5.88 
    8.005 sec          0       2.254         3.332 
    8.01 sec      -2.058       2.352        -0.392 
    8.015 sec     -4.018       2.352        -4.116 
    8.02 sec      -6.076        2.45        -7.742 
    8.025 sec     -8.036       2.548       -11.466 
    8.03 sec     -10.094       2.548          -9.8 
    8.035 sec     -8.232       2.646        -8.134 
    8.04 sec       -6.37       2.646        -6.566 
    8.045 sec     -4.508       2.744          -4.9 
    8.05 sec      -2.646       2.842        -3.234 
    8.055 sec     -0.784       2.842        -1.568 
    8.06 sec       1.078       2.548         0.098 
    8.065 sec       2.94       2.254         1.764 
    8.07 sec       4.802        1.96          3.43 
    8.075 sec      6.664       1.666         4.998 
      ⋮

Visualize the three acceleration variables on three separate axes. To create a tiled layout with the three plots in a 3-by-1 tile arrangement, use the tiledlayout function.

tiledlayout(3,1)
% Tile 1
nexttile
plot(quakeData8to15.Time,quakeData8to15.EastWest)
ylabel("East-West")
title("Acceleration")

% tile 2
nexttile
plot(quakeData8to15.Time,quakeData8to15.NorthSouth)
ylabel("North-South")

% Tile 3
nexttile
plot(quakeData8to15.Time,quakeData8to15.Vertical)
ylabel("Vertical")

Figure contains 3 axes objects. Axes object 1 with title Acceleration, ylabel East-West contains an object of type line. Axes object 2 with ylabel North-South contains an object of type line. Axes object 3 with ylabel Vertical contains an object of type line.

Calculate Summary Statistics

To display statistical information about the data use the summary function.

summary(quakeData8to15)

quakeData8to15: 1400x3 timetable

Row Times:

    Time: duration

Variables:

    EastWest: double
    NorthSouth: double
    Vertical: double

Statistics for applicable variables and row times:

                  NumMissing         Min             Median             Max               Mean              Std      

    Time              0           8 sec            11.498 sec        14.995 sec        11.498 sec        2.0214 sec  
    EastWest          0           -255.0940           -0.0980          244.5100            0.9338           66.1188  
    NorthSouth        0           -198.5480            1.0780          204.3300           -0.1028           42.7980  
    Vertical          0           -157.8780            0.9800          134.4560           -0.5254           35.8554

Additional statistical information about the data can be calculated using varfun. This is useful for applying functions to each variable in a table or timetable. The function to apply is passed to varfun as a function handle. Apply the mean function to all three variables and output the result in format of a table, because the time is not meaningful after computing the temporal means.

mn = varfun(@mean,quakeData8to15,OutputFormat="table")

mn=1×3 table
    mean_EastWest    mean_NorthSouth    mean_Vertical
    _____________    _______________    _____________

       0.9338           -0.10276          -0.52542

Calculate Velocity and Position

To identify the speed of propagation of the shockwave, integrate the accelerations once. Use cumulative sums along the time variable to calculate the velocity of the wave front.

edot = (1/200)*cumsum(quakeData8to15.EastWest);
edot = edot - mean(edot);

Next perform the integration on all three variables to calculate the velocity. It is convenient to create a function and apply it to the variables in the timetable with varfun. In this example, the function is included at the end of this example and is named velFun.

vel = varfun(@velFun,quakeData8to15)

vel=1400×3 timetable
      Time       velFun_EastWest    velFun_NorthSouth    velFun_Vertical
    _________    _______________    _________________    _______________

    8 sec           -0.56831             0.44642             1.8173     
    8.005 sec       -0.56831             0.45769              1.834     
    8.01 sec         -0.5786             0.46945              1.832     
    8.015 sec       -0.59869             0.48121             1.8114     
    8.02 sec        -0.62907             0.49346             1.7727     
    8.025 sec       -0.66925              0.5062             1.7154     
    8.03 sec        -0.71972             0.51894             1.6664     
    8.035 sec       -0.76088             0.53217             1.6257     
    8.04 sec        -0.79273              0.5454             1.5929     
    8.045 sec       -0.81527             0.55912             1.5684     
    8.05 sec         -0.8285             0.57333             1.5522     
    8.055 sec       -0.83242             0.58754             1.5444     
    8.06 sec        -0.82703             0.60028             1.5449     
    8.065 sec       -0.81233             0.61155             1.5537     
    8.07 sec        -0.78832             0.62135             1.5708     
    8.075 sec         -0.755             0.62968             1.5958     
      ⋮

Apply the same function velFun to the velocities to determine the position.

pos = varfun(@velFun,vel)

pos=1400×3 timetable
      Time       velFun_velFun_EastWest    velFun_velFun_NorthSouth    velFun_velFun_Vertical
    _________    ______________________    ________________________    ______________________

    8 sec                2.1189                    -2.1793                    -3.0821        
    8.005 sec            2.1161                     -2.177                    -3.0729        
    8.01 sec             2.1132                    -2.1746                    -3.0638        
    8.015 sec            2.1102                    -2.1722                    -3.0547        
    8.02 sec              2.107                    -2.1698                    -3.0458        
    8.025 sec            2.1037                    -2.1672                    -3.0373        
    8.03 sec             2.1001                    -2.1646                    -3.0289        
    8.035 sec            2.0963                     -2.162                    -3.0208        
    8.04 sec             2.0923                    -2.1592                    -3.0128        
    8.045 sec            2.0883                    -2.1564                     -3.005        
    8.05 sec             2.0841                    -2.1536                    -2.9972        
    8.055 sec              2.08                    -2.1506                    -2.9895        
    8.06 sec             2.0758                    -2.1476                    -2.9818        
    8.065 sec            2.0718                    -2.1446                     -2.974        
    8.07 sec             2.0678                    -2.1415                    -2.9662        
    8.075 sec             2.064                    -2.1383                    -2.9582        
      ⋮

Notice how the variable names in the timetable created by varfun include the name of the function used. It is useful to track the operations that have been performed on the original data. Adjust the variable names back to their original values using dot notation.

pos.Properties.VariableNames = varNames;
vel.Properties.VariableNames = varNames;

Plot the three components of the velocity and position for the time interval of interest. Plot velocity and position in a 2-by-1 tile arrangement.

tiledlayout(2,1)

% Tile 1
nexttile
plot(vel.Time,vel.Variables)
legend(vel.Properties.VariableNames,Location="bestoutside")
title("Velocity")

% Tile 2
nexttile
plot(pos.Time,pos.Variables)
title("Position")

Figure contains 2 axes objects. Axes object 1 with title Velocity contains 3 objects of type line. These objects represent EastWest, NorthSouth, Vertical. Axes object 2 with title Position contains 3 objects of type line.

Visualize Trajectories

The trajectories can be plotted in 2D or 3D by using the component value. This plot shows different ways of visualizing this data.

Begin with 2-dimensional projections. Here is the first with a few values of time annotated.

figure
plot(pos.NorthSouth,pos.Vertical)
xlabel("North-South")
ylabel("Vertical")
% Select locations and label
nt = ceil((max(pos.Time) - min(pos.Time))/6);
idx = find(fix(pos.Time/nt) == (pos.Time/nt))';
text(pos.NorthSouth(idx),pos.Vertical(idx),char(pos.Time(idx)))

Figure contains an axes object. The axes object with xlabel North-South, ylabel Vertical contains 5 objects of type line, text.

Use plotmatrix to visualize a grid of scatter plots of all variables against one another and histograms of each variable on the diagonal. The output variable Ax, represents each axes in the grid and can be used to identify which axes to label using xlabel and ylabel.

figure
[S,Ax] = plotmatrix(pos.Variables);

for ii = 1:length(varNames)
    xlabel(Ax(end,ii),varNames{ii})
    ylabel(Ax(ii,1),varNames{ii})
end

MATLAB figure

Plot a 3-D view of the trajectory and plot a vertical line at every tenth position point. The spacing between vertical lines indicates the velocity.

step = 10;
figure
plot3(pos.NorthSouth,pos.EastWest,pos.Vertical,"r")
hold on
plot3(pos.NorthSouth(1:step:end),pos.EastWest(1:step:end),pos.Vertical(1:step:end),"|")
hold off
box
axis tight
xlabel("North-South")
ylabel("East-West")
zlabel("Vertical")
title("Position")

Figure contains an axes object. The axes object with title Position, xlabel North-South, ylabel East-West contains 2 objects of type line. One or more of the lines displays its values using only markers

Supporting Functions

Functions are defined below.

function y = velFun(x)
    y = (1/200)*cumsum(x);
    y = y - mean(y);
end