Advertisements

Archive

Tag Archives: Python

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!

Advertisements

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
        this.id = 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]) 
            this.id = int(line[3])
        
            if this.prev_id:
                if this.prev_id == this.id:
                    # continue this road id
                    change = elev - this.prev_elevation 
                    if change > 0:
                        this.increase += change
                    else:
                        this.decrease += change   
                else:
                    # end of road id, save current data and reset
                    this.addOutputLine()
                    this.reset()
                    
            this.prev_id = this.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')
        file.write('road_id,sum_up,sum_down\n'+this.output) 
        file.close()
        # create .csvt file describing the columns
        file = open(file_name+'t','w')
        file.write('"Integer","Integer","Integer"')
        file.close()
        
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(datetime.now()))
    grass.run_command("db.connect", driver='sqlite', database=db_file)
    
    # import files    
    print('import files: '+str(datetime.now()))
    grass.run_command("r.in.gdal", input=elevation_file, output='srtm_dem', overwrite=True)
    grass.run_command("v.in.ogr", 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(datetime.now()))
    grass.run_command("g.region", rast='srtm_dem')

    # convert to points
    print('convert to points: '+str(datetime.now()))
    grass.run_command("v.to.points",  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(datetime.now()))
    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(datetime.now()))
    grass.run_command("v.out.ascii", input='osm_road_pts_3d', output=point_file, overwrite=True)
    
    # calculate height differences
    print('calculate height differences: '+str(datetime.now()))
    calculator = AltitudeCalculator(point_file)
    calculator.calculate()
    calculator.save(alt_diff_file)
    
    # import height differences into grass
    print('import height differeces into grass: '+str(datetime.now()))
    grass.run_command("db.droptable", table='alt_diff', flags='f')
    grass.run_command("db.in.ogr", dsn=alt_diff_file, output='alt_diff')
    
    # create an index
    print('create an index: '+str(datetime.now()))
    connection = sqlite.connect(db_file) 
    cursor = connection.cursor()
    sql = 'create unique index alt_diff_road_id on alt_diff (road_id)'
    cursor.execute(sql)
    connection.close()
    
    # join
    print('join: '+str(datetime.now()))
    grass.run_command("v.db.join", map='osm_roads', layer=1, column='cat', otable='alt_diff', ocolumn='road_id')
    
    print('Done: '+str(datetime.now()))

if __name__ == "__main__":
    main()

http://www.sethoscope.net 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])
  else:
   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
BEFORE INSERT ON user_data
FORE EACH ROW
EXECUTE PROCEDURE update_keyword_count();

… Wasn’t that bad ;)

For some plugins it’s necessary to save the settings inside the QGIS project file. Saving is done with this simple one-liner:

QgsProject.instance().writeEntry(pluginName, setting, value)

Then you just need to save the project.

Reading is performed with one of the following functions:

QgsProject.instance().readEntry (pluginName, setting) # for strings
QgsProject.instance().readNumEntry (pluginName, setting)
QgsProject.instance().readDoubleEntry (pluginName, setting)
QgsProject.instance().readBoolEntry (pluginName, setting)
QgsProject.instance().readListEntry (pluginName, setting)

You’ll find the corresponding API documentation at: http://doc.qgis.org/stable/classQgsProject.html. As you can see, you can only read/write simple data types. To allow the plugin developer to save more complex plugin settings, I filed an enhancement request.

To handle all those different read functions in a convenient way, I created the following functions:

def readSetting(self,func,setting):
    """read a plugin setting from QgsProject instance"""
    value,ok = func('pluginName',setting)
    if ok:
        return value
    else:
        return None
            
def readSettings(self,setting,value):
    """read plugin settings from QgsProject instance"""
    # map data types to function names
    prj = QgsProject.instance()
    functions = { 'str' : prj.readEntry,
                  'int' : prj.readNumEntry,
                  'float' : prj.readDoubleEntry,
                  'bool' : prj.readBoolEntry,
                  'pyqtWrapperType' : prj.readListEntry # QStringList
                }
        
    dataType = type(value).__name__
    return = self.readSetting(functions[dataType],setting)

readSettings() has to be supplied with the name of the setting and an example or default value for the setting (from which we can determine the data type of the setting). Of course this can be done in many different ways. In Time Manager plugin, readSettings() receives a dictionary of all settings that have to be read. The function then loops through the dictionary and reads the available settings.

%d bloggers like this: