A calculator for polynomialsCS109 Programming ProjectsBloxorzPlotting a trajectory on a map

Plotting a trajectory on a map

When I go hiking or travelling, I often keep my GPS device running to generate a log. The device saves this log in a GPX-file, which looks like this:
<?xml version="1.0"?>
<gpx version="1.0" creator="maemo-mapper" xmlns="http://www.topografix.com/GPX/1/0">
      <trkpt lat="36.315860" lon="127.419857">
      <trkpt lat="36.315849" lon="127.419865">
      <trkpt lat="36.315833" lon="127.419867">
      // ... many lines omitted ...
      <trkpt lat="36.285810" lon="127.388883">
As you can see, the file records the location in latitude and longitude in degrees, as well as the elevation in meters, and the time the location was measured.

Here is zip-file gpx.zip with a few sample GPX files.

In this project we will write a program to analyze such a GPS log in GPX format.

Location class

We start by defining a class Location. It stores latitude, longitude, and elevation of a point on Earth. It also implements the minus operator - to compute the distance between two locations in meters:

data class Location(val lat: Double, val lon: Double, val ele: Double) {
  operator fun minus(rhs: Location): Double = 
    haversine(lat, lon, rhs.lat, rhs.lon)

To compute the distance, we need to use the Haversine formula. To save you time, I have already implemented this function:

fun haversine(lat1: Double, lon1: Double,
              lat2: Double, lon2: Double): Double {
  val EARTH_RADIUS_IN_METERS = 6372797.560856
  val rlat1 = Math.toRadians(lat1)
  val rlat2 = Math.toRadians(lat2)
  val sdlat = Math.sin(Math.toRadians(lat2-lat1)/2)
  val sdlon = Math.sin(Math.toRadians(lon2-lon1)/2)
  val a = sdlat * sdlat + Math.cos(rlat1) * Math.cos(rlat2) * sdlon * sdlon 
  val c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1.0 - a))

You can use the Location class like this:

>>> val p1 = Location(36.315860, 127.419857, 95.0)
>>> val p2 = Location(36.296373, 127.417900, 397.0)
>>> p1
Location(lat=36.31586, lon=127.419857, ele=95.0)
>>> p2
Location(lat=36.296373, lon=127.4179, ele=397.0)
>>> p1 - p2
Note how we can use subtraction for Location objects to compute the distance.

Reading the GPX file

We now need a function readGpx(fname: String): List<Location> that reads a GPX-file with the given file name, and returns the locations in a list.

A GPX file is in XML-format, so the "correct" way of reading it would be to use an XML-parser. That goes far beyond the material of CS109, so we will simply use some string matching to get the job done:

Simple analysis

Write functions to compute the total length of the track in your GPX file, the minimum elevation of the track, and the maximum elevation. Here is a tip for getting the minimum and maximum elevation:

>>> val gpx = readGpx("GPX/Bomunsan1.gpx")
>>> val heights = gpx.map { it.ele }
>>> val minHeight = heights.min()
>>> val maxHeight = heights.max()
>>> minHeight
>>> maxHeight

Creating an elevation profile

In the next step, we want to create a plot of the elevation as a function of the distance traveled along the track, like this:

Plot of elevation profile of Bomunsan1.gpx

You can learn about drawing an image and saving it to a file in the tutorial.

I suggest a size of about 800x600 pixels for the plot (or make it a command line parameter). Plot the minimum and maximum elevation near the left corners of the plot.

If you have time left (but not enough to do the next task), you can also plot vertical lines at 1000m intervals, or mark the extreme points by plotting the height there.

Plotting the track

The final task is a bit more advanced. Only start this if you have enough time. Our goal is to actually plot the track on a map, like this:

Track of Bomunsan1.gpx on openstreetmap map

We will use maps from openstreetmap.org. (We could also use Naver or Google maps, but those require you to register as a developer and get an "API key" to download maps. The principle would be the same.)

Openstreetmap and Google maps use the Mercator projection (Naver only covers Korea, so it can use a different projection that is more accurate for a smaller region, a transverse Mercator projection.)

I've written two functions that convert between latitute/longitude and Mercator coordinates:

fun mercator(lat: Double, lon: Double): Pair<Double, Double> {
  val b = Math.toRadians(lat)
  return Pair(Math.toRadians(lon), Math.log(Math.tan(b) + 1/Math.cos(b)))

fun latlong(x: Double, y: Double): Pair<Double, Double> {
  return Pair(Math.toDegrees(Math.atan(Math.sinh(y))), Math.toDegrees(x))
Here, longitude is in degrees between -180 (west) and +180 (east). Latitude is in degrees between -85 (south) and +85 (north).

The Mercator coordinates \(x\) and \(y\) both range from \(-\pi\) to \(+\pi\), where \((0,0)\) is again the point of latitude and longitude zero.

There are 19 zoom levels, numbered 0 to 18. At zoom level 0, the whole earth fits on a 256x256 pixel square. The point with coordinates \((-\pi,-\pi)\) is the bottom left corner, the point with coordinates \((\pi,\pi)\) the top right corner. At zoom level 1, the whole earth fits on a 512x512 pixel square. Here is the earth at zoom level 1:

Quarter earth Quarter earth
Quarter earth Quarter earth
Increasing the zoom level by one doubles the resolution, so at the final zoom level 18 the whole earth would need \(2^{26}\times 2^{26}\) pixels.

Fortunately, Openstreetmap provides the functionality to download a map of exactly the region you want of the size you want. The following function returns the URL where you can get this map:

fun mapurl(lat: Double, lon: Double, zoom: Int, 
	   width: Int, height: Int): String = 
  "http://staticmap.openstreetmap.de/staticmap.php?center=" +
  "%.4f,%.4f&zoom=%d&size=%dx%d".format(lat, lon, zoom, width, height) +
You have to provide latitude and longitude of the center of the desired map in degrees, the zoom level between 0 and 18, and the width and height of the desired map in pixels. Accessing this URL will return a PNG file containing the desired map.

You can download a URL and store it in a file using the following function:

fun download(url: String, fname: String) {
  System.err.print("Downloading map ... ")
  val stream = java.net.URL(url)
  val rbc = java.nio.channels.Channels.newChannel(stream.openStream())
  val fos = java.io.FileOutputStream(fname)
  fos.getChannel().transferFrom(rbc, 0, 0x1000000)
Please write your program such that it only downloads the map once. Once you have the map on your local disk, use it from there while debugging your program.

Your strategy should be the following:

Bonus task

Write your program such that it does not use the staticmap.openstreetmap.de service, but instead downloads map tiles directly from tile.openstreetmap.org and assembles them to a map itself. See here for details on the tile numbering.

A calculator for polynomialsCS109 Programming ProjectsBloxorzPlotting a trajectory on a map