DEV Community

Cover image for Advent of Code 2020: Day 12 using QGIS and Python
Yuan Gao
Yuan Gao

Posted on • Edited on

Advent of Code 2020: Day 12 using QGIS and Python

Today's challenge involves a navigating ferry. Since it's dealing with coordinates, we can do some fun visualization in QGIS

Things mentioned in this post: geographic information systems (GIS), transformation matrix, affine transformation

The Challenge Part 1

Link to challenge on Advent of Code 2020 website

The challenge is another one in which you read in instructions, and execute them. In this case, you are a ferry, navigating from instructions. Each instruction contains a command, and a number.

Four of the commands move the ferry in compass directions: N, E, S, W; two of the commands alter the heading of the ferry: L, R: and one command tells the ferry to move in its heading direction rather than a compass direction: F.

This means our ferry has internal state: its position in space, as well as its heading.

QGIS

For no reason other than that I want this blog series to be about exploring different applications of Python, I'm going to use QGIS, an open-source GIS (geographic information system) program, that can be scripted with Python.

GIS is a whole industry sector, skillset, and set of tools, and python is not a bad choice thanks to its image and data processing libraries. Common libraries in python used for GIS work include the GDAL bindings, and shapely. However, today, we're using none of those, and instead using the embeded Python inside QGIS, which comes with its own SDK/API.

Screenshot 2020-12-18 101000

I've used QGIS a lot for my work with satellite and drone imagery, but it's mostly been for loading and inspecting data, rather than scripting. So it's interesting to explore the python scripting capabilities of a program I've used a lot but haven't delved too deep into.

The code

The actual Python code is pretty straightforward: import the data, loop over it, and update the state and move the ferry, to track where it goes.

To do this, I'm going to write a function that reads the command and the heading of the ship (the only intrinsic state of the ship), and output a movement vector and the new heading, which an outer loop will use to keep track of the extrinsic state of the ship.

There isn't a strong reason for splitting the state into an intrinsic and extrinsic state, other than trying to reduce how much state has to pass through this function.

from math import cos, sin, radians

vectors = {
    "N": (0, 1, 0, 0),
    "E": (1, 0, 0, 0),
    "S": (0, -1, 0, 0),
    "W": (-1, 0, 0, 0),
    "L": (0, 0, -1, 0),
    "R": (0, 0, 1, 0),
    "F": (0, 0, 0, 1)
}

def process_command(command, heading):
    code = command[0]
    amount = int(command[1:])

    easting, northing, rotate, forward = vectors[code]
    north_delta = amount*(northing + forward*cos(radians(heading)))
    east_delta = amount*(easting + forward*sin(radians(heading)))
    heading += rotate * amount

    return east_delta, north_delta, heading 
Enter fullscreen mode Exit fullscreen mode

