Voyager

I run, bike, hike, and ski. I do many other things, as well, but these are the primary ways that I am able to explore the world around me. For me, that means finding new roads and trails, either far from home or hiding in my own backyard. To this end, I spend an unreasonable amount of time perusing maps, usually on OpenStreetMap, and creating routes with On The Go Map. However, roads are sometimes too remote or too infrequently travelled to be marked on any map, in which case I turn to satellite imagery and add them myself.

Another interest of mine is quantifying things. This has taken the form of analyzing years of my Spotify listening history, crunching race results to accurately predict the winner of the 2019 MSHSL Cross Country Championship, or creating tools to enhance split-taking for track meets. All of my favorite projects have centered around extracting data from the world and deriving meaning (though it can certainly be said that all software is just data manipulation at the end of the day).

This project is the natural intersection of these two interests, an attempt to quantify my exploration of the physical world.

I had quite a lot of inspiration for this project. Most significantly was the Pac Tom project by Tom Murphy, in which he ran the length of every street in Pittsburgh, making some interesting maps along the way. Another inspiration are Strava’s heatmaps, which display a summary of all users’ GPS traces. Both of these fall short of what I’m interested in, unfortunately. For the Pac Tom project, he manually marked every completed road using Google Earth. Not exactly something I have the patience for. The Strava heatmaps beat this by being automatically generated for paying users, but they have no options for varying how fine or coarse your GPS traces are on the map. This means no satisfying ‘completion’ of any areas, unless you are to walk back and forth across an area, slowly ‘painting in’ a region.

This also isn’t the first time I’ve tried to do something like this. A few years back, I created a little tool that allowed me to view which of St. Paul’s roads I had run on. Using publicly available descriptions of every road within the St. Paul city limits, I was able to generate GPS coordinates at an approximately 10 meter interval along each road. Comparing these points to every GPS coordinate I had recorded on a run, I marked a road segment as complete if I had run within 10 meters of it. This worked great as an automated version of the Pac Tom project, but shared most of the same pitfalls. It would only work to record activity on marked roads. No unmarked roads, no trails, and certainly no bushwacking. Additionally, I’d have to add data for every new jurisdiction that I’d want to record progress for. These were two big deal breakers for me that meant that I didn’t move any further with this project.

Requirements

With this past attempt under my belt, and knowing the shortcomings of some other similar projects, I could put together some requirements. In order to be successful, this tool needs to do two things. It must operate independently of marked roads, and it must create some sense of completion.

Operating independently of marked road

This project can only accurately describe the GPS traces I generate during my activities if it can do so without the context of known and recorded geographic features. It doesn’t happen especially often, but I don’t want there to be any loss of information if the road or trail I am on isn’t documented on any map, or if there is no road or trail at all.

Creating a sense of completion

I need to have some way of saying that I have travelled the entire length of some road, or I have ran the entirety of some region. Essentially, I just need to be able to mark off regions as being ‘done’. The root case for this can be that I have recorded at least one GPS point within some region. This works for small areas, but it hardly makes sense to say that you have ‘completed’ an entire county after running along a single road. Instead, you could define a region recursively, and say that you have completed a region when you have completed all of its component regions. In my previous attempt, this could have meant that a city was completed when I had covered all of the roads. This idea adds another requirement to the project: I need a scheme for dividing the entire planet into recursively defined regions. Also, it would be good if I can have regions of mostly similar size and shape, so that they can have a consistent ‘weight’ in my mind.

Dividing a sphere

Without using features like roads or jurisdiction bounds, I need to be able to recursively divide the world into regions which have a similar area on each level within the tree. Another way of formulating this problem is that I need to be able to divide up the surface of a sphere (the Earth) into regions. This is commonly done using spherical polyhedra, polygons that approximate a sphere as you add more and more faces. These are great for this application, as they can be constructed with really any number of faces, with 4 at the very minimum, and no real upper bound on their complexity. Class I spherical polyhedra would be the way to go for me, as they mostly consist of identically-sized equilateral triangles.

While researching my options, I stumbled upon a paper called Indexing the Sphere with the Hierarchical Triangular Mesh, which perfectly describes the problem I am trying to solve. The paper describes the Hierarchical Triangle Mesh (or HTM), a data structure that allows spherical data to be mapped to a recursively-defined spherical polyhedra. Starting with an octahedron, the faces are divided into 4 self-similar triangles (called trixels) up to some target depth as data points are inserted. The HTM is used for applications in astronomy, where stars and other objects can be described with their position on the celestial sphere, or earth sciences, where mapping complex datasets with geographic context is common. This is the route that I ended up taking.

Implementation

