# 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!

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)
#y = float(line)
elev = float(line)
this.id = int(line)

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.reset()

this.prev_id = this.id
this.prev_elevation = elev

this.addOutputLine() # last entry

def reset(this):
this.increase = 0
this.decrease = 0

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.close()
# create .csvt file describing the columns
file = open(file_name+'t','w')
file.write('"Integer","Integer","Integer"')
file.close()

def main():
elevation_file = '/home/user/maps/N48E016.hgt'
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)

# 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()))

print('Done: '+str(datetime.now()))

if __name__ == "__main__":
main()
```

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)
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["cnt"] == 0:
rv = plpy.execute(insert_plan, [keyword, 1, user])
else:
rv = plpy.execute(update_plan, [user, keyword, results["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 ;)

Interesting news from the gispython.org community mailing list: There is a new Python library out there called PyKML.

PyKML allows parsing and authoring of KML documents based on the lxml.objectify API which provides Pythonic access to XML documents.

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
```

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

"""read plugin settings from QgsProject instance"""
# map data types to function names
prj = QgsProject.instance()
functions = { 'str' : prj.readEntry,