Wolfram Computation Meets Knowledge

Visualizing the Recent Yellowstone Earthquakes

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 webpage, 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 webpage:

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"]

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 webpages, geodetic conversions, graphic generation, PNG image export and finally calling an external command to construct a movie.

Comments

Join the discussion

!Please enter your comment (at least 5 characters).

!Please enter your name.

!Please enter a valid email address.

7 comments

  1. You might also be interested to see the more conventional view on my website, but less conventional than usual though!. Of course it does not have the animation that mathematica can give.

    Reply
  2. if you are interested in the new madrid fault zone please visit my site below for more information.

    New Madrid Earthquake

    Reply
  3. You should share your findings with a seismologist as I’m sure your combined knowledge could help us all.

    Reply
  4. Could you plot depth against longitude or latitude? This would give a view to see if the stress was moving vertically.
    Thanks, Brad Johnson

    Reply
  5. Very interesting – from a non-mathematician viewpoint the first thing that I wanted to know after watching the movie was/is: WHERE IS THAT! Sure would be good to have one or two landmarks on the map. Or an interactive cursor which would display the lat/lon in a box as you moved the mouse so we could plug the coordinates into Google Earth. Or, better yet, an option to download a kmz file which would put it as an overlay to Google Earth!

    Is this the routine rumblings of Old Faithful or a fault under the lake or true slippages somewhere that you could see on the surface??? I want to know more!

    Reply
  6. I find this stuff fascinating! I am not a mathematidian or any kind of scientist. But I live in Tucson, Az. I have heard of the earthquakes in Yellowstone and the terrible destruction that could occur if it explodes!! I read that the entire Yellowstone basin has reisen about 8 inches and that new geisers are showing up all over the area? Can you confirm any of this?

    Reply
  7. Looking for prediction of seismic excitation of water basins creating Seiche actions.

    larry.slotta@gmail.com

    Reply