I use a Garmin Forerunner 245 to record GPS data for all of my activities. This data gets uploaded to the Garmin Connect service, and then on to Strava. Handily, Strava provides an option to download a full archive of all of your data (found under Settings, My Account, Download or Delete Your Account). This includes the raw files from every uploaded activity. In my case this was a pile of some GPX files, but mostly FIT files. I wrote a script in Golang to read all 2,000 activities, and dump the ~5 million GPS points into a PostgreSQL database. I also added a record of each activity, including the date and type of the activity. With the points, I added an activity ID and the order of the point within it’s activity, in case I want to view per-activity data in the future.

With all the data readily available, I could start working on constructing an HTM from the points. To do this, I needed to be able to determine if a point lies within a certain trixel. While the paper I referenced earlier described methods for querying what trixels lie within some arbitrary bounds, it didn’t describe with much specificity any way to insert data into the data structure (or if it did, it eluded my understanding). At first glance, it seems like it’s the same as determining if a point lies within a triangle. In reality, what you are really doing is creating a conical region, defined by the Earth’s center, and rays cast through the trixel’s vertices from the Earth’s center. I’ll spare both you and myself from more of the math, but I ended up with the following function:

func (t *Trixel) ContainsPoint(point []float64) bool {
  if dot(t.CrossAB, point)*t.D < -Tolerance {
    return false
  }
  
  ct2 := cross(t.A, point)
  if dot(ct2, t.C)*t.D < -Tolerance {
    return false
  }
  
  ct1 := cross(point, t.B)
  return dot(ct1, t.C)*t.D >= -Tolerance
}

t.CrossAB is the cross product of vertices A and B, t.D is the dot product of CrossAB and vertex C, point is in cartesian form. CrossAB and D are calculated when a trixel is created, as those values are not dependent on the point which is being checked.

With this finally working reliably, I could write a function to generate an HTM from an array of cartesian coordinates. When creating an HTM, you start from a set of 8 trixels, which represent the faces of the octahedron that is the initial state of an HTM. To add a point to the HTM, I check it against all of these 8 trixels. When I find the trixel that contains the point, I generate the children for that trixel, and check against those trixels. This is repeated to some defined maximum depth. In my case, I also record the number of points that each trixel contains, which I increment as I traverse the tree.

Initially, to get the children of a trixel, I attempted to read them out of the database, creating them if they weren’t present. I would then increment the point count, and update or insert the trixels when I had finished adding a batch of points. I did this based on the assumption that it would be faster to read the trixels from the database than it would be to calculate their vertices, but it turned out to be the complete opposite. This makes sense when you consider that a trixel is divided into its 4 children by finding the midpoint of each side, so a total of 9 additions and 9 divisions. This significantly shortened the amount of time required to generate the trixels, and simplified the entire process.

Something about the HTM is that there is no inherent stopping point for creating child trixels. Some papers suggest using a variable limit to the depth of the tree, based on factors like the density of data points in a certain region. For simplicity, I chose to have a fixed limit of 17. In theory, this would allow a maximum of 34 billion trixels, each covering approximately 14 thousand square meters. As a square region, it would be approximately 120 meters to a side. I considered pushing the limit to 18, but 17 gets the trixels small enough that I’d have to run down most streets of a normal city block in order to reach them all.

Accuracy Issues

With dividing the Earth into so many regions, I was having difficulty accurately determining which trixels contained which points. Early on, I was displaying all the GPS points alongside the trixels that were generated from them, and I could plainly see that there were trixels that contained no points, and points that were not contained by any trixel. This inaccuracy also manifested when displaying trixels of varying depths. Visible gaps or overlaps appeared between adjacent trixels, showing that this was even affecting my ability to accurately display the trixels in their correct position. This is caused by the fact that I initialized the HTM using what you might call a ‘unit octahedron’, where the vertices had coordinates like (1, 0, 0), (0, 1, 0), (-1, 0, 0), and so on. The issue with this arises when you begin dividing the faces to produce higher depths in the HTM. At each depth, the coordinates for the vertices would be a number in the form

\[\dfrac{n}{2^{d-1}}\]

where \(n\) is an integer in the range 0 to \(2^{d-1}\) inclusive, and \(d\) is the depth, ranging from 1 to 17 in my case.

I first tried to correct this by using Go’s math/big library, which provides types and methods for high-precision math. While I was able to resolve the problem while keeping all the data in memory, whenever I read or wrote to the database, I was losing precision again. This was probably user error, but along the way I came up with a much more elegant solution. Looking back at the equation for the coordinates, we have an integer denominator and an integer numerator. For us, that means we can multiply by the denominator, and we’ll be left with an integer! And what’s even better is that for different levels, we know that multiplying by the higher level’s denominator will also leave us with an integer value. Looking at the simplified formulas for values of coordinates for levels 4, 5, and 6 we have:

