GIS Quest - Part I

For me, in computer science, there exists silos of knowledge, opaque and appearing like magic when looked upon from inexperienced eyes. Eventually the intrigue becomes too much and I find an excuse to start looking under the hood. GIS is such a silo, and this is my excuse.

Over the years I’ve done a lot of running and hiking, and via that I’ve collected a lot of data. What I’ve collected isn’t homogenous in any way and as I’m not a fan of uploading personal data to third parties, especially not GPS data, I’ve done little visualization or analysis of the aggregate data.

Self-hosting tiles

A good place to start is drawing a self-hosted map and plotting a route. Thankfully, Protomaps has made this easy, it’s an awesome open-source project that allows self-hosting of map tiles using a single static file.

In my case, since I live in Bergen, Norway, the tiles I need can be retrieved with:

$ pmtiles extract https://build.protomaps.com/20231120.pmtiles bergen.pmtiles --bbox=5.282364,60.375606,5.408535,60.417412

The protomaps website provides instructions on how to visualize the downloaded tiles using popular APIs such as leaftlet or maplibregl.

Data processing

Normalization

The location data comes from various sources, in different formats such as gpx, GeoJSON, and tcx. To normalize I first converted everything to geojson using togeojson.

Filtering

To start off the dataset will be limited to an area on a nearby mountain massif called Byfjellsmassiven, around Fløyen. To constrain the data to this region I wrote a small tool to draw a polygon around the area to include:

Left click to add a point, right click to remove.

and using a slightly modified version of James Halliday’s point in polygon algorithm, the set is filtered using the following criteria:

• A trip does not start until a point is found within the area of interest
• A trip ends when there are no further points within the area of interest

const pointInPolygon = (point, polygon) => {
    const [x, y] = point
  
    let inside = false;
  
    for (let i = 0, j = polygon.length - 1; i < polygon.length; j = i++) {
        const [xi, yi] = polygon[i];
        const [xj, yj] = polygon[j];
    
        const intersect =
              ((yi > y) != (yj > y)) && (x < (xj - xi) * (y - yi) / (yj - yi) + xi);

        inside = !!(intersect ^ inside);
    }
  
    return inside;
}

Decimating

Most tracking software grabs coordinates periodically, often resulting in a lot more points than what’s needed for a fairly accurate representation. We can save a lot of space and rendering speed by reducing the number of points. For example, given the coordinates:

[[5.36482,60.38294],[5.36485,60.38295],[5.36484,60.38292],[5.364855,60.382925],[5.36487,60.38293],[5.36488,60.38293]]

The array of points can be decimated by calculating the angle between three points, and points with angles exceeding a threshold angle removed.

const decimate = (points, threshold) => {
    let decimated = [...points];
    let i = 0;
    while (i < decimated.length-2) {
        const [a, b, c] = decimated.slice(i, i+3);
        const v = [a[0] - b[0], a[1] - b[1]];
        const w = [c[0] - b[0], c[1] - b[1]];
        const angle = Math.atan2(w[1] * v[0] - w[0] * v[1],    w[0] * v[0] + w[1] * v[1]);
        if (Math.abs(180/Math.PI * angle) >= threshold) {
            decimated.splice(i+1, 1);
        } else {
            i+=1;
        }
    }
    return decimated;
}

Vanity Poster

Putting together our map and processed data, as well as some metrics collected from the data, we can generate our first visualization, the vanity poster: