Wolfram Computation Meets Knowledge

Mapping GPS Data

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.

Defining the directory where the data is

Importing the data

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

Extracting the track points

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

Sample trkpts data

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.

Code that pulls out the latitudes, longitudes, and elevations

Code that pulls out 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.

Partitioning the list of points into a list of adjacent pairs of points

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.

Converting relative distances to absolute 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.

Converting meters to miles and feet

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

Elevation plot profile

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.

Installing TerraServer with WebServices`

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.

Gathering the latitudes and longitudes

Increasing the background area

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

Figuring out which tiles to request

Figuring out which tiles to request

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

Getting a sample TileId

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

The downloaded tiles

We can easily put these tiles together with ImageAssemble.

Assembling the tiles

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

A topographic map downloaded from TerraServer

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:

Function to convert latitudes and longitudes to 2D Cartesian coordinates

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.

Figuring out the area the tiles cover

Figuring out the area the tiles cover

Now we’ll convert these to UTM.

Converting to UTM

Converting to UTM

Converting 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.

Figuring out what scaling factor to use

Figuring out what scaling factor 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.

Single function that converts latitude and longitude values

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

Mapping the function across all GPS points

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

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.

Trimming the image to the bounding rectangle

Trimming the image to the bounding rectangle

Trimming the image to the bounding rectangle

We can even make it interactive.

Adding Manipulate to make the image interactive

Here’s the track on a topographic map.

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.

Importing the GPX files

Extracting the latitude and longitude points

Extracting the latitude and longitude points

Extracting the latitude and longitude 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.

Putting the code into a function

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.

Map image and pure function

Let’s take a look.

A map with two paths marked in red and blue

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.

Getting the time stamp from each point

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.

Building an InterpolatingFunction

Building an InterpolatingFunction

Building an InterpolatingFunction

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

Determining the earliest and latest times

Now we can make the interactive map.

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.

Zooming in and cropping 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.

Tracking the larger size

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

Choosing the size of 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.

Centering the image at the midpoint between the two tracks' current points

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).

Export to QuickTime

Click to view exported movie

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