Contents

  1. Getting PyQGIS to work
    1. Manual compilation
    2. Troubleshooting
  2. Using PyQGIS
    1. Using raster layer
    2. Using vector layer
    3. Geometry handling
    4. Spatial references
    5. Coordinate transformations
    6. Using spatial index
    7. Rendering
    8. Measuring
    9. Writing shapefiles
    10. Memory provider
  3. Creating bindings

Starting from 0.9 release, QGIS has optional scripting support using Python language. We've decided for Python as it's one of the most favourite languages for scripting. PyQGIS bindings depend on SIP and PyQt4. The reason for using SIP instead of more widely used SWIG is that the whole QGIS code depends on Qt libraries. Python bindings for Qt (PyQt) are done also using SIP and this allows seamless integration of PyQGIS with PyQt.

Getting PyQGIS to work

When using binary packages from official QGIS download page you should get Python bindings automatically installed with QGIS >= 0.9.

If you open QGIS you should find in menu Plugins an item "Python console" - if not, you don't have Python support installed.

Manual compilation

PyQGIS has following dependencies:

In case your distribution doesn't contain packages fresh enough, you'll need to compile them manually. Windows users can use PyQt binary package provided on PyQt download page, it also contains SIP python module.

Now the compilation: first you need to make sure that CMake variable WITH_BINDINGS turned on (it's on by default) and SIP and PyQt4 must have been successfully detected. During configuration CMake will tell whether there are any problems.

By default it installs PyQGIS to QGIS_PREFIX/share/qgis/python (or QGIS_PREFIX/python on windows). So if you plan to use bindings in some custom scripts, you'll need to set PYTHONPATH environment variable to the one where bindings reside. In case you want python bindings to be installed to global python site-packages directory, set BINDINGS_GLOBAL_INSTALL variable in CMake to TRUE. For this type of installation might need root privileges, however setting PYTHONPATH won't be necessary.

Troubleshooting

Using PyQGIS

There are several ways how to use QGIS python bindings:

This wiki page is intended mainly for creating 3rd party apps by providing some useful code snippets. There's also complete QGIS API reference which also includes classes that are not part of QGIS libraries. Pythonic QGIS API is 99% identical to the API in C++.

There are some resources about programming with PyQGIS on QGIS blog. See QGIS tutorial ported to Python for some examples of simple 3rd party apps. A good resource when dealing with plugins is to download some plugins from plugin repository and examine their code.

Notes:

First of all you have to import qgis module, set QGIS path where to search for resources - database of projections, providers etc. When you set prefix path with second argument set as True, QGIS will initialize all paths with standard dir under the prefix directory. Calling initQgis() function is important to let QGIS search for the available providers. 

   1 from qgis.core import *
   2 
   3 # supply path to where is your qgis installed
   4 QgsApplication.setPrefixPath("/path/to/qgis/installation", True)
   5 
   6 # load providers
   7 QgsApplication.initQgis()

When QGIS library is initialized, let's create some layers. If you would like to use them for rendering, don't forget to add them to map layer registry, which takes ownership of them. Map layer registry is later used for retreiving layers by their ID.

It's good to check whether vector/raster layer loaded correctly as it's not possible to find this out when creating the object: 

   1 if vlayer.isValid():
   2    print "Layer loaded."
   3 else:
   4    print "Failed to load layer!"

In case you plan to do some rendering, don't forget to add layer to the registry. If you want to just get read some data from it or write some data, you don't need to add it to map layer registry. 

   1 QgsMapLayerRegistry.instance().addMapLayer(vlayer)

When you're done with using QGIS library, call exitQgis() to make sure that everything is cleaned up (e.g. clear map layer registry and delete layers): 

   1 QgsApplication.exitQgis()

Using raster layer

/!\ Note: QgsRasterLayer.identify() is available from QGIS 0.11.

To do a query on value of bands of raster layer at some specified point: 

   1 ident = rlayer.identify(QgsPoint(15.30,40.98))
   2 for (k,v) in ident.iteritems():
   3         print str(k),":",str(v)

The identify function returns a dictionary - keys are band names, values are the values at chosen point. Both key and value are QString instances so to see actual value you'll need to convert them to python strings (as shown in code snippet).

Using vector layer

Below is an example how to go through the features of the layer. To read features from layer, use getNextFeature() calls. 

   1 provider = vlayer.getDataProvider()
   2 
   3 feat = QgsFeature()
   4 allAttrs = provider.allAttributesList()
   5 
   6 # start data retreival: fetch geometry and all attributes for each feature
   7 provider.select(allAttrs)
   8 
   9 # retreive every feature with its geometry and attributes
  10 while provider.getNextFeature(feat):
  11 
  12    # fetch geometry
  13    geom = feat.geometry()
  14    print "Feature ID %d: " % feat.featureId() ,
  15 
  16    # show some information about the feature
  17    if geom.vectorType() == QGis.Point:
  18      x = geom.asPoint()
  19      print "Point: " + str(x)
  20    elif geom.vectorType() == QGis.Line:
  21      x = geom.asPolyline()
  22      print "Line: %d points" % len(x)
  23    elif geom.vectorType() == QGis.Polygon:
  24      x = geom.asPolygon()
  25      numPts = 0
  26      for ring in x:
  27        numPts += len(ring)
  28      print "Polygon: %d rings with %d points" % (len(x), numPts)
  29    else:
  30      print "Unknown"
  31 
  32    # fetch map of attributes
  33    attrs = feat.attributeMap()
  34 
  35    # attrs is a dictionary: key = field index, value = QgsFeatureAttribute
  36    # show all attributes and their values
  37    for (k,attr) in attrs.iteritems():
  38       print "%d: %s" % (k, attr.toString())

select() gives you flexibility in what data will be fetched. It can get 4 arguments, all of them are optional:

No.

Name

Explanation

1

fetchAttributes

List of attributes which should be fetched. Default: empty list

2

rect

Spatial filter. If empty rect is given (QgsRect(), all features are fetched. Default: empty rect

3

fetchGeometry

Whether geometry of the feature should be fetched. Default: True

4

useIntersect

When using spatial filter, this argument says whether accurate test for intersection should be done or whether test on bounding box suffices. This is needed e.g. for feature identification or selection. Default: False

Some examples: 

   1 # fetch features with geometry and only first two fields
   2 provider.select([0,1])
   3 
   4 # fetch features with geometry which are in specified rect, attributes won't be retreived
   5 provider.select([], QgsRect(23.5, -10, 24.2, -7))
   6 
   7 # fetch features without geometry, with all attributes
   8 allAtt = provider.allAttributesList()
   9 provider.select(allAtt, QgsRect(), False)

To obtain field index from its name, use provider's indexFromFieldName function: 

   1 fldDesc = provider.indexFromFieldName("DESCRIPTION")
   2 if fldDesc == -1:
   3   print "Field not found!"

Geometry handling

To extract information from geometry there are accessor functions for every vector type. How do accessors work:

asPoint()

returns QgsPoint

asPolyline()

returns a list of QgsPoints

asPolygon()

returns a list of rings, every ring consists of a list of QgsPoints. First ring is outer, subsequent rings are holes

You can use QgsGeometry.isMultipart() to find out whether the feature is multipart. For multipart features there are similar accessor functions: asMultiPoint(), asMultiPolyline(), asMultiPolygon().

There are several options how to create a geometry:

Spatial references

Spatial reference systems (SRS) are encapsulated by class QgsSpatialRefSys. Instances of this class can be created by several different ways:

It's wise to check whether creation (i.e. lookup in the database) of the SRS has been successful: isValid() must return True.

Note that for initialization of spatial reference systems QGIS needs to lookup appropriate values in its internal database srs.db. Thus in case you create an independent application you need to set paths correctly with QgsApplication.setPrefixPath() otherwise it will fail to find the database. If you're developping a plugin you don't care: everything is already set up for you.

Accessing spatial reference system information: 

   1 print "QGIS SRS ID:", srs.srsid()
   2 print "PostGIS SRID:", srs.srid()
   3 print "EPSG ID:", srs.epsg()
   4 print "Description:", srs.description()
   5 print "Projection Acronym:", srs.projectionAcronym()
   6 print "Ellipsoid Acronym:", srs.ellipsoidAcronym()
   7 print "Proj4 String:", srs.proj4String()
   8 # check whether it's geographic or projected coordinate system
   9 print "Is geographic:", srs.geographicFlag()
  10 # check type of map units in this SRS (values defined in QGis::units enum)
  11 print "Map units:", srs.mapUnits()

Coordinate transformations

You can do transformation between different spatial reference systems by using QgsCoordinateTransform class. The easiest way to use it is to create source and destination SRS and construct QgsCoordinateTransform instance with them. Then just repeatedly call transform() function to do the transformation. By default it does forward transformation, but it's capable to do also inverse transformation. 

   1 srsSrc = QgsSpatialRefSys(4326)    # WGS 84
   2 srsDest = QgsSpatialRefSys(32633)  # WGS 84 / UTM zone 33N
   3 xform = QgsCoordinateTransform(srsSrc,srsDest)
   4 
   5 # forward transformation: src -> dest
   6 pt1 = xform.transform(QgsPoint(18,5))
   7 print "Transformed point:", pt1
   8 
   9 # inverse transformation: dest -> src
  10 pt2 = xform.transform(pt1, QgsCoordinateTransform.INVERSE)
  11 print "Transformed back:", pt2

Using spatial index

  1. create spatial index - the following code creates an empty index 
    •    1 index = QgsSpatialIndex()

  2. add features to index - index takes QgsFeature object and adds it to the internal data structure. You can create the object manually or use one from previous call to provider's getNextFeature() 

    •    1 index.insertFeature(feat)
  3. once spatial index is filled with some values, you can do some queries 
    •    1 # returns array of feature IDs of five nearest features
   2 nearest = index.nearestNeighbor(QgsPoint(25.4, 12.7), 5)
   3 
   4 # returns array of IDs of features which intersect the rectangle
   5 intersect = index.intersects(QgsRect(22.5, 15.3, 23.1, 17.2))

Rendering

Render some layers using QgsMapRender - create destination paint device (QImage, QPainter etc.), set up layer set, extent, output size and do the rendering: 

   1 
   2 # create image
   3 img = QImage(QSize(800,600), QImage.Format_ARGB32_Premultiplied)
   4 
   5 # set image's background color
   6 color = QColor(255,255,255)
   7 img.fill(color.rgb())
   8 
   9 # create painter
  10 p = QPainter()
  11 p.begin(img)
  12 p.setRenderHint(QPainter.Antialiasing)
  13 
  14 render = QgsMapRender()
  15 
  16 # set layer set
  17 lst = [ layer.getLayerID() ]  # add ID of every layer
  18 render.setLayerSet(lst)
  19 
  20 # set extent
  21 rect = QgsRect(render.fullExtent())
  22 rect.scale(1.1)
  23 render.setExtent(rect)
  24 
  25 # set output size
  26 render.setOutputSize(img.size(), img.logicalDpiX())
  27 
  28 # do the rendering
  29 render.render(p)
  30 
  31 p.end()
  32 
  33 # save image
  34 img.save("render.png","png")

Measuring

To compute distances or areas, use QgsDistanceArea class. If projections are turned off, calculations will be planar, otherwise they'll be done on ellipsoid. When ellipsoid is not set explicitly it uses WGS84 parameters for calculations. 

   1 d = QgsDistanceArea()
   2 d.setProjectionsEnabled(True)
   3 
   4 print "distance in meters: ", d.measureLine(QgsPoint(10,10),QgsPoint(11,11))

Writing shapefiles

You can write shapefiles using QgsVectorFileWriter class. At some point it should support any kind of vector file that OGR supports, so far it's limited only to shapefiles. There are two possibilities how to export a shapefile:

Memory provider

/!\ This functionality has been added in QGIS 0.11

Memory provider is intended to be used mainly by plugin or 3rd party app developers. It doesn't store data on disk, allowing developers to use it as a fast backend for some temporary layers (until now one had to use e.g. OGR provider). You can use it by passing "memory" provider string to vector layer constructor.

Following URIs are allowed: "Point" / "LineString" / "Polygon" / "MultiPoint" / "MultiLineString" / "MultiPolygon" - for different types of data stored in the layer.

When adding fields you can use following field type names: "string" / "int" / "double"

Following example should be self-explaining: 

   1 # create layer
   2 vl = QgsVectorLayer("Point", "name for legend", "memory")
   3 pr = vl.getDataProvider()
   4 
   5 # add fields
   6 # to preserve correct order they must be added one-by-one (d'oh)
   7 pr.addAttributes( { "name" : "string" } )
   8 pr.addAttributes( { "age" : "int" } )
   9 pr.addAttributes( { "size" : "double" } )
  10 
  11 # add a feature
  12 fet = QgsFeature()
  13 fet.setGeometry(QgsGeometry.fromPoint(QgsPoint(10,10)))
  14 fet.setAttributeMap( { 0 : QVariant("Johny"), 1 : QVariant(20), 2 :
  15 QVariant(0.3) } )
  16 pr.addFeatures( [ fet ] )
  17 
  18 # update layer's extent when new features have been added
  19 # because change of extent in provider is not propagated to the layer
  20 vl.updateExtents()
  21 
  22 # show some stats
  23 print "fields:", pr.fieldCount()
  24 print "features:", pr.featureCount()
  25 e = pr.extent()
  26 print "extent:", e.xMin(),e.yMin(),e.xMax(),e.yMax()
  27 
  28 # iterate over features
  29 f = QgsFeature()
  30 pr.select()
  31 while pr.getNextFeature(f):
  32       print "F:",f.featureId(), f.attributeMap(), f.geometry().asPoint()

Memory provider also supports spatial indexing. This means you can call provider's createSpatialIndex() function. Once the spatial index is created (using QgsSpatialIndex class), you'll be able to iterate over features within smaller regions faster (since it's not necessary to traverse all the features, only those in specified rectangle).

Creating bindings

It's good to read the SIP's reference first to get the overall idea: http://www.riverbankcomputing.co.uk/static/Docs/sip4/sipref.html

From looking in the *.sip files you can see that they're just a bit modified header files which follow some rules and some special directives starting with percent sign. There are not many examples available on net, but PyQt4 sip files are a good example.

PythonBindings (last edited 2008-12-19 08:01:00 by MartinDobias)