I’m a GPS addict. I have a handheld GPS, a computer GPS, a GPS phone, two GPS watches, two GPS cameras, and maybe some others. Wherever I go, chances are pretty good I have at least one GPS with me. Anytime I run/bike/hike/walk/ski I keep a record of where I went using GPS. Now that I have all this data I want to make use of it in a meaningful way. Mathematica is a fantastic tool to analyze all my geographic data.

Here’s an example from a recent trail run I did at a nearby park. The data is stored in a GPX (GPS exchange format) file, which is a specific type of XML. We can bring the data into Mathematica using Import.

Now we can extract all the track points (the XML elements named “trkpt”).

Each of these track points has a latitude, longitude, elevation, and time.

The first thing I’d like to do is make an elevation profile plot. This plot will show the distance on the x axis and the elevation on the y axis. The elevation values are measured quantities that come straight from the GPX file, while the distance values are quantities we can compute from the latitude/longitude pairs.

Let’s pull out all the latitudes, longitudes, and elevations.

Next, let’s partition the list of points into a list of adjacent pairs of points. Then we can compute the distance between the two points in each pair using the GeoDistance function. GPS latitude and longitude coordinates are measured in the WGS-84 reference datum, so we’ll pass that information along to GeoPosition and GeoDistance to make sure the computed distances are as accurate as possible.

Now we have a list of relative distances between each successive point. Let’s convert these relative distances to absolute distances from the beginning of the track to each point by summing up all the previous relative distances.

The resulting distance and elevation values are in meters, but in the United States we are accustomed to measuring distance in miles and elevation in feet, so let’s convert the units using the Convert function.

Now that the data is ready, we can finally make our elevation profile plot.

Those three 100-foot hills in the last couple of miles were tough.

So that was interesting, but what I’m really more interested in is seeing my GPS track on a map. This will be a two-step process. First, we need to draw a suitable map. Second, we need to overlay the GPS track on top of it.

For the map I’m going to use USGS aerial photography from TerraServer, a site run by Microsoft that provides a web service to download USGS images. We’ll use InstallService from the WebServices` package to make it easy to get the images we need.

We have a rectangle (the bounding box of the GPS track) that we’d like to map, so let’s determine what that rectangle is and ask TerraServer for the AreaBoundingBox. This will tell us which images we need to download. For good measure we’ll go ahead and make sure the background is at least 10% larger than the track on each side.

Now we need to create a two-dimensional grid of map tile images. Let’s figure out which tiles we need to request.

We’ll get a sample TileId value for later modification.

Now we can download the tiles in the specified X and Y ranges using the GetTile web service method.

We can easily put these tiles together with ImageAssemble.

In addition to aerial photography, TerraServer also provides topographic maps.

Now comes the challenging part. We have GPS points given in latitude and longitude (3D spherical coordinates). We have an image drawn in some arbitrary 2D Cartesian coordinate system. How do we correctly combine these? Well, we need to make sure our latitude and longitude values are projected from 3D spherical to 2D Cartesian coordinates using the same method as the source images.

We can see from the TerraServer web service API documentation that it converts latitude and longitude values to 2D Cartesian coordinates using the Universal Transverse Mercator coordinate system. TerraServer also provides a function to do the conversion for you. This function would work quite well for a handful of points, but this GPS track has thousands of points. It would take an awfully long time and put an unnecessary burden on TerraServer’s resources to make this API request every time we want to convert a single point.

Fortunately, this conversion is well defined, so we can implement it directly in Mathematica. Unfortunately, the code is a little messy and beyond the scope of this blog post. Here’s the magical function:

I should note that the image tiles we downloaded from TerraServer cover an area equal to or greater than that which we requested. They almost certainly cover a greater area, so let’s figure out the area they do cover.

Now we’ll convert these to UTM.

The image itself isn’t in UTM coordinates, it’s in some other coordinate system. So we will need to scale all UTM values to the image’s coordinate system. The next step is to figure out what scaling factors to use.

With this information we can now define a single function that will convert latitude and longitude values into the coordinate space of our image.

Now we can map this function across all of our GPS points.

Finally, we can show the track on top of the map image.

The map is good, but not perfect. Let’s trim the image down to the original bounding rectangle we requested from TerraServer.

We can even make it interactive.

Here’s the track on a topographic map.

Let’s try a more elaborate example to see if we can take this mapmaking to the next level.

Last November my wife ran the Indianapolis Monumental Marathon. I rode my bike to several different locations along the course to cheer her on and take photos. We both wore GPS devices to record our tracks. The funny thing about this was that she ran faster than either of us expected and I missed her several times because she arrived at a location on the course before I did. Let’s find out what this looked like.

For this example we’ll be dealing with two GPS tracks, the Rob one and the Melissa one. Let’s begin by importing the GPX files and extracting the latitudes and longitudes from the track points.

Next, we’ll get the map. Let’s take all the code we used earlier and stick it into a function to make it easier to reuse.

This function returns the map image and a pure function that can be used to translate latitude/longitude pairs into the coordinate space of the map.

Let’s take a look.

Now we have a map of the race course (in blue) and my path on the bike (in red). This would make a neat interactive map, but we can’t just use the same trick as before where we just added another track point to the line for each step of the animation, because the two tracks won’t be properly synchronized.

To solve the synchronization problem, we’ll need to get the original time stamp of each track point from the GPX data.

Next we will build an InterpolatingFunction that when given a time will return the approximate location for the track. Since the track points are spaced one or more seconds apart, any time between track points will return a value interpolated between the two nearest points.

We need to determine the earliest and latest times of either track.

Now we can make the interactive map.

Because the course is almost entirely north/south, this map has a very inconvenient aspect ratio. It’s much too tall and not very wide. We could improve it by focusing on a small section (with a reasonable aspect ratio) at any given time. So let’s zoom in a bit and crop the image.

Here we have the cropped PlotRange and ImageSize, as well as the larger, tile-aligned PlotRange and ImageSize. Let’s keep track of the larger size so we have a larger work area.

Let’s choose an arbitrary size for the final image.

As we want the focus at any given time to be the latest points in the two tracks, let’s center the image at the midpoint between the latest track points from each track.

We can just as easily export this as a sequence of images and assemble them into a movie using QuickTime Player (with QuickTime Pro features enabled).

I counted at least four times when I (red dot) arrived at a location on the marathon course after my wife (blue dot) had already passed. Next time I may have to tell her to slow down.

These are just a few ways to analyze GPS data with Mathematica. We could just as easily plot speed or pace. We could request map images from other services. We could zoom in or out on the map. Mathematica‘s integration of symbolic computation, graphics, import/export, web connectivity, and computable data (e.g. GeodesyData) make it the ideal tool for these types of data explorations.

Download this notebook