Adding FrameTicks to Maps with Mathematica

As I’ve illustrated in a number of posts, Mathematica is developing a number of useful and relatively easy way to construct maps. One of the most frustrating issues with the implementation is that the default latitude tick marks are not in degrees (along with the ever-present problem that they have no idea how long a tick mark should be in a quality plot). In this post I describe a simple function to create tick marks for the Mercator projection maps, and the code is relatively easily generalized to work with other projections.

GeoGraphics – Incorrect Ticks Marks

Let’s start with a simple map of the location of the 2015 Nepal (Gorkha) Earthquake.

(* set up the earthquake information *)
cat = {607208674, "CSEM", "2015-04-25", "06:11:26.30", 28.28, 84.79, 
       "     ", "", "GCMT", "MW", 7.9};
eloc = {#[[5]], #[[6]], #[[11]]} & /@ {cat}
eqradius[m_] := N[10^(0.5*m - 2.25)]
scaleFactor = 1;
(* plot the map *)
GeoGraphics[{Polygon[Entity["Country", "Nepal"]], 
  GeoStyling[Opacity[0.25], EdgeForm[Black], FaceForm[Red]],
  GeoDisk[{#[[1]], #[[2]]}, (scaleFactor*Quantity[eqradius[#[[3]]], "km"]) ] & /@ eloc},
  GeoRange -> {{26.25, 29.75}, {82.5, 88.5}}, Frame -> True,
  FrameTicks -> {Automatic, Automatic},
  BaseStyle -> {18, FontFamily -> "Helvetica"},
  PlotLabel -> Style["Gorkha, Nepal Earthquake, 2015-04-25 06:11:26.30 UTC, Mw 7.9", {16}],
  GeoZoomLevel -> 8, FrameStyle -> AbsoluteThickness[1.2],
  GeoScaleBar -> {"Imperial", "Metric"},
  ImageSize -> 600] // Print

Here is the output

Mathematica map of the 2015 Nepal earthquake epicenter plotted with the default FrameTicks. Those are not latitudes on the vertical axis.

The map looks pretty good, but those are not latitude tick marks. Worse yet, in this example, they happen to be close enough to actual latitudes that they could fool you. The latitude of the center of the red circle is 28.28N, the tick marks lead you to believe the location is at roughly 29.5N. Technically, the tick marks are not wrong, they are simply in the practically useless transformation coordinate system. To get the correct tick marks, you need to specify them yourself, convert to projection coordinates, and label appropriately.

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

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