\[\dfrac{n_1}{2^3}, \dfrac{n_2}{2^4}, \dfrac{n_3}{2^5}\]

Multiplying all by \(2^5\) yields

\[2^2n_1, 2n_2, n_3\]

Which are all integers! Using this, we can ensure that all the coordinates for our HTM are integer values, and thus stored with perfect precision (assuming our values stay within the 64-bit integer range).

With this new scheme, I initialized the HTM with the same coordinates as before, just multiplied by a factor of \(2^{16}\), given my chosen maximum depth of 17. Now, I can store the coordinates as int64, and update the database to reflect this.

Displaying the results

While everything up to this point was all fun and I learned about some interesting topics, I didn’t exactly do all of this to end up with a pile of numbers. I wanted to be able to share this with others, so making it available through a website was a must. Going off of that, I found CesiumJS, a seemingly very robust JavaScript library for GIS applications. It made it pretty painless to render all the trixels, which I could group together into ‘primitives’ to accelerate rendering. I could then simply toggle between levels once it had been rendered for the first time. I added buttons to allow the user to select a specific depth of trixel to display. The keen eyed might notice that these buttons start at 4. I encountered a strange error with Cesium where if I attempted to render a polygon that had a vertex at one of the poles, it crashed the renderer. I probably could have fixed this by checking if a coordinate is at a pole, and offsetting it a small amount, but I decided to just remove the option to view those levels. This will probably come back to bite me the next time I try to display the data from my most recent arctic expedition, but I’ll burn that bridge when I get to it.

Something I found while working on this step was that the server was spending a lot of time reading from the database and formatting the trixels whenever the client requested a new level. For the highest level, it was taking multiple seconds to perform this operation. There were a few steps of this that could be optimized. Firstly, I had stored the trixels using cartesian coordinates. This was important, as that was the form required for the equations used to check if a trixel contains a point, and that was the form generated when creating the HTM. The downside to this is that I needed latitude/longitude pairs to display the trixels with Cesium. This is simple to implement, but requires two ATan2 operations and a square root, all expensive operations. I resolved this by also storing the latitude/longitude pair for each coordinate when a trixel was created. I further optimized these requests by returning only the latitude and longitude to the client in an array, whereas before I had simply been dumping an array of trixels into JSON.

Again I could optimize these requests by reducing the precision for the coordinates I was returning to the client. The default behavior for when converting a float to JSON using Go’s json.Marshall function is to retain the least number of decimal digits required to reproduce the value exactly should it be converted back to a float. This is certainly desirable behavior, but by reducing the precision I could significantly reduce the size of the response to the client with no perceivable loss of accuracy. This might sound like I am back-pedalling on my earlier efforts to increase the accuracy of the trixel coordinates, but in both cases I have an abundance of precision, which I can now sacrifice a little of given that I have achieved excellent accuracy. If you didn’t quite follow that last sentence, it’s probably because I didn’t do a great job explaining it, or you need to touch up on accuracy vs precision.

You can get a feel for what sort of precision I might want to keep by thinking about the distance between two latitude/longitude coordinates. For example, the distance between 0N0W and 0N1W is 60 nautical miles, or 111.12 kilometers. With that in mind, we can see that a step of \(10^{-7}\) degrees westward on the equator is equal to 0.000011112 kilometers, or approximately 1.11 centimeters. This means that if we retain 7 decimal digits, the coordinates will be at most 7.8 millimeters from the real location. I, for one, consider this to be more than enough accuracy, given that our smallest trixels are measured in thousands of square meters.

This optimization results in a meaningless performance improvement, but it was a fun experiment, nonetheless. However, a very real improvement can be found in just caching the response as a byte array for each depth of the HTM that can be requested by the client. Unfortunately, I did this in a lazy manner that requires me to restart the server whenever I upload new activities. Oh well, maybe I’ll fix this at some point.

Post-Mortem

I came into this project with a few goals that I went over earlier. I wanted to create a tool that I could use to map out my activities, regardless of roads or trails, and to create some sense of completion as I explored new areas. I believe that I was absolutely able to accomplish those stated goals, and some others that appeared along the way. I was able to convert over 5 million GPS points from 2 thousand activities into a map that I can check from anywhere, and easily view most levels of the generated HTM. Since completing this project, I find myself digging through the map, noticing places that I had been years ago but have since forgotten. There are a few things that I’d like to add to this at some point. I’d like to be able to see what trixels I can visit that would have a cascading effect, completing their parent trixels, and maybe even their parent’s parent. I also need to get around to adding Strava integration, so that activities are automatically added as they are synced to my account.

Overall, I’m really happy with how this turned out, and I think this will be a tool I’ll be playing with for some time to come. You can find the completed project here, and the source is available on my GitLab.