Adding elevation data to a point layer from a DEM

I’ve put together a Python script that uses pyqgis to add elevation data from a DEM to the attribute table of a point layer.  I’m going to use it to compare the elevations recorded by my Edge 500 to elevations derived from local DEMs.  I pieced the script together with tons of help from the PyQGIS Cookbook.

Here’s what I came up with.  It requires both the point layer and DEM to be in the same projection.

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
#

from qgis.core import *
import qgis.utils

#initializes QGIS and points python to where it can find QGIS's stuff
QgsApplication.setPrefixPath("/usr/", True)
QgsApplication.initQgis()

vfileName = "/path/to/point_layer"

#QgsVectorLayer(file name, display name, data provider) 
vlayer = QgsVectorLayer(vfileName, "elev_pts", "ogr")
if not vlayer.isValid():
    print "Vector did not load"

rfileName = "/path/to/dem_raster"

#QgsRasterLayer(file name, display name)
rlayer = QgsRasterLayer(rfileName, "RedBrushDEM")
if not rlayer.isValid():
    print "Raster did not load"

provider = vlayer.dataProvider()
feat = QgsFeature()
allAttrs = provider.attributeIndexes()

#fetch all geometries and all attributes.  Can be filtered, if needed
provider.select(allAttrs)

#Loop through each vector feature
while provider.nextFeature(feat):
    geom = feat.geometry()
    x = geom.asPoint()

    #Get the raster value of the cell under the vector point
    res, ident = rlayer.identify(QgsPoint(x))

    #The key is the raster band, and the value is the cell's raster value
    for (k,v) in ident.iteritems():
        elevation = float(v)

    fid = int(feat.id())
    vlayer.startEditing()
    #changeAttributeValue(feature id, index of field to be changed, new attribute value)
    vlayer.changeAttributeValue(fid, 2, elevation)
    vlayer.commitChanges()
#Clean up
QgsApplication.exitQgis()
Advertisements
    • Kirubaharan
    • July 11th, 2013

    Can you please suggest a good IDE/editor compatible with pyqgis in linux/ubuntu ?

    • I mainly used PyScripter on Windows and gvim on Linux with the PyQGIS cookbook and the API open in a browser on both.

  1. No trackbacks yet.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: