EGB13

Data on thousands of earthquakes in the New Madrid Seismic Zone since 1974 is available from the Center for Earthquake Research and Information at the University of Memphis.  Most of those were too small to feel.  But they’re real interested in that zone at Memphis because they’re awfully close to it, and it has unleashed some monster quakes in the past.

The data includes the estimated depth of the quakes as well as the geographic coordinates.  Below is my first effort to use that data to visualize the fault zone in 3 dimensions.  There is a lot of data that gets in its own way, so an animation that tilts and rotates was chosen as a a means to see it from different angles.

In order to view this animation you need a modern browser that supports HTML5 videos in Ogg Theora format.  The current versions of FireFox and Chrome work., though if you insist there is another version on YouTube.

I’m fairly well pleased with how it turned out, but I’m concerned the reliability of some of the data.  That is discussed below along with notes on how it was produced.

Depth data

Historgram of depths of 8,557 central US quakes to nearest kilometer

Did you notice that there seems to be a plateau of sorts in the data at exactly 5 km?  That’s suspicious.  Likewise the lighter grouping at 10 km.   I can take a bucket of data and manipulate it into a visualization.  But I’m afraid that I don’t know enough about this particular data, how it is collected, and what issues there are in its accuracy and reliability.  The page titled Description of Catalog Data doesn’t help much in that regard, but it does make me think I should comb the data I’m using for zero depths.

So, as interesting and useful as this looks, I don’t think I’m going anywhere with it until I know why the data bunches up at those two depths.  If anyone comes across this who can tell me how to cull potentially unreliable entries from the dataset, I’d love to hear from you.

If that data is true, then those groupings at 5 and 10 km depths would be very interesting, indeed.

Update: Here is a histogram of depths rounded to the nearest kilometer.  Remove data at zero, 5, and 10, and the smoothing curve will be even lower than shown.

Historgram of depths of 8,557 central US quakes to nearest 0.1 kilometer

Update #2: The curator of this data kindly responded to a query, and he reports that the USGS normally fixes the depth at 5 km in order to get a more accurate epicenter and origin time. The data is from a variety of sources, though, and that’s why it doesn’t all describe a level plain right at 5 km.

He also states that the USGS sometimes fixes the depth of events at 10 km, but normally not in this part of the country. This second histogram rounds all depths to 0.1 km. There is an obvious spike at exactly 10.0 that certainly suggests artificial coersion of some class of data; if it were due to actual activity at around 10 km, it would more likely create a hump in the histogram including 9.8, 9.9, 10.1, 10.2, etc.

Production

I used gnuplot for this.  It’s a great graphing and plotting program, but there is precious little scripting power in its language.

Scripts

One file sets up the environment and then calls a second file to render a frame, adjust the view, and then call itself to render the next.  It iterates until the z axis has rotated a full 360°.  That’s a pretty crude way to build a loop.

First, file “nm.gp”:

#!/usr/bin/gnuplot

set        xlabel "Easting (km)\nUTM 15N"
set        ylabel "Northing (km)\nUTM 15N"
set        zlabel "Depth (km)"
set xrange [751:849]
set yrange [3960:]
set zrange [-30:0]
set size 0.67,1
set grid
xrot = 0
zrot = 0
xdelta = 0.5
zdelta = 0.5
call  "nm2.gp" 1
exit

The size is set as it is to induce the x and y axes into compatible scales so the map comes out right — gnuplot isn’t really a GIS program.  The parameter to nm2.gp is an iteration count used to build the frame’s file name.  There is no scoping for variables in this utility, so that could have been done with “global” variables in the manner of xrot and zrot.  Here is file “nm2.gp”:

#!/usr/bin/gnuplot

iter = $0
ofile = sprintf("frames/nm-%d.png", iter)

if (xrot <= 85) xrot2 = xrot ; else xrot2 = 85 ;
print "xrot=", xrot2, " xdelta=", xdelta, " zrot=", zrot, " zdelta=", zdelta, " iter=", iter, " ofile=", ofile
set        view xrot2, zrot, 1.12268, 0.84487

splot 'quakes1.txt'     using 1:2:3 title "" with dots 11       \
, 'quakes2.txt'         using 1:2:3 title "" with dots 18       \
, 'mo3.txt'             using 2:3:4 title "" with lines 30      \
, 'tn1.txt'             using 1:2:3 title "" with lines 30      \

set terminal png small
set output ofile
replot

if (zrot < 90) xrot = xrot + xdelta
if (zrot > 270) xrot = xrot - xdelta
if (xrot < 0) xrot = 0
if (xrot > 360) xrot = 360

zrot = zrot + zdelta
iter = iter + 1

if (zrot <= 360) call  "nm2.gp" iter

exit

The rotation about the x-axis (tilt), stored in xrot, increases with each step as the vertial rotation about the z-axis, in zrot, increases from 0° to 90°, then stays steady until the z rotation reaches 270°.  The tilt then decreases steadily for the last 90° until the starting point is reached.  Note xrot2 implements a ceiling function limiting the actual tilt to 85°.  This is pretty simple stuff.

The work of drawing the frames is done with the splot function, of course.  The file “quakes1.txt” contains all quakes below magnitude 2, and “quakes2.txt” contains data on all quakes magnitude 2 and above.  The larger magnitude quakes are plotted with a darker dot.  Data in a gnuplot-suitable format only for those sections of the state lines in the constrained view is in the files “mo3.txt” and “tn1.txt”.

(You may have noticed that I put the comma after the line break between data sources in the plot command, not before.  If that looks odd, look at it again and see how it makes adding or removing data sources during the iterative development process more straightforward.  A fifth source can be added, or the last one can be moved higher up, without modifying any of the lines’ contents.  I don’t recall where I picked that up, but I like it and use similar formatting in Perl and C, etc.)

Data

All geographic data is in UTM coordinates for zone 15 North with units of kilometers instead of meters.  That makes the range of numbers on the x and y axes have reasonable values.  The original data has geographic coordinates and had to be projected with this command:

proj +proj=utm +zone=15 +ellps=GRS80 +datum=NAD83 +units=km +no_defs

That’s identical to EPSG:26915 with the units set to km.  In that coordinate system, only quakes within these geographic limits were included:

Lower left: x=748.6 y=3963.52

Upper right: x=849.06 y=4079.34

That cooresponds to these geographic coordinates, though the projected area doesn’t describe a level rectangle in geographic coordinates:

90.21266°W 35.7538°N
89.13400°W 36.8272°N

Drawing the lines for Missouri and for Tennessee at this location actually gives us parts of Missouri, Arkansas, Tennessee, and Kentucky.  I wrote a little Perl script to read shapefiles containing line or polygon data and print the line segements out in a text format suitable for gnuplot.  A lot of manual work with sort and vi resulted in the two files containing just those parts of the states’ boarders in the area to be plotted.  The states lines were duplicated then with a depth of 30 km to create a map reference below the quake data points as well as at the surface.

Rendering

The gnuplot script(s) generated 730 frames, but left too much empty space at the top, so that was cropped off with ImageMagick:

mogrify -crop 428x$((480-60))+0+60 frames/nm-*.png

Those were rendered into an MPEG-2 video:

ffmpeg -r 20 -i frames/nm-%d.png -vcodec mpeg2video -r 20 -b 6000k -an nm.mpg

That was then converted into an Ogg-Theora video:

ffmpeg2theora -o nm.ogg --videobitrate 1200 --optimize nm.mpg

I’d rather have ffmpeg render the ogg file directly, but the version I have doesn’t do that, apparently.  Some loss of fidelity is involved in doing it this way.  That’s why I render the original video file at a rather high bitrate for its size.