Arnoud Buzing
Visualizing the Recent Yellowstone Earthquakes
January 1, 2009
Arnoud Buzing, Quality Assurance Manager

A couple of days ago I read about an unusual “swarm” of earthquakes at Yellowstone National Park. After reading up on this topic a bit (and determining that my home state of Illinois would not be obliterated immediately by a supervolcano outburst), I decided to make an animation about it in Mathematica. First I searched for “yellowstone map sdts” on Google and downloaded this geological map of Yellowstone from the U.S. Geological Survey website. After uncompressing the zip file, I simply pointed Import to the top directory containing the SDTS files:

In[1]:= yellowstone=Import["sdts","Graphics"]

The resulting graphic contains a lot of distracting detail, so I decided to extract just the polygons and give them a muted gray background color. What remains are the outline polygons for each geological layer as observed by the USGS. Also, I set the image size to 1280×720, which makes it suitable for a 720p high-definition video stream:

In[2]:= yellowstone=Graphics[{EdgeForm[Gray],FaceForm[GrayLevel[.8]],Cases[yellowstone,_Polygon,\[Infinity]]},ImageSize->{1280,720}];Show[yellowstone,ImageSize->400]

Next I obtained the full list of earthquakes from the University of Utah Seismographic Stations website. For this I simply copied and pasted the table part of this page and used:

In[3]:= quakes=ImportString["data...","Table"];

to obtain a Mathematica table of the quake data. Another way to do this would be to import the entire web page, but in this case copying and pasting was quick and simple. A typical line with earthquake information came out like this:

{2.8`,"2008/12/31","08:02:11","44.540N","110.370W",0.1`,60,"km","(37","mi)","ESE","of","West","Yellowstone,","MT"}

The data points of interest to me were the magnitude, latitude, and longitude, which in this table are contained in columns 1, 4, and 5 respectively. The latitude and longitude columns also needed to be turned into numbers, which I achieved by dropping the N and W (the last character) from the string and applying ToExpression on the remainder. Finally, I needed to have all the earthquakes in chronological order, which was the reverse of the order given on the web page:

In[4]:= quakes=Reverse[Map[{ToExpression[StringDrop[First[#],-1]],ToExpression[StringDrop[#[[2]],-1]],Last[#]}&,quakes[[All,{4,5,1}]]]];

Now a typical single data point is simply latitude, longitude, and magnitude:

{44.54`, 110.37`, 2.8`}

Because the imported geological map was using a coordinate system based on the transverse Mercator mapping, I had to write a conversion function to take the earthquake latitude and longitude and convert them to coordinates that could be projected onto the map. Several new functions in Mathematica 7 made this relatively straightforward. Most of the conversion options were provided by the Import command, when using the “Projection” element for import:

In[5]:= projection=Import["sdts","Projection"]

In[6]:= convert[{lat_,long_}]:=First[GeoGridPosition[GeoPosition[{lat,-long}],projection]]

For example, this maps directly onto the geological map background:

In[7]:= convert[{44.54,110.37}]

Now it’s time to actually combine the map with the earthquake data. For this I wrote a helper function that shows the background map, the current earthquake, and a list of previous earthquakes:

In[8]:= show[{lat_,long_,mag_},rest_]:=Show[{yellowstone,Graphics[Table[{Black,Disk[convert[Take[r,2]],400]},{r,rest}]],Graphics[Table[{RGBColor[i/4,0,0],Disk[convert[{lat,long}],250 i mag]},{i,4,1,-1}]]}]

(At this point one could use ListAnimate to make an animation, but I was interested in creating a movie I could easily share on a website like Vimeo or YouTube.)

The last step was to export every earthquake as a keyframe PNG image. I used IntegerString[i,10,4] as a convenient built-in function to generate the 0001, 0002, 0003, … sequence. For simplicity in this post I am leaving out the generation of the initial and final key frames:

In[9]:= Do[Print[i]; Export[FileNameJoin[{$HomeDirectory,"quakes/test"<>IntegerString[i,10,4]<>".png"}],show[Part[quakes,i],Take[quakes,i-1]]], {i,215}]

With the PNG files all numbered and saved into a single directory, I used an external tool, FFmpeg, to generate a movie. The -r option specifies the frame rate and the -sameq option preserves the high quality 1280×720 resolution the best:

In[10]:= Run["ffmpeg -r 5 -sameq -i test%04d.png out.mov"]

Earthquake movie--Click to view full movie

This makes a pretty interesting animation. The initial quakes are clustered close together and their intensity ebbs and flows. Then the occasional outliers in the immediate area start to appear. After a while, an elongated ellipse-shaped cluster of earthquakes starts to appear. At the end of this sequence a number of far outliers in the northwest part of the park appear as well, which seem unrelated to the main swarm.

I don’t know enough about seismology to be able to make a conclusion from all of this. But it’s great to be able to take this kind of data and so quickly be able to visualize it. Plus, in my role as Quality Assurance Manager for Mathematica, it’s good to get a chance to test out so many new features in Mathematica 7 in a real-life setting and see how different parts of the system integrate with one another: SDTS import, data imported from tabulated web pages, geodetic conversions, graphic generation, PNG image export, and finally calling an external command to construct a movie.