The slight trick here is to avoid doing a massive if/elif loop and handling each case separately. Instead, we just use a linear-algebra-like method to reduce the branches in our code, and apply the rules as a transformation (I actually blame yesterday's solution on pushing my thinking in this direction).

That's really the meat of the code. What's more interesting however, is having QGIS visualize it.

Visualization

To do that, here's the full code that includes inserting these data points into QGIS features:

from math import cos, sin, radians
from PyQt5.QtWidgets import QFileDialog
from qgis.PyQt.QtCore import QVariant
from qgis.core import (
    Qgis,
    QgsField,
    QgsProject,
    QgsVectorLayer,
    QgsFeature,
    QgsGeometry,
    QgsPoint,
    QgsPointXY,
    QgsCoordinateReferenceSystem,
    QgsMessageLog
)

vectors = {
    "N": (0, 1, 0, 0),
    "E": (1, 0, 0, 0),
    "S": (0, -1, 0, 0),
    "W": (-1, 0, 0, 0),
    "L": (0, 0, -1, 0),
    "R": (0, 0, 1, 0),
    "F": (0, 0, 0, 1)
}

def process_command(command, heading):
    code = command[0]
    amount = int(command[1:])

    easting, northing, rotate, forward = vectors[code]
    north_delta = amount*(northing + forward*cos(radians(heading)))
    east_delta = amount*(easting + forward*sin(radians(heading)))
    heading += rotate * amount

    return east_delta, north_delta, heading


def main():
    data_file = QFileDialog.getOpenFileName(None, "Open data file")
    if not data_file:
        return

    QgsMessageLog.logMessage(f"Opening file {data_file}", "AoC", Qgis.Info)

    # create vector layers
    vl_points = QgsVectorLayer("Point", "AoC Points",  "memory")
    pr_points = vl_points.dataProvider()
    pr_points.addAttributes([
        QgsField("idx", QVariant.Int),
        QgsField("cmd", QVariant.String),
        QgsField("heading", QVariant.Double),
    ])
    vl_points.updateFields()

    vl_lines = QgsVectorLayer("LineString", "AoC Lines",  "memory")
    pr_lines = vl_lines.dataProvider()

    # read data
    data = open(data_file[0]).readlines()
    QgsMessageLog.logMessage(f"{data}", "AoC", Qgis.Info)

    heading = 90 # ship starts facign east
    northing = 0
    easting = 0
    line = [QgsPoint(0, 0)]
    for idx, entry in enumerate(data):
        east_delta, north_delta, heading = process_command(entry, heading)
        northing += north_delta
        easting += east_delta
        QgsMessageLog.logMessage(f"cmd[{idx}]:{entry} E{easting:.0f}, N{northing:.0f}, H{heading}", "AoC", Qgis.Info)

        # add feature
        line.append(QgsPoint(easting, northing))

        feat = QgsFeature()
        feat.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(easting, northing)))
        feat.setAttributes([idx, entry, heading])
        pr_points.addFeatures([feat])
    vl_points.updateExtents()

    # construct linestring
    feat = QgsFeature()
    feat.setGeometry(QgsGeometry.fromPolyline(line))
    pr_lines.addFeatures([feat])
    vl_lines.updateExtents()

    QgsProject.instance().addMapLayer(vl_points)
    QgsProject.instance().addMapLayer(vl_lines)

    # disable crs (we need flat world)
    crs = QgsCoordinateReferenceSystem(None)
    QgsProject.instance().setCrs(crs)

main()
Enter fullscreen mode Exit fullscreen mode

This creates a couple of QGIS vector layers in our project, and also sets the CRS to "None", since we want a flat earth. Then, rather than accept the defaults, we can adjust the symbology a little to give us some nice colours and shapes, including setting it to use the heading property set up in the code as the rotation for the symbol

Screenshot 2020-12-18 114907

And we produce this lovely diagram of the movement of the ferry:

Screenshot 2020-12-18 115355

The arrows show the heading of the ferry; the text labels are the commands, and the colours use the Viridis colour scheme, going from dark blue to yellow (viridis does better for red-green colourblindness than some other scales).

While it's possible to do the customization of the visualization in code, for the purpose of this exercise, it was just easier to configure the rest of the visualization using the UI.

The Challenge Part 2

In the second part of the challenge, things are switched up a bit. Now, instead of directly moving the ship, a virtual waypoint is moved instead, and then the ship moved towards the waypoint. The waypoint here is effectively acting as another intrinsic state that defines a movement vector.

Now, four of the commands move the waypoints in compass directions: N, E, S, W; two of the commands rotate the waypoint around its origin: L, R: and one command tells the ferry to move in the vector defined by the waypoint.

So I adjust the code, dropping the "forward" component from the vector, here are the main changes:


vectors = {
    "N": (0, 1, 0),
    "E": (1, 0, 0),
    "S": (0, -1, 0),
    "W": (-1, 0, 0),
    "L": (0, 0, -1),
    "R": (0, 0, 1),
}

def process_command(command, ship, waypoint):
    code = command[0]
    amount = int(command[1:])

    if code == "F":
        ship[0] += amount*waypoint[0]
        ship[1] += amount*waypoint[1]
    else:
        easting, northing, rotate = vectors[code]
        angle = radians(-amount*rotate)
        old_e, old_n = waypoint
        waypoint[0] = amount*easting + old_e*cos(angle) - old_n*sin(angle)
        waypoint[1] = amount*northing + old_e*sin(angle) + old_n*cos(angle)
Enter fullscreen mode Exit fullscreen mode

Here, I've updated the function to accept the whole ship and waypoint states, and split out the "F" state to update the ship location using the waypoint vector. The remaining code now acts only to update the waypoint state. The code itself effectively implements an affine transformation matrix to achieve both the rotation and translation.

The rest of the code updates the way the points are stored for visualization:

from math import cos, sin, radians, degrees, atan2
from PyQt5.QtWidgets import QFileDialog
from qgis.PyQt.QtCore import QVariant
from qgis.core import (
    Qgis,
    QgsField,
    QgsProject,
    QgsVectorLayer,
    QgsFeature,
    QgsGeometry,
    QgsPoint,
    QgsPointXY,
    QgsCoordinateReferenceSystem,
    QgsMessageLog
)

vectors = {
    "N": (0, 1, 0),
    "E": (1, 0, 0),
    "S": (0, -1, 0),
    "W": (-1, 0, 0),
    "L": (0, 0, -1),
    "R": (0, 0, 1),
}

def process_command(command, ship, waypoint):
    code = command[0]
    amount = int(command[1:])

    if code == "F":
        ship[0] += amount*waypoint[0]
        ship[1] += amount*waypoint[1]
    else:
        easting, northing, rotate = vectors[code]
        angle = radians(-amount*rotate)
        old_e, old_n = waypoint
        waypoint[0] = amount*easting + old_e*cos(angle) - old_n*sin(angle)
        waypoint[1] = amount*northing + old_e*sin(angle) + old_n*cos(angle)

def main():
    data_file = QFileDialog.getOpenFileName(None, "Open data file")
    if not data_file:
        return

    QgsMessageLog.logMessage(f"Opening file {data_file}", "AoC", Qgis.Info)

    # create vector layers
    vl_points = QgsVectorLayer("Point", "AoC Points",  "memory")
    pr_points = vl_points.dataProvider()
    pr_points.addAttributes([
        QgsField("idx", QVariant.Int),
        QgsField("cmd", QVariant.String),
        QgsField("heading", QVariant.Double),
    ])
    vl_points.updateFields()

    vl_lines = QgsVectorLayer("LineString", "AoC Lines",  "memory")
    pr_lines = vl_lines.dataProvider()

    # read data
    data = open(data_file[0]).readlines()
    QgsMessageLog.logMessage(f"{data}", "AoC", Qgis.Info)

    heading = 90 # ship starts facign east
    ship = [0, 0]
    waypoint = [10, 1]
    line = [QgsPoint(ship[0], ship[1])]

    for idx, entry in enumerate(data):
        # plot feature
        if entry[0] == "F":
            line.append(QgsPoint(ship[0], ship[1]))

            feat = QgsFeature()
            feat.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(ship[0], ship[1])))
            heading = degrees(atan2(waypoint[0], waypoint[1]))
            feat.setAttributes([idx, entry, heading])
            pr_points.addFeatures([feat])

        # add feature
        process_command(entry, ship, waypoint)
        QgsMessageLog.logMessage(f"cmd[{idx}]:{entry} E{ship[0]:.0f}, N{ship[1]:.0f}, WE{waypoint[0]:.0f}, WN{waypoint[1]:.0f}", "AoC", Qgis.Info)

    # last point
    line.append(QgsPoint(ship[0], ship[1]))

    feat = QgsFeature()
    feat.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(ship[0], ship[1])))
    heading = degrees(atan2(waypoint[0], waypoint[1]))
    feat.setAttributes([idx, entry, heading])
    pr_points.addFeatures([feat])

    vl_points.updateExtents()

    # construct linestring
    feat = QgsFeature()
    feat.setGeometry(QgsGeometry.fromPolyline(line))
    pr_lines.addFeatures([feat])
    vl_lines.updateExtents()

    QgsProject.instance().addMapLayer(vl_points)
    QgsProject.instance().addMapLayer(vl_lines)

    # disable crs (we need flat world)
    crs = QgsCoordinateReferenceSystem(None)
    QgsProject.instance().setCrs(crs)

main()
Enter fullscreen mode Exit fullscreen mode

This generates the following visualization:

Screenshot 2020-12-18 132959

The arrows now point towards the waypoint vector; and the text label is just the command number.

The answers needed for the challenges can simply be found from the final position of the ship as output in the log messages.

Onward!

Top comments (0)