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.

Some Slab Model Convenience Lists and Functions

Each slab is in its own layer, so you can get a list of slabs using the “LayerNames” option.

slabnames = Import["/usr/local/data/slab_model_1.0.kmz", "LayerNames"]

Out[99]= {"alu", "van", "sum", "sol", "sco", "sam", "ryu", "phi", "mex", "kur",
 "ker", "izu"}

The function GetSlabContours lets you extract the contours for a particular slab (asked for by name) from the kmz data. The function depends on the slabnames list defined immediately above. The contour interval in the files is 20 km and the shallowest contour corresponds to 20 km deep.

GetSlabContours[slabname_, slabkmzdata_] := 
 Module[{slabindex, slablines, slabcontours},
  slabindex = Flatten[Position[slabnames, slabname]][[1]];
  slablines = "Geometry" /. slabkmzdata[[slabindex]];
  slabcontours = Cases[#, _GeoPosition, Infinity] & /@ slablines
  ]

Extracting the contour depth took some experimentation to discover the index level required. I ended up using the four-level index [[1]][[1]][[1]][[3]]. A more elegant command probably exists, but this works, so I’ll just keep it.

The next functions provides a set of Mathematica GeoGraphics commands to draw the contours in Mathematica. In the first, I map the colors using the Hue command and normalize the depths by a value of 700 km, which is near the depths of the deepest earthquakes. In the second the contours are all black, but contours that are multiples of 100 km are shown as thicker lines. The input for either of these functions is a list of slab contours that you obtain from GetSlabContours to either of these functions.

(* Returns graphics commands for plotting color slab contours *)
GetSlabGeoGraphicsColor[slabcontours_] := Module[{},
  {{ColorData["BrightBands"][-#[[1]][[1]][[1]][[3]]/700], Line[#]} & /@ slabcontours}
  ]
(* Returns graphics commands for plotting black slab contours with thicker lines
  for each 100 km of slab depth *)
GetSlabGeoGraphicsBlack[slabcontours_] := Module[{},
  {{If[Mod[-#[[1]][[1]][[1]][[3]], 100] == 0, AbsoluteThickness[2], 
       AbsoluteThickness[1.0]], Line[#]} & /@ slabcontours}
  ]

An Example Display

Now that we have those utility functions, we can extract the slab contours by slab name, and then use the basic Mathematica GeoGraphics commands to add the information to a map. Here’s an example for the Middle-American Trench slab, which has the id “mex”.

Extract Contours for a Slab from the kmz Data List

We start by extracting the contours for one of the slabs by name.

mexslab = GetSlabContours["mex", slabkmzdata];

Show Color-Coded Depth Contours for the Slab

The following commands will create the graphics commands and show them on a standard relief map using colors (red for shallow and blue for deep). The contour interval in the files is 20 km.

slabgeographics = GetSlabGeoGraphicsColor[mexslab];
reliefmap = GeoGraphics[{AbsoluteThickness[2.28], slabgeographics}, 
    GeoBackground -> "ReliefMap", ImageSize -> 500]

The output is

Mid-American Trench slab contours from model Slab 1.0. Contour interval is 20 km and red indicates lesser depths, blues greater depths. The shallowest contour corresponds to 20 km below sea level.

Show Thickness Varying Contours for the Slab

The following commands will create the graphics commands and show them on a standard relief map using black lines and thicker lines for contours that are multiples of 100 km. The contour interval in the files is 20 km.

slabgeographics = GetSlabGeoGraphicsBlack[mexslab];
reliefmap =  GeoGraphics[{AbsoluteThickness[2.28], slabgeographics,
    GeoBackground -> "ReliefMap", ImageSize -> 500]

The output is

Mid-American Trench slab contours from model Slab 1.0. Contour interval is 20 km. Thick lines indicate contour that is a multiple of 100 km.

Summary

Mathematica provides a simple way to parse kmz files and to extract the information on latitude and longitude. The Slab 1.0 Model kmz is relatively straight forward since it contains one slab per layer and the slab coordinates (lat, lon, depth) for each contour are stored as kml line objects. During the import, Mathematica converts to these elements into Line objects that can be easily manipulated in Mathematica.

Credits

I used the syntax highlighting tool at https://andrewsun.com/tools/syntax-highlighter/ to typeset the Mathematica code.