Tag Archives: Python

PyQGIS 101: Introduction to QGIS Python programming for non-programmers has now reached the part 10 milestone!

Beyond the obligatory Hello world! example, the contents so far include:

If you’ve been thinking about learning Python programming, but never got around to actually start doing it, give PyQGIS101 a try.

I’d like to thank everyone who has already provided feedback to the exercises. Every comment is important to help me understand the pain points of learning Python for QGIS.

I recently read an article – unfortunately I forgot to bookmark it and cannot locate it anymore – that described the problems with learning to program very well: in the beginning, it’s rather slow going, you don’t know the right terminology and therefore don’t know what to google for when you run into issues. But there comes this point, when you finally get it, when the terminology becomes clearer, when you start thinking “that might work” and it actually does! I hope that PyQGIS101 will be a help along the way.


We’ve seen a lot of explorative movement data analysis in the Movement data in GIS series so far. Beyond exploration, predictive analysis is another major topic in movement data analysis. One of the most obvious movement prediction use cases is trajectory prediction, i.e. trying to predict where a moving object will be in the future. The two main categories of trajectory prediction methods I see are those that try to predict the actual path that a moving object will take versus those that only try to predict the next destination.

Today, I want to focus on prediction methods that predict the path that a moving object is going to take. There are many different approaches from simple linear prediction to very sophisticated application-dependent methods. Regardless of the prediction method though, there is the question of how to evaluate the prediction results when these methods are applied to real-life data.

As long as we work with nice, densely, and regularly updated movement data, extracting evaluation samples is rather straightforward. To predict future movement, we need some information about past movement. Based on that past movement, we can then try to predict future positions. For example, given a trajectory that is twenty minutes long, we can extract a sample that provides five minutes of past movement, as well as the actually observed position five minutes into the future:

But what if the trajectory is irregularly updated? Do we interpolate the positions at the desired five minute timestamps? Do we try to shift the sample until – by chance – we find a section along the trajectory where the updates match our desired pattern? What if location timestamps include seconds or milliseconds and we therefore cannot find exact matches? Should we introduce a tolerance parameter that would allow us to match locations with approximately the same timestamp?

Depending on the duration of observation gaps in our trajectory, it might not be a good idea to simply interpolate locations since these interpolated locations could systematically bias our evaluation. Therefore, the safest approach may be to shift the sample pattern along the trajectory until a close match (within the specified tolerance) is found. This approach is now implemented in MovingPandas’ TrajectorySampler.

def test_sample_irregular_updates(self):
    df = pd.DataFrame([
        {'geometry':Point(0,0), 't':datetime(2018,1,1,12,0,1)},
        {'geometry':Point(0,3), 't':datetime(2018,1,1,12,3,2)},
        {'geometry':Point(0,6), 't':datetime(2018,1,1,12,6,1)},
        {'geometry':Point(0,9), 't':datetime(2018,1,1,12,9,2)},
        {'geometry':Point(0,10), 't':datetime(2018,1,1,12,10,2)},
        {'geometry':Point(0,14), 't':datetime(2018,1,1,12,14,3)},
        {'geometry':Point(0,19), 't':datetime(2018,1,1,12,19,4)},
        {'geometry':Point(0,20), 't':datetime(2018,1,1,12,20,0)}
    geo_df = GeoDataFrame(df, crs={'init': '4326'})
    traj = Trajectory(1,geo_df)
    sampler = TrajectorySampler(traj, timedelta(seconds=5))
    past_timedelta = timedelta(minutes=5)
    future_timedelta = timedelta(minutes=5)
    sample = sampler.get_sample(past_timedelta, future_timedelta)
    result = sample.future_pos.wkt
    expected_result = "POINT (0 19)"
    self.assertEqual(result, expected_result)
    result = sample.past_traj.to_linestring().wkt
    expected_result = "LINESTRING (0 9, 0 10, 0 14)"
    self.assertEqual(result, expected_result)

The repository also includes a demo that illustrates how to split trajectories using a grid and finally extract samples:


Need to geocode some addresses? Here’s a five-lines-of-code solution based on “An A-Z of useful Python tricks” by Peter Gleeson:

from geopy import GoogleV3
place = "Krems an der Donau"
location = GoogleV3().geocode(place)

For more info, check out geopy:

geopy is a Python 2 and 3 client for several popular geocoding web services.
geopy includes geocoder classes for the OpenStreetMap Nominatim, ESRI ArcGIS, Google Geocoding API (V3), Baidu Maps, Bing Maps API, Yandex, IGN France, GeoNames, Pelias,, OpenMapQuest, PickPoint, What3Words, OpenCage, SmartyStreets, GeocodeFarm, and Here geocoder services.

If you’re are following me on Twitter, you’ve certainly already read that I’m working on PyQGIS 101 a tutorial to help GIS users to get started with Python programming for QGIS.

I’ve often been asked to recommend Python tutorials for beginners and I’ve been surprised how difficult it can be to find an engaging tutorial for Python 3 that does not assume that the reader already knows all kinds of programming concepts.

It’s been a while since I started programming, but I do teach QGIS and Python programming for QGIS to university students and therefore have some ideas of which concepts are challenging. Nonetheless, it’s well possible that I overlook something that is not self explanatory. If you’re using PyQGIS 101 and find that some points could use further explanations, please leave a comment on the corresponding page.

PyQGIS 101 is a work in progress. I’d appreciate any feedback, particularly from beginners!

This post describes how to calculate raster profile information (altitude changes to be exact) for a road vector layer from a DEM. DEM source is NASA’s free SRTM data.

The following script is based on scw’s answer to my related question on gis.stackexchange. This script is supposed to be run from inside a GRASS console (as I haven’t figured out yet how to import grass into python otherwise).

The idea behind this is to create points along the input road vectors and then sample the DEM values at those points. The resulting 3D points are then evaluated in a python function and the results inserted back into GRASS. There, the new values (altitude changes up and down hill) are joined back to the road elements.

import grass.script as grass
from pysqlite2 import dbapi2 as sqlite
from datetime import datetime

class AltitudeCalculator():
    def __init__(this, point_file):
        this.file = open(point_file,'r')
        this.prev_elevation = None
        this.prev_id = None
        this.increase = 0
        this.decrease = 0 = None
        this.output = ''
    def calculate(this):
        for line in this.file:
            line = line.split('|')
            #x = float(line[0])
            #y = float(line[1])
            elev = float(line[2]) 
   = int(line[3])
            if this.prev_id:
                if this.prev_id ==
                    # continue this road id
                    change = elev - this.prev_elevation 
                    if change > 0:
                        this.increase += change
                        this.decrease += change   
                    # end of road id, save current data and reset
            this.prev_id =
            this.prev_elevation = elev
        this.addOutputLine() # last entry  

    def reset(this):
        this.increase = 0
        this.decrease = 0 
    def addOutputLine(this):
        this.output += "%i,%i,%i\n" % (this.prev_id, int(round(this.increase)), int(round(this.decrease)))
    def save(this, file_name):
        # create .csv file
        file = open(file_name,'w')
        # create .csvt file describing the columns
        file = open(file_name+'t','w')
def main():
    road_file = '/home/user/maps/roads.shp'
    elevation_file = '/home/user/maps/N48E016.hgt'
    point_file = '/home/user/maps/osm_road_pts_3d.txt'
    alt_diff_file = '/home/user/maps/alt_diff.csv'
    db_file = '/home/user/grass/newLocation/test/sqlite.db'
    print('db.connect: '+str(
    grass.run_command("db.connect", driver='sqlite', database=db_file)
    # import files    
    print('import files: '+str(
    grass.run_command("", input=elevation_file, output='srtm_dem', overwrite=True)
    grass.run_command("", dsn=road_file, output='osm_roads', min_area=0.0001, snap=-1, overwrite=True)
    # set resolution to match DEM
    print('set resolution to match DEM: '+str(
    grass.run_command("g.region", rast='srtm_dem')

    # convert to points
    print('convert to points: '+str(
    grass.run_command("",  input='osm_roads', output='osm_road_pts', type='line', dmax=0.001, overwrite=True, flags='ivt')
    # extract elevation at each point
    print('extract elevation at each point: '+str(
    grass.run_command("v.drape", input='osm_road_pts', output='osm_road_pts_3d', type='point', rast='srtm_dem', method='cubic', overwrite=True)
    # export to text files
    print('export to text files: '+str(
    grass.run_command("v.out.ascii", input='osm_road_pts_3d', output=point_file, overwrite=True)
    # calculate height differences
    print('calculate height differences: '+str(
    calculator = AltitudeCalculator(point_file)
    # import height differences into grass
    print('import height differeces into grass: '+str(
    grass.run_command("db.droptable", table='alt_diff', flags='f')
    grass.run_command("", dsn=alt_diff_file, output='alt_diff')
    # create an index
    print('create an index: '+str(
    connection = sqlite.connect(db_file) 
    cursor = connection.cursor()
    sql = 'create unique index alt_diff_road_id on alt_diff (road_id)'
    # join
    print('join: '+str(
    grass.run_command("v.db.join", map='osm_roads', layer=1, column='cat', otable='alt_diff', ocolumn='road_id')
    print('Done: '+str(

if __name__ == "__main__":
    main() hosts a really neat Python script for heatmap creation.

The script takes data points with latitude and longitude coordinates (list of points or a GPX track log) and plots them on a map so areas where points are dense show brightly.

Nice extras include: GPX tracks can be rendered as line segments instead of disconnected points. You can generate animations using ffmpeg. The heatmaps can be put on top of OpenStreetMap tiles.

The script is released under AGPL making it free to improve and share.

Maybe this is a bit off-topic, but I just spent quite some time on this and I need to write it down so I can look it up again later :)

These are instructions for Ubuntu running Postgres 8.4. By default, Postgres ships without PL/Python so we need to get it first:

sudo apt-get install postgresql-plpython-8.4

Next, we need to create the language for our database. I’m using PgAdmin3. From there, I ran:

CREATE PROCEDURAL LANGUAGE 'plpython' HANDLER plpython_call_handler;

This should have been it. Let’s try with a simple function:

CREATE FUNCTION replace_e_to_a(text) RETURNS text AS
import re
Text1 = re.sub(''e'', ''a'',args[0])
return Text1
LANGUAGE 'plpython';

SELECT replace_e_to_a('Meee');

… should return ‘Maaa’.

Now for the juicy part: Let’s create an INSERT trigger function!

First, let’s have a look at the corresponding table structure. We have two tables “user_data” and “user_keywords”. “user_data” is the table that’s being filled with information from external functions. “user_keywords” has to be kept up-to-date. It is supposed to count the appearance of keywords on a per-user base.

user_data                                   user_keywords
user_id, event_id, keywords                 user_id, keyword,   count
1,       1,        'music,rock'             1,       'music',   2
1,       2,        'music,classic'          1,       'rock',    1
                                            1,       'classic', 1

First, the keyword list has to be split. Then a row has to be inserted for new keywords (compare insert_plan) and the counter has to be increased for existing keywords (update_plan).

The values that are about to be inserted can be accessed via TD[“new”][“column_name”].

CREATE FUNCTION update_keyword_count() RETURNS trigger AS '

keywords = TD["new"]["keywords"]
user = TD["new"]["user_id"]

insert_plan = plpy.prepare("INSERT INTO user_keywords (keyword, count, user_id) VALUES ($1, $2, $3)", ["text", "int", "int"])

update_plan = plpy.prepare("UPDATE user_keywords SET count = $3 WHERE user_id = $1 AND keyword = $2", ["int", "text", "int"])

for keyword in keywords.split(","):
  select_cnt_rows = plpy.prepare("SELECT count(*) AS cnt FROM user_keywords WHERE user_id = $1 AND keyword = $2", ["int", "text"])
  cnt_rows = plpy.execute(select_cnt_rows, [user, keyword])

  select_plan = plpy.prepare("SELECT count AS cnt FROM user_keywords WHERE user_id = $1 AND keyword = $2", ["int", "text"])
  results = plpy.execute(select_plan, [user, keyword])

  if cnt_rows[0]["cnt"] == 0:
   rv = plpy.execute(insert_plan, [keyword, 1, user])
   rv = plpy.execute(update_plan, [user, keyword, results[0]["cnt"]+1])

' LANGUAGE plpython;

Now, we need to wire it up by defining the trigger:

CREATE TRIGGER update_keywords
EXECUTE PROCEDURE update_keyword_count();

… Wasn’t that bad ;)

%d bloggers like this: