Converting waypoints to text file with GPSBabel

GPSBabel is a great command-line app to convert just about any kind of GPS file to any other. You can use it to convert, say, a GPX tracklog to a KML file in one step:

gpsbabel -t -i gpx -f "2012-10-01 track.gpx" -o kml -F "2012-10-01 track.kml"
  • -t : convert tracks
  • -i gpx : input file is GPX format
  • -f “…” : input filename
  • -o kml : output filetype kml
  • -F “…” : output filename

Bash script to convert GPX waypoints to CSV

Put this in /usr/local/bin or other $PATH directory, and chmod +x to make it executable:

# gpsbabel script to convert GPX waypoints to CSV
LEN1=${#1}-3  # length of the filename
OUT=${1:0:$LEN1}csv  # convert the input filename to *.csv
echo Converting "$1" to $OUT
gpsbabel -w -i gpx -f "$1" -o unicsv,grid=4 -F "$OUT"


  • GPSBabel in a script will fail if there are spaces in the filename UNLESS you use quotes around the command argument $1 and $OUT.

uTrack data to KML: parsing GPS data in Python

uTrack lets you download your vehicle GPS data as a CSV file. However, the CSV doesn’t include nicely-formatted lat/lon columns – instead, you get a fifth row with an address and latlon pair:

Monduli,  (36.9005,-2.9013)

My first attempt at a real Python script took me into regex, reading command line arguments, and basic CSV parsing, all of which took me about twelve times longer than they should have.

If you have Python installed, this script (save as or similar) will parse the file and save a KML file in the same directory. Usage:

python VehicleHistory_AAA123_20121001.csv

### Read uTrack report file into KML
# Load csv, system and getopt libraries
import csv, sys, getopt, re

## Read the file
fname = sys.argv[1]
with open(fname, 'rU') as csvfile:
data = csv.reader(csvfile)
#Skip the 1st header row.

## Create the KML file
f = open(fname[:-3] + 'kml', 'w')
# Writing the kml file.
The following lines create the headers for a KML file, and set up a style called 'msn_placemark_circle' which will be used for all the points.

f.write("<?xml version='1.0' encoding='UTF-8'?>\n")
f.write("<kml xmlns=''>\n")
f.write(" <name>" + fname + '.kml' +"</name>\n")
f.write('<StyleMap id="msn_placemark_circle"><Pair><key>normal</key><styleUrl>#sn_placemark_circle</styleUrl></Pair><Pair><key>highlight</key><styleUrl>#sh_placemark_circle_highlight</styleUrl></Pair></StyleMap><Style id="sh_placemark_circle_highlight"><IconStyle><color>ff0000ff</color><scale>0.590909</scale><Icon><href></href></Icon></IconStyle><ListStyle></ListStyle></Style><Style id="sn_placemark_circle"><IconStyle><color>ff0000ff</color><scale>0.5</scale><Icon><href></href></Icon></IconStyle><ListStyle></ListStyle></Style>')
# Loop through all the rows of the CSV file and create a point for each
for row in data:
f.write(" <Placemark>\n")
# f.write(" <name>" + str(row[0]) + "</name>\n") # best to leave <name> out as it clutters the view in GE
f.write(" <description>" + str(row[0]) + "-" + str(row[1]) + "-" + str(row[2]) + "</description>\n")
f.write(" <styleUrl>#msn_placemark_circle</styleUrl>")
f.write(" <Point>\n")

# Regex to find lat/lon numbers:
p = re.compile('-?[0-9]+.[0-9]+')
m = p.findall(str(row[5]))
lat = m[1]
lon = m[0]

# Write that pair as coordinates for KML:
f.write(" <coordinates>" + lon + "," + lat + "," + str(0) + "</coordinates>\n")
f.write(" </Point>\n")
f.write(" </Placemark>\n")

Using FieldPyculator for QGIS

Python power for your calculations

It’s available in the QGIS plugin repository and thus installable with ‘Fetch Plugins’; source available on github. According to their webpage, FieldPyculator was developed by GIS-Lab. They have a Russian-language introduction to the plugin on their webpage, with an excellent Google translation.


Use the toolbar button or the Plugins|FieldPyculator menu to open the dialog.

  • Enter your Python code in the “Field expression” text box. Values are assigned by the value = statement – this must exist in some form in the final code.
  • If you need other libraries, or to set up variables or any other code that runs once before the field expression, tick the ‘Global expression’ button and enter the code in that text box.
    • For example, to use the datetime functions, you would enter import datetime as the Global Expression and could then use a command like value = as the Field Expression.
    • You could also set a variable in the Global Expression (such as f = 123.45‘) and use that in calculations.


  1. There doesn’t seem to be a way to create a new field in the dialog, unlike the usual field calculator. You have to add a new column somewhere else and then use it here.
  2. I can’t seem to get ‘IF’ statements working …

Building a GPX Spatialite database

Every aerial survey generates dozens of tracklogs and waypoints – and these tracks are a very important part of the data for the survey. We can use tracklog data to confirm ground speed, actual coverage of the survey, off-track measurements and other metadata about the transects flown.

I’m playing with storing all the survey data in a Spatialite database – RSO and FSO data, photographs, and the spatial files for the planning and reporting. While it’s straightforward to import a shapefile directly into Spatialite, the shapefiles I have in a QGIS project have usually been edited, and it would be good to keep the original tracks straight from the GPS for complete records.

GPX to Spatialite

Luckily, OGR (command line) can work with both GPX and Spatialite databases. A great post from ‘Scratching Surfaces’ gave me the initial code for this:

ogr2ogr -append -f "SQLite" -dsco SPATIALITE=yes
    -dsco INIT_WITH_EPSG=yes -t_srs epsg:4326
    gpx.db GPXFile.gpx tracks -nln TrackTable
ogr2ogr command line for ogr utilities; on OSX, installed with the GDAL complete package from Kyngchaos.
-append Add data to the the table
-f “SQlite” -dsco SPATIALITE=yes Use SQLite as the base format (the basis of Spatialite) with Spatialite extensions loaded.
 -dsco INIT_WITH_EPSG=yes Initialise table with CRS? Yes!
-t_srs epsg:4326 … and use 4326 (WGS84) as the CRS.
gpx.db The Spatialite database
GPXFile.gpx tracks The GPX file and the geometry to import – just tracks in this case.
-nln TrackTable ‘New layer name’ and the name of the table (TrackTable). Because of ‘append’, this will be ignored if there’s already data.

 Loading a directory

It would be nice to keep track of the original filename for each track that gets imported. I added a new column to the table after first importing:


A bash script can then iterate through all the GPX files in a directory:


DB=$"spatial_db2.sqlite"; # Database name

T=$"TrackTable"; # Table name in SQLite

for file in $(ls *.gpx); do  # Loop through all the GPX files in the directory
ROW=$(spatialite $DB "SELECT COUNT(ROWID) FROM $T"); # What was the last row?
ogr2ogr -append -f "SQLite" -dsco SPATIALITE=yes # add to table with Spatialite filetype
    -dsco INIT_WITH_EPSG=yes -t_srs epsg:4326 # using WGS84
    $DB "$file" tracks -nln $T; # adding tracks
spatialite $DB "UPDATE $T SET trackname='$file' WHERE ROWID &gt; $ROW"; # update trackname for all the new rows

This now gives you a Spatialite table with all the tracks – including saved tracks, not just the ‘ACTIVE LOG’ – as rows in the database. What else do we need?

  • It would be good to do the same thing for all the waypoints (as a separate table);
  • Perhaps save the GPX file itself as a BLOB.

QGIS with R: Working with the SEXTANTE plugin

SEXTANTE toolboxI’ve tried to get QGIS and R playing together in the past, with only occasional success. The new SEXTANTE plugin for QGIS, which provides links into a number of applications including R, SAGA and GDAL, looks really promising.

Getting it working on Mac OSX was a little tricky, as there is a (known) bug in the plugin that hasn’t been resolved yet.

The fix:

  1. Open up the python plugin directory for SEXTANTE: ~/.qgis/python/plugins/sextante/r/
  2. Make a backup copy of (command-D in the Finder);
  3. Open up in a text or code editor (I’m using Komodo edit these days and it’s working really well);
  4. On line 58, change:
    command = ["R", "CMD", "BATCH", "--vanilla", RUtils.getRScriptFilename(), RUtils.getConsoleOutputFilename()
    command = "R CMD BATCH --vanilla " + RUtils.getRScriptFilename() + " "+ RUtils.getConsoleOutputFilename()
  5. Save and restart QGIS.

While you can’t simply run an R script directly, you can create a .rsx file either in an editor or directly with the SEXTANTE toolbox and use it. This is the Raster_histogram.rsx file from the examples:

##[Example scripts]=group
##layer = raster

Giving, from a DEM layer in QGIS:

Summarising QGIS project layers

With sometimes dozens of layers in a QGIS project, consolidating everything into one directory can be really handy (Plugins | QConsolidate). I’m finding this really useful when sending a project file to another device (like an Android tablet running QGIS) or to a colleague.

Multiple layers with varying projections is a big performance killer in QGIS. Plus, it looks like there’s a serious bug using a GPS connection with non-WGS84 projections – so, I started looking at a current project (Luangwa Valley aerial survey) and removing or reprojecting layers so that I had a single set of projections and could turn off on-the-fly projections.

Tedious. Rick-click and look at CRS for every one of 25 layers? I really need an overview of layers, data sources and CRS for every layer.

I’ve tried to use the XML package in R before without much success – this time I made sure to go through the excellent XPath tutorial on and found some worked examples on stackoverflow.

The result, a function that takes a QGIS project file (an XML document) and spits out a dataframe with columns for the layer name, CRS and location of the data source.


QGIS.summary <- function (fn) {
# Get the QGIS project file and parse it:
file <- xmlParse(fn)

# Main doc called 'qgis' and the layers are in the 'maplayer' level;
# we need the srid and the layername xmlValues:
file.CRS <- xpathApply(file, "/qgis//maplayer/srs/spatialrefsys/description", xmlValue, "description")
file.layers <- xpathApply(file, "/qgis//maplayer/layername", xmlValue, "layername")
file.datasource <- xpathApply(file, "/qgis//maplayer/datasource", xmlValue, "datasource")

# Return a dataframe:
return(data.frame(layers = as.character(file.layers), CRS = as.character(file.CRS), sourcefiles = as.character(file.datasource)))

QGIS on Android phone: Installing

There is a more up-to-date version of this page detailing the installation on an Asus TF300T tablet here.

I decided to try out QGIS on my Samsung Galaxy Advance – a little 4-inch screen, but why not? I’d run through the same process on an Asus TF300 tablet previously (much larger screen and more appropriate for GIS use – but more on that later) and wanted to document how it went.

You’ll need an active and fast network connection – you’ll download a bit over 100MB of files during the installation. (For more screenshots and installation details see here from


  1. To install ensure you’ve selected “install from unknown sources” in the Settings;
  2. Download the installer .apk and open it (on most Androids it will download to a “Downloads” directory, but may be elsewhere; you can just click to open the file from your file browser and it should start the installer).
  3. Run the installer and select ‘download and install’. 83 MB of data will get downloaded. Luckily, it is a resumable download – if it fails, you can start the process off again where you left off.
  4. Confirm that you want to install (standard Android dialogue) – takes about 2 minutes to finish this stage.
  5. Run QGIS: “Unpacking post-install data” … 10 seconds.
  6. QGIS requires a supporting service, available from the Play store “This application requires Ministro service. Would you like to install it?” Confirm and install Ministro II, only 523KB.
  7. Looks like it needs MORE libraries to run! “Qgis needs extra libraries to run. Do you want to download them now?” Confirm and install the QtCore libraries – a lot bigger at 31MB.

Then QGIS just seems to run! Tough to use the little tiny buttons / menus, but workable. You can even use the GPS connection – took me a while to realise it was under View | Panels | GPS – click ‘connect’ on the panel after displaying it.



East Africa’s latest country: 坦桑尼亚

Google Maps seems to have a bit of a labelling issue this morning. Here’s Tanzania:

Well, it had to happen eventually. Let’s just hope it’s pronounced the same way.

Even better:

Aerial photo examples

OpenStreetMap development in Tanzania: Tandale

OpenStreetMap is a great alternative to Google Maps – community developed and maintained, and the data are freely available for use in your own mapping projects. Good recent discussion about OSM and GMaps on GIS.Stackexchange.

I favour OSM:

  1. Editing is far more straightforward, particularly with a large area that is previously unmapped, or (worse yet) badly mapped. JOSM, QGIS and other tools let you update large areas, correct or add data in bulk. Working on areas here in Tanzania that are poorly covered by OSM and Google (click links to see example of same area), OSM is far, far better for editing.
  2. With this comes better access to data for offline mapping. While both APIs give you great access to maps (and Google’s gives you more options and are prettier), OSM can be downloaded monthly for any area of interest. Google’s given us access to road data (shapefiles) for Tanzania before but it takes time.

(1) in particular has been annoying recently. I have battled with Google editors for approval of known bad data, or not had my GPS-derived roads approved because it doesn’t match with the 8-year old satellite data! The Google editors who approve changes to the Mapmaker data don’t have any ground experience in the local areas you might edit, relying on what they can see in the GMaps satellite layer.

You can edit OSM quickly and to great effect. The ‘ramaniTanzania’ project recently mapped Tandale. Here’s what OSM had before:

Tandale map on 9th August 2011 - bare-bones map

And after eleven days of mapping and editing:

Tandale map after eleven days of editing, with far more detail

I doubt you could do anything like that, in that amount of time, with Google Maps. Great job!

The project also has an Ushahidi installation for community reports:

This is similar to the tracker TZGISUG is trialling for Parthenium hysterophorus.

For comparison, here’s a screenshot of Google Maps for the area: