Parsing USGS QuakeML files with Python

The USGS provides earthquake lists in a number of formats including kml for Google Earth, easy-to-read comma-separated-value (csv) files that are quite simple to parse with any spreadsheet, matlab, Mathematica, and many unix tools. Some more information is available in the JSON formatted files, which are the preferred methods for downloading the data feeds. The XML formatted data, stored in the community-standard quakeml dialect, contain the most information. XML is not a human-readable format, but is compatible with the database centric workflows often used in modern earthquake monitoring activities. Accessing the data directly from the databases makes it more likely to have fewer translation errors (or at least consistent data).

My goal in these notes is not to explain the QuakeML (qml) format (I couldn’t if I wanted). You can read about QML at https://quake.ethz.ch/quakeml/, it’s a complex beast of a data format that can be used to store many of the different types of information used in earthquake analyses. Information items range from seismic event hypocenters and origin times (which are collectively called origins) to individual seismic wave arrival time and amplitude measurements. QML contains many schema and allows much of the seismic data to be encoded as they are measured. For example an arrival time pick can be created and stored before the earthquake location associated with that arrival is estimated. The same is true for the seismic wave amplitudes. Storing the information as it is measured is more robust and more conducive for near-real time monitoring processing flows. This helps produce quick locations and magnitudes, which lead in turn to quick earthquake alerts.

At the analysis end of the flow, the information organization can be less convenient than many old formats. For example, information that a seismologist might like to list together may be split apart in theQML file. An arrival time pick and/or measured wave amplitude may be stored as part of the event object, but a travel time residual must be associated with an event origin object. And a single event may have multiple origins – in the simplest case, one origin could be the hypocenter location estimated from body-wave arrival times, another origin could be a rupture centroid estimated using a moment-tensor inversion of long-period waveforms. These files can be complicated. Another issue with USGS QML is that sometimes you get more information, other times you get less. For example, a recent earthquake file might Include information on the station, network, channel, and location code for the waveform used for a measurement, but a file corresponding to an event from a few years older may just provide the station ID. Changing content can have you pulling your hair out when trying to extract information.

QuakeML can also simply be confusing, look at this encoding of an arrival time pick

.
.
.
        <arrival publicID="quakeml:us.anss.org/arrival/200081v8/us_3395_bpid-2001275574">
          <pickID>quakeml:us.anss.org/pick/200081v8/us_3395_bpid-2001275574</pickID>
          <phase>Pn</phase>
          <azimuth>226.495</azimuth>
          <distance>13.8121</distance>
          <takeoffAngle>
            <value>71.6</value>
          </takeoffAngle>
          <timeResidual>1.30</timeResidual>
          <timeWeight>0.6300</timeWeight>
          <creationInfo>
            <agencyID>us</agencyID>
            <agencyURI>smi:anss.org/metadata/agencyid/us</agencyURI>
            <author>manual</author>
            <creationTime>2016-12-17T10:59:33Z</creationTime>
          </creationInfo>
        </arrival>
.
.
.

Look carefully and you will see that the azimuth, distance, and time items all are specific values, the takeoffAngle item has a sublevel to store a value. I don’t know why. Variability in data storage practice adds flexibility, but also makes codes hard to maintain.

Continue reading

Recovering Missing Inkscape Windows

I am evaluating Inkscape (http://inkscape.org) as a replacement for Adobe Illustrator. Inkscape is an open-source, X-Windows-based illustration tool. The interface is nowhere near as nice as that in Illlustrator, but even with steep educational discounts, the Adobe suite costs me roughly $240 per year if I want to have it on my office iMac and my six-year-old Macbook Pro. So far, I have been able to import most of the illustrations that I made with Illustrator into Inkscape, with the exception of some files in which I used some nonstandard fonts. So, I plan to see how it goes without an Illustrator license.

I have two monitors connected to my office iMac, and when I load a file into Inkscape, the main window moves out of view (offscreen). Inkscape is running and the window exists, it just off the screen. You can see it, and even select it, using Mission Control (ctrl-up-arrow), you just can’t do anything with the file. Spending your time finding a recently launched program’s windows is not a good start for any software, but we are talking about an interesting tool, so a work around may be worth a little effort in work-arounds.

Continue reading

Seismic Reflection/Transmission Coefficients with Matlab

A common computation in seismology is the calculation of reflection and transmission coefficients that describe the partitioning of energy when a seismic wave strikes a boundary between elastic materials. The coefficients depend on the seismic wave speeds and densities on either side of the boundary and on the incident wave horizontal slowness, or ray parameter. The formulas for the values can be found in Table 3.1 of Lay and Wallace (1995) (note that an error in the sign of the second term in the computation of b) or in equations 5.38 and 5.39 of Aki and Richards (1980). For a welded elastic boundary, an incident P or SV wave results in four waves, two reflected and two transmitted or refracted waves. Thus for each case we need four coefficients – notation for each is illustrated in the figure below.

Notation used to identify reflection and transmission coefficients. We must consider two cases, that for an incident P wave and that for an incident SV wave. R identifies “Reflected” waves, T represents “Transmitted” or “refracted” waves.

Continue reading

Elevation Profiles with Mathematica

Introduction

One of the activities that I use in my 100-level earthquakes and society class is an exploration of Earth’s topography using Mathematica. Mathematica provides easy access to elevation data and has some relatively simple to use commands for making rudimentary maps. They are not competitive with a nice GMT (Generic Mapping Tools) product, but they are useful for exploring concepts. I keep things simple and have students simply execute the notebook commands to produce topographic profiles along the paths following a constant latitude (e.g. 30N and 30S are interesting). This keeps the commands short and limits their focus on the code. For this note I will use a profile along the latitude of Penn State, which is located near the geographic center of Pennsylvania (latitude is 40.795N).

The Code

We start by creating a reference map that shows the location of the latitudinal profile using a red line. Here’s the Mathematica code:

(* create a map and show where the latitude "slide" is *)
latitude = 40.795;
p0 = Framed[
  GeoGraphics[GeoRange -> "World", GeoProjection -> "CylindricalEquidistant",
  GeoZoomLevel -> 2, GeoBackground -> "ReliefMap", ImageSize -> Full,
  GeoGridLinesStyle -> Directive[Red, Thick], 
  GeoGridLines -> {{latitude}, None}]]

and here is the map:

Mathematica world map with a red reference line at a latitude equal to 40.795N, the latitude University Park, PA.

Continue reading

Converting Between Moment Magnitude and Seismic Moment with Python

Introduction

The historical tradition of using magnitude to communicate the size of an earthquake is deeply entrenched both within seismology and throughout the public. So it makes sense to have a way of mapping seismic moment into a magnitude. Kanamori (1977) developed just such a conversion, which we call moment magnitude and represent as \(M_W\). The relationship was defined in CGS units as:

$$M_w=\frac{2}{3} \log _{10}\left(M_0\right)-10.7~,$$

where \(M_0\) is the moment in dyne-cm. To convert a moment from dyne-cm to N-m (newton meters), divide by \(10^7\). We can invert the equation to convert \(M_W\) to seismic moment:

$$M_0=10^{1.5 \left(M_w+10.7\right)}~.$$

Converting back and forth from the two measures is a common task, and long ago I wrote a simple command-line driven fortran code to perform the conversion. In this note I describe two python scripts that will do the job. Like the original Fortran, both use command-line arguments, and part of the reason for this note is to provide an example of using command-line arguments in python. I use the convenient argparse package to handle all the details.

Python Code

The first version uses two command line arguments -mw and -m0 to indicate whether you want to go from \(M_W\) to \(M_0\) or from \(M_0\) to \(M_W\). Here is the code:

#!/usr/local/anaconda/bin/python
#
#====================================================================
#   MAIN PROGRAM
#====================================================================
if __name__ == '__main__':
    #
    import argparse, math
    #
    parser = argparse.ArgumentParser( \
        description='Convert between seismic moment and moment magnitude.')
    #
    parser.add_argument('floats', metavar='Value', type=float, nargs='+',\
        help='A seismic moment in dyne cm or an Mw value.')
    #
    args = parser.parse_args()
    #
    for value in args.floats:
        if value &lt;= 12:
            print "Mw = %.2f, M0 = %9.3e dyne cm" % (value,math.pow(10,1.5*(value+10.7)))
        else:
            print "M0 = %9.3e dyne cm, Mw = %.2f" % (value,0.66667*math.log10(value)-10.7)
#

Continue reading

Quick Lists of Recent Large Earthquakes with Python

Introduction

I am often interested in a quick summary of recent global seismic activity. A number of good web sites and apps can provide such information, but I often want a simple, quick summary. I wrote the following python script to execute from a command line (as in the OS X Terminal App) and to provide basic information on recent large earthuake activity. I place the script in ~/Dropbox/bin and can run it from all my Macs. I have also installed it on my iPads using the Pythonista App.

The script does nothing fancy, it simply performs a few searches on the U.S. Geological Survey’s earthquake feeds and summarizes information in three categories – great (M ≥ 8) earthquake occurrence in the last three years, major and great (M ≥ 7) earthquake activity in the last year, and strong, major, and great (M ≥ 6) earthquakes in the last 30 days.

Continue reading

Word Clouds with Mathematica

One of the examples in my introductory classes is to create a word cloud from a text file containing the text of one of their reading assignments. A word cloud is a graphical representation of the frequency with which words occur in a section of text. The underlying assumption is that the more often a word is used, the more important it is to the topic. To meet that  assumption, we have to remove common words such as “the”, “and”, “a”, etc. as well as combine different forms of a word such as plurals and singular versions – at least for the important terms.

The example used below includes the text from an introductory chapter on Plate Tectonics, that is assigned early in the course. The overall goal is to show students a use of computational thinking (and computing) that helps in the analysis of something other than mathematical computations.

Before making the word cloud we have to scrub the text to remove non-words such as numbers, and to remove words that carry little information about the topic such as “the” – Mathematica calls these terms stop words. I use the following code to make a word cloud and then adapts some of the scrubbing functions to better clean the specific text that I am processing.

Import the text file

I created the text file using the scanned images of the textbook and the Mac OSX tool pdfPenPro. I copied and pasted the text into an editor and saved the file. We start by reading
the text into a string variable called d0, which we can use to manipulate the text.

(* you will have to replace the path to the text file using 
   the the menu command: Insert->File Path...
   We load and store the text into an item we'll call d0 *)
d0 = Import["Chapter02.txt"];

Continue reading

Working with Slab Model 1.0 in Mathematica

The GeoGraphics features in Mathematica that were introduced some time around Version 9 and later allow some useful exploration of geophysical data, models, etc. In this note I illustrate how to read and display (in a simple form) the tectonic Slab 1.0 Model of Hayes et al [2012].

Slab Model 1.0 is available at https://earthquake.usgs.gov/data/slab/, and was constructed to be able to compare earthquake locations and faulting geometry with historical spatial patterns of seismicity, which are believed to outline the shapes of recently subducted lithospheric slabs.

You can read about the model in Hayes, G. P., D. J. Wald, and R. L. Johnson (2012), Slab1.0: A three-dimensional model of global subduction zone geometries, J. Geophys. Res., 117, B01302, doi:10.1029/2011JB008524, and if you use the model, be sure to cite that paper.

Parsing the Slab 1.0 kmz File

The simplest way to import the model is to parse the kmz file available from the US Geological Survey’s site (see the link above). You can read about the kmz file format here. Essentially kmz is a compressed form of kml, which is short for keyhole markup language.

KMZ files can contain a number of constructs and in this case, the slab contours are stored in graphics primitives for drawing lines. What we ultimately want to extract from the file are the latitudes and longitudes of contours of slab depth. The contours are stored in increments of 20 km, starting at 20 km beneath sea level.

The information can be imported into Mathematica using the standard Import command. We can see what’s in the kmz file by asking Import for the elements of the file.

Import["/usr/local/data/slab_model_1.0.kmz", "Elements"]

Out[97]= {"Data", "Graphics", "GraphicsList", "ImageOverlays", "LayerNames", 
    "LayerTypes", "SpatialRange"}

We can import the kmz data using

slabkmzdata = Import["/usr/local/data/slab_model_1.0.kmz","Data"];

where you will of course have to change the path name to match that of your kmz file.

The list slabkmzdata contains a list of the graphics primitatives in the kml file. We need to unpack the list and pull the graphics commands apart. I’ll do it with some short simple functions to extract what we need in a way that makes it easy to display (to check it correct) and the manipulate.

Continue reading