Animated Watercolour Map for Cycle Tour/Race Videos

Introduction
I always record GPS traces of my cycle expeditions, as well as often recording videos, so when I go somewhere interesting I like to combine the two when I get back home. After coming across Stamen's beautiful watercolour maps, I knew that I had found the missing ingredient. My process now is twofold, I use a Python script to retrieve the map and plot the GPS trace to an image, then I overlay and animate the images in a video editor. I have tried to write a guide that will be easy to follow for both the novice python programmer as well as the novice video editor.

The python code is definitely not the most streamlined, because I have focused on readability and cutting the process into discrete steps, once you know how it's done it will be very easy to modify it to your own needs. Everything is done using standard modules.

If you already know how to generate a map image and an trace image, then jump ahead to step 8 (or if you just want the complete script, head to the previous step)

Licensing
All of the software that I used is free/open source. The maps from Stamen are licensed under creative commons. Please don't abuse them by trying to download the entire world at maximum zoom. Obviously you need to license you video accordingly, contact Stamen before using them commercially.

Software Used

• Endomondo (for recording GPS tracks)
• Scripting: Python
• Text Editor
• Image Editor: GIMP
• Video Editing: Kdenlive (can anybody recommend a good free windows option that supports masking?)

Example

Create the Script
First of all, make sure that you have python installed by typing "python" in a command line window. If you don't, head here. I have tested the script with version 2.7.1

Create a project folder a make a new text file in it that will be your script, I called mine "overlay.py". Open your file with a text editor like gedit (Ubuntu) or Notepad++ (Windows). Regular notepad will work too, but you won't have any syntax highlighting.

Import the Data

Our first objective is to import the GPS trace as a simple list of lat,lon pairs. There are number of formats for storing GPS traces, but I will just deal with comma-separated values (CSV)  and GPS eXchange Format (GPX) because there are plenty of tools to convert almost any format to any other one (GPSVisualiser, is my favourite online converter).

A CSV file can then be easily read in to a python list with the following simple function which makes use of Python's built in "csv" module.

def traceImportCSV(fname):

    import csv
    trace = []
    for lat,lon in csv.reader(open(fname,'r')):
        trace.append([float(lat),float(lon)])
    return trace

Most sites/devices have the option to save routes as GPX (such as Endomondo's export function), so if we can read in a GPX file directly that would be handy. The following bit of code reads the file line by line and uses regular expressions to look for text that says "lat=" or "lon=" and  retrieves whatever non-whitespace characters fall between the following quotation marks.

If both are found on one line, then that lat,lon pair is added to the list. For simplicity the only unusual condition that I have handled is having lat,lon appear as lon,lat instead. If the flags are on different lines, or the file contains items that are not track points, strange results may occur.

def traceImportGPX(fname):

    import re
    trace = []

    for line in open(fname,'r'):
        matchLat = re.search(r'.* lat=\"(\S*)\".*',line)
        matchLon = re.search(r'.* lon=\"(\S*)\".*',line)               
        if matchLat != None and matchLon != None:
            lat = matchLat.group(1)
            lon = matchLon.group(1)
            trace.append([float(lat),float(lon)])

    return trace 

As with most web based maps, the image that you see on your screen is made of of dozens of tiles which are retrieved from the map servers and arranged in the web browser. The tiles retrieved are determined by the geographic coordinates as well as the zoom level. In order for us to create a nice big map that covers the whole GPS route, we need to retrieve all of the relevant tiles at an appropriate zoom and arrange them correctly.

The tiles naming convention is described in depth here. Using that definition the following functions are defined to convert from lat,lon to x,y and visa versa.

### Function to convert lat,lon degrees to tile x/y number ###
def deg2num(lat_deg, lon_deg, zoom):

  lat_rad = math.radians(lat_deg)
  n = 2.0 ** zoom
  xtile = int((lon_deg + 180.0) / 360.0 * n)
  ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
  return (xtile, ytile)

### Function to convert xy tile to NW corner of tile ###
def num2deg(xtile, ytile, zoom):

  n = 2.0 ** zoom
  lon_deg = xtile / n * 360.0 - 180.0
  lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
  lat_deg = math.degrees(lat_rad)
  return (lat_deg, lon_deg)

In order to know which tiles to retrieve, we need to determine the boundaries of our GPS trace. Since latitude (West - East) ranges from -180 to 180 and longitude (South - North) ranges from -90 to 90, the North and East boundaries of the trace correspond to the maximum longitude and latitude respectively. The South and West borders correspond to the minimum values of latitude and longitude respectively.

The following function transposes the list of lat,lon pairs, effectively creating a list of latitudes and a list of longitudes so that it is easy to get the minimum and maximum value of each.

Once the trace boundaries are known we can use the previously defined function to determine the x,y values of the tiles on each corner.

### Function to get trace boundries ###
def traceBoundaries(trace):

    #transpose 2d list so that we have a row of lats and a row of lons
    lat = zip(*trace)[0]
    lon = zip(*trace)[1]

    return {"north":max(lat),"south":min(lat),"east":max(lon),"west":min(lon)}Since we now have a method of determining the boundaries of the GPS trace in lat,lon coordinates and a function to convert those lat,lon coordinates to x,y tile numbers it is a simple matter to determine the range of tiles to be downloaded for a given zoom level. The tiles are numbered with 0,0 in the top left (North-West) corner, so the NW boundary corresponds to the top left tile and the SE boundary corresponds to the bottom right tile.

### Determine tile range given boundaries and zoom ###
def determineTileRange(boundaries,zoom):
    Xmax,Ymin = deg2num(boundaries["north"],boundaries["east"],zoom)
    Xmin,Ymax = deg2num(boundaries["south"],boundaries["west"],zoom)
    return {"xMin":Xmin,"xMax":Xmax,"yMin":Ymin,"yMax":Ymax}

The tiles are available from Stamen's map server (or whichever other map you choose to use) with the URL http://tile.stamen.com/watercolor/ZOOM/X/Y.jpg

By Using Python's "os.path" tool we can make the script cross-compatible between Windows and Linux systems. After a simple file check to determine whether we have downloaded that particular tile in the past, we use Python's "urllib" module to retrieve the image and store it in a local directory with the same structure as the web server.

It is worth noting that the extension of the watercolour map tiles is ".jpg" which is unusual, maps are normally saved as ".png" because their limited palette and sharp edges are better compressed that way.

### Take a tile range and download them (if not locally present) ###

def getTiles(xyRange,zoom):

    #set acive directory to that of the script
    currentdir = os.curdir
    tileDir = os.path.join(currentdir,"tiles")
    #define the tile server
    tileServerUrl = "http://tile.stamen.com/watercolor/"

    #create a list of all the x and y coordinates to download
    xRange = range(xyRange["xMin"],xyRange["xMax"]+1)
    yRange = range(xyRange["yMin"],xyRange["yMax"]+1)

    for x in xRange:
        for y in yRange:
            #define the file name
            tileFileName = str(y)+".jpg"

            #define the local path as well as the complete path to the local and remote files
            localPath = os.path.join(tileDir,str(zoom),str(x))
            localFile = os.path.join(localPath,tileFileName)
            remoteFile = tileServerUrl+str(zoom)+"/"+str(x)+"/"+str(y)+".jpg"

            #check if the file exists locally
            if not os.path.isfile(localFile):               
                print "retrieving "+remoteFile
                #if local directory doesn't yet exist, create it
                if not os.path.isdir(localPath):
                    os.makedirs(localPath)
                #retrieve the file from the server and save it   
                urllib.urlretrieve(remoteFile,localFile)

Now that we have all of the map tiles saved locally it is a simple task to automate merging them all into one image. This is the point when you will be glad that you are scripting the solution instead of trying to do it by hand.

Python proves some basic image manipulation tools which we will use. The tile size is defined within the function as 256x256 which something of a standard, but worth confirming if you use tiles from a different source.

Simply loop through the required range of the tiles, reading in each tile and pasting it into the image. Since the "paste" command places the top left corner of the image at the specified coordinates, each tile will be placed at (X*tilesize,Y*tilesize) where X,Y ranges from zero to number of tiles in X or Y direction.

### Merge tiles into one image ###
def mergeTiles(xyRange,zoom,filename):
    import Image
    tileSize = 256
    tileDir = os.path.join(os.curdir,"tiles",str(zoom))

    out = Image.new( 'RGB', ((xyRange["xMax"]-xyRange["xMin"]+1) * tileSize, (xyRange["yMax"]-xyRange["yMin"]+1) * tileSize) )

    imx = 0;
    for x in range(xyRange["xMin"], xyRange["xMax"]+1):
        imy = 0
        for y in range(xyRange["yMin"], xyRange["yMax"]+1):
            tileFile = os.path.join(tileDir,str(x),str(y)+".jpg")
            tile = Image.open(tileFile)
            out.paste( tile, (imx, imy) )
            imy += tileSize
        imx += tileSize

    out.save(os.path.join(os.curdir,filename))

The final step in the GPS process is to draw the path onto an image of the same size as the map.

Drawing the path requires us to convert GPS coordinates to something that we can relate to the x,y coordinates of an image. Since the image that we generated is based on the Mercator Projection, any horizontal/vertical line in it has an associated latitude or longitude.

Subtracting the lat,lon of the Northwest corner from each coordinate effectively makes the Northwest coordinate correlate with pixel (0,0).

The values then need to be scaled so that a coordinate on the Southwest corner correlates to the (X resolution,Y resolution) pixel. To do this we determine the lat,lon of the point (relative to the Northwest corner) as a percentage of the image's lat,lon range and multiply by the resolution in that direction. The lat,lon range is simply (North border - South border, East border - West border).

### Draw Path Image ###

def drawTraceMask(trace,xResolution,yResolution,traceBoundaries,zoom,filename):
    import Image
    import ImageDraw

    # Get XY number of NW and SE corner tiles
    xy_nw = deg2num(traceBoundaries["north"],traceBoundaries["west"],zoom)
    xy_se = deg2num(traceBoundaries["south"],traceBoundaries["east"],zoom)

    # get lat,lon of corners
    # (since the function returns the NW corner of a tile,
    # we need lat,lon of X+1,Y+1 for the SE corner)
    NW = num2deg(xy_nw[0],xy_nw[1],zoom)
    SE = num2deg(xy_se[0]+1,xy_se[1]+1,zoom)

    # The image boundaries are actually different, because
    # they are the boundaries of the tiles, not the trace
    # define the new boundaries

    mapBoundaries = {}
    mapBoundaries["north"] = NW[0]
    mapBoundaries["south"] = SE[0]
    mapBoundaries["west"] = NW[1]
    mapBoundaries["east"] = SE[1]

    # Offset to ensure that NW corner is 0,0
    latOffset = -(mapBoundaries["north"])
    latDivisor = mapBoundaries["north"]-mapBoundaries["south"]
    lonOffset = -(mapBoundaries["west"])
    lonDivisor = mapBoundaries["east"]-mapBoundaries["west"]

    out = Image.new( 'RGB', (xResolution, yResolution) )
    draw = ImageDraw.Draw(out)

    firstRun = True
    for lat,lon in trace:
        # Convert zeroed lat,lon into x,y coordinates
        # this will need correction for northern hemisphere       
        x = abs(int(xResolution*((lon + lonOffset)/lonDivisor)))
        y = abs(int(yResolution*((lat + latOffset)/latDivisor)))
        if firstRun:
            firstRun = False
        else:
            draw.line((x,y,xPrev,yPrev),fill="white",width=10)
        xPrev = x
        yPrev = y
    del draw

    out.save(os.path.join(os.curdir,filename))

All that is left is to use the functions that we have defined. I have attached the completed script as well.

# define parameters
zoom = 12
trace = traceImportCSV("test_trace_4.csv")

# determine the boundaries of the trace
boundaries_trace = traceBoundaries(trace)

# determine xy numbers of boundary tiles
tileRange = determineTileRange(boundaries_trace,zoom)

# count number of tiles in x and y direction
xTiles = tileRange["xMax"]-tileRange["xMin"]
yTiles = tileRange["yMax"]-tileRange["yMin"]
numTiles = xTiles*yTiles

# download tiles if needed not locally available
getTiles(tileRange,zoom)

# merge tiles into oneimage
mergeTiles(tileRange,zoom,"output-map.jpg")

# draw the path
# Note: the range in "tilerange" refers to the NW corner, but our image extends on block further
drawTraceMask(trace,256*(xTiles+1),256*(yTiles+1),boundaries_trace,zoom,"output-mask.jpg")

Now that we have an image of the map and an image of the GPS trace of exactly the same size, we can overlay the two and make the image or video that we wanted.

The most flexible way to use the trace image is as a mask. When a black and white image is applied to a picture as a mask the black portions become transparent and the white portions opaque. Applying the black and white path image to a solid colour layer as a mask will result in a path of that colour.

The black and white mask can be modified as desired, I opened the path in GIMP, selected the white section with the magic wand and grew then shrunk the selection to fill in any gaps, after that I applied a small amount of the "oil painting" filter to blur the edges slightly and achieve a look that was more consistent with the watercolour theme.

Any fully featured video editor should be able to use the mask image in the same way that one uses it in GIMP/Photoshop. The only difference in video editing is that you we want to define an additional mask that moves over time, revealing more of the GPS trace.

This is a brief overview of how to achieve the effect in the free and open source Kdenlive, the same concepts will be applicable (but probably easier to implement) in Adobe Premier, Sony Vegas etc.

Since I am using Kdenlive, which has all sorts of bizarre idiosyncrasies, I find it easier to save a version of the image which is just the red trace on a transparent background (using .png format because it supports transparency). When using images with transparency in Kdenlive make sure that the "transparent background" flag is checked in the clip's properties.

Once I was happy with my two images I brought them into Kdenlive and placed them on two video layers, the one with the path above the one with the map.

The "rotoscoping" effect is used to make a layer mask that can be moved around with time. First add "rotoscoping" to the effect tree of the path clip, then draw a polygon in your project window (as seen in the images). When the "alpha operation" of the rotoscoping effect is set to "subtract" then anything inside of the polygon will be invisible. By changing the shape of the polygon at various keyframes throughout the clip the rotoscoping mask can be made to move and reveal more of the trace.

Since the rotoscoping effect only affects the masking layer, a "composite" transition needs to be applied between the two layers in order for the effect to be visible.

rotoscoping in Kdenlive is explained fully here

Example