A Soul Trip across Northern Europe. GPX file generated by namereducegpx.py and visualised with Qgis
While working on Soul Walker, I’ve been revisiting the files we made for the Nottingham University research phase last year and wondering about those files generated by NetLogo of a ‘soul’ patiently trudging across the territory I flew over, oblivious to puddles, hedgerow, supermarkets.
The NetLogo script took a GPX file as input and using a flocking algorithm (a flock of two), generated a second GPX file derived from following the points in the first file at 6kph.
The output GPX file from such an operation is very featureless and straight. As you see from the illustration above, the ‘soul’ had to turn around in the North Sea as I flew back to Berlin before it had reached me.
I started to think about a way of extracting the places this ‘soul’ might travel through and get an idea of this fly-over territory and thought I would generate a placename poem from points along its path.
Hence the script below. It takes a GPX file as input and reduces the number of trackpoints to those that are a definable threshold apart (default is 20km) in order to ensure the list contains unique placenames and spare the reverse geocoding server.
The output of namereducegpx.py is a second GPX file with a suffix before the file extension (‘_named’ by default) which you can visualise in Qgis as I have done above. You can also use the script to simply reduce the trackpoints to the distance threshold by running with the ‘-o’ or ‘–offline’ option in which case the points are not sent to a geocoding service. There is also the option of using the Google reverse geocode api instead of the default Mapquest Open reverse geocoder (in Beta at time of writing).
Now to actually write a script that generates the poem from the GPX file.
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
namereducegpx.py
Tue 11 Dec 2012 10:24:43 CET
Loads a gpx file of a journey (preferably over many kilometres) and
reduces the track points so that only ones of a certain threshold
apart remain. Then the points are iterated through and a place name is
derived from two possible geocoding services and copied to the track
point name sub element. The final gpx file is saved with a suffix.
Argument:
path to a gpx file
Options:
-d --distance Set a threshold distance in km (default 20)
n-s --suffix Specify a suffix to be appended to the file name before
the extension (default '_named')
-o --offline Do not use Internet geocoding services. Has the effect of
reducing a gpx track only
-g --google Use Google reverse geocoder instead of default Mapquest
open api
Copyright 2012 Daniel Belasco Rogers <https://planbperformance.net/dan>
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or (at
your option) any later version.
This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
<https://www.gnu.org/licenses/>.
"""
try:
from lxml import etree
except ImportError:
print """
*************************************************
You do not have a module that this script needs
Please install python-lxml from the repositories
*************************************************
"""
sys.exit(2)
try:
from geopy import distance
except ImportError:
print 'geopy not installed - download from https://code.google.com/p/geopy/'
sys.exit(2)
import sys
import simplejson
from urllib2 import urlopen
from time import sleep
from os import path
from optparse import OptionParser
def parseargs():
"""
Get the path to the gpx file
"""
usage = """
%prog path/to/gpx/file.gpx
"""
parser = OptionParser(usage, version="%prog 0.1")
parser.add_option("-d",
"--distance",
dest="gapdist",
action="store",
type="int",
default=20,
help="Distance between points to get placenames from in km (default is 20)")
parser.add_option("-s",
"--suffix",
dest="suffix",
action="store",
type="string",
default="_named",
help="Suffix to append to output gpx file before the extension. Default is '_named'")
parser.add_option("-o",
"--offline",
dest="offline",
action="store_true",
default=False,
help="Work in Offline mode (default False). If chosen, does not look up placenames.")
parser.add_option("-g",
"--google",
dest="google",
action="store_true",
default=False,
help="Use Google api as the reverse geocoder. If not chosen, defaults to Mapquest open api")
options, args = parser.parse_args()
# warn the user they haven't passed an argument
if len(args) != 1:
parser.error("""
Please enter a path to a gpx file
e.g. python namereducegpx.py ~/dansdocs/soul_walker/mrl/gpsForAnimation/FinalSoulGPXforDB/2011-01-12-15ApartBrumNottsDan.gpx
""")
# args is a list, even if it contains only 1 item
return args[0], options.gapdist, options.suffix, options.offline, options.google
def findelement(xmlns, trkpoint, element):
"""
get text for an element.
return 'NULL' if not found
"""
try:
element = trkpoint.find(('{%s}%s' % (xmlns, element))).text.strip()
except AttributeError:
element = 'NULL'
return element
def getpoints(gpxpath):
"""
read the gpx file and return a list of points in a dictionary- try
and get everything the file has to offer i.e. speed and course if
possible
"""
pointlist = []
# parse gpxfile
tree = etree.parse(gpxpath)
root = tree.getroot()
# get namespace and define variables for the loop
xmlns = root.nsmap[None]
trkseg = '{%s}trkseg' % xmlns
for trackseg in root.iter(trkseg):
# filter out empty tracks by searching for first time element
try:
trktime = trackseg.find(('{%s}trkpt/{%s}time' % (xmlns, xmlns))).text
except AttributeError:
# time element isn't there, which means its an empty track
continue
# create a dictionary for each line in list (each track point)
for trkpoint in trackseg:
pointlist.append({
'lat': trkpoint.get('lat'),
'lon': trkpoint.get('lon'),
'ele': findelement(xmlns, trkpoint, 'ele'),
'time': findelement(xmlns, trkpoint, 'time'),
'course': findelement(xmlns, trkpoint, 'course'),
'speed': findelement(xmlns, trkpoint, 'speed'),
'name': findelement(xmlns, trkpoint, 'name'),
'fix': findelement(xmlns, trkpoint, 'fix'),
'hdop': findelement(xmlns, trkpoint, 'hdop')
})
return pointlist
def reducepoints(pointlist, gapdist):
"""
take the pointlist, calculate distances and reject points that are
nearer than defined threshold apart (gapdist)
"""
# set up the previous trackpoint (1st in list) for comparison
lat1 = pointlist[0]['lat']
lon1 = pointlist[0]['lon']
# start the output list with the first trackpoint
reducedlist = [{
'lat': lat1,
'lon': lon1,
'ele': pointlist[0]['ele'],
'time': pointlist[0]['time'],
'course': pointlist[0]['course'],
'speed': pointlist[0]['speed'],
'name': pointlist[0]['name'],
'fix': pointlist[0]['fix'],
'hdop': pointlist[0]['hdop']
}]
# start iterating through the list from 2nd track point
for x in range(1, len(pointlist) - 1):
lat2 = pointlist[x]['lat']
lon2 = pointlist[x]['lon']
gap = (distance.distance((lat1, lon1), (lat2, lon2)).m) / 1000
if gap >= gapdist:
reducedlist.append({
'lat': lat2,
'lon': lon2,
'ele': pointlist[x]['ele'],
'time': pointlist[x]['time'],
'course': pointlist[x]['course'],
'speed': pointlist[x]['speed'],
'name': pointlist[x]['name'],
'fix': pointlist[x]['fix'],
'hdop': pointlist[x]['hdop']
})
# this trackpoint will be the 'previous' trackpoint in the
# next iteration
lat1 = lat2
lon1 = lon2
return reducedlist
def gpxfrompoints(reducedlist, offline, google):
"""
prepare a gpx xml tree from the reduced list of point
"""
# prepare the gpx xml tree we're going to write
newgpxfile = etree.Element('gpx', attrib={"creator": "reducegpx.py",
"version": "1.0",
"xmlns": "https://www.topografix.com/GPX/1/0"
})
trk = etree.SubElement(newgpxfile, 'trk')
trkname = etree.SubElement(trk, 'name')
if google:
trkname.text = 'namereducegpx.py using Google reverse geocode api'
elif not google:
trkname.text = 'namereducegpx.py using Mapquest open reverse geocode api'
trkseg = etree.SubElement(trk, 'trkseg')
# iterate through the points, writing the attributes from each
# point dictionary in the reduced list
for pointdict in reducedlist:
trkpoint = etree.SubElement(trkseg, 'trkpt')
for key in pointdict:
if key in ('lat', 'lon'):
trkpoint.set(key, pointdict[key])
else:
if pointdict[key] != 'NULL':
element = etree.SubElement(trkpoint, key)
if key == 'name' and not offline:
element.text = lookuppoint(google, (pointdict['lat'], pointdict['lon']))
else:
element.text = pointdict[key]
if not 'name' in pointdict and not offline:
element = ettree.SubElement(trkpoint, 'name')
element.text = lookuppoint(google, (pointdict['lat'], pointdict['lon']))
else:
continue
return newgpxfile
def lookuppoint(google, query):
"""
using open.mapquestapi (not sure about its accuracy...)
or
https://maps.googleapis.com/maps/api/geocode/json?latlng=40.714224,-73.961452&sensor=false
"""
lat = float(query[0])
lon = float(query[1])
if not google:
url = "https://open.mapquestapi.com/nominatim/v1/reverse?format=json&lat=%f&lon=%f" % (lat, lon)
elif google:
url = "https://maps.googleapis.com/maps/api/geocode/json?latlng=%f,%f&sensor=false" % (lat, lon)
else:
print "service %s not implemented" % service
print "Looking up lat %f lon %f" % (lat, lon)
result = urlopen(url)
jsonstring = result.read()
json = simplejson.loads(jsonstring)
if not google:
place = json.get('display_name', [])
elif google:
place = json.get('results')
if place == []:
place = 'Unknown'
return place
place = place[0]['formatted_address']
print place
return place
def makenewgpxpath(gpxpath, suffix):
"""
derive a new gpxpath from current one. If a destination path is
specified, use this in constructing newgpxpath
"""
head, tail = path.split(gpxpath)
gpxpath, extension = path.splitext(tail)
newgpxpath = '%s%s%s' % (gpxpath, suffix, extension)
newfilepath = path.join(head, newgpxpath)
return newfilepath
def writegpxfile(newfilepath, newgpxfile):
"""
writes the gpx file to the new file path
"""
with open(newfilepath, 'w') as writefile:
writefile.write(etree.tostring(newgpxfile,
encoding="utf-8",
pretty_print=True,
xml_declaration=True))
return
def main():
"""
"""
gpxpath, gapdist, suffix, offline, google = parseargs()
# read points from gpx file into a list
print "Reading gpx file..."
pointlist = getpoints(gpxpath)
print "Found %d points" % len(pointlist)
# reduce the list by the defined threshold
print "Finding points that are at least %dkm apart" % gapdist
reducedlist = reducepoints(pointlist, gapdist)
print "Reduced points to %d" % len(reducedlist)
# produce a new gpx file from the reduced points list
print "Preparing new gpx file"
if not offline:
if not google:
print "Using mapquest open reverse geocode api"
elif google:
print "Using google reverse geocode api"
newgpxfile = gpxfrompoints(reducedlist, offline, google)
# generate a new file name from old
newfilepath = makenewgpxpath(gpxpath, suffix)
# print the tree to a new gpx file in the same directory and
# finish up
writegpxfile(newfilepath, newgpxfile)
print "Wrote", newfilepath
print "Script ends here sucessfully"
if __name__ == '__main__':
sys.exit(main())