In this article, we will take a closer look at the Python libraries available for the QGIS Python developer, and also look at the various ways in which we can use  these libraries to perform useful tasks within QGIS.

In particular, you will learn:

• How the QGIS Python libraries are based on the underlying C++ APIs
• How to use the C++ API documentation as a reference to work with the Python APIs
• How the PyQGIS libraries are organized
• The most important concepts and classes within the PyQGIS libraries and how to use them
• Some practical examples of performing useful tasks using PyQGIS

About the QGIS Python APIs

The QGIS system itself is written in C++ and has its own set of APIs that are also written in C++. The Python APIs are implemented as wrappers around these C++ APIs. For example, there is a Python class named QgisInterface that acts as a wrapper around a C++ class of the same name. All the methods, class variables, and the like, which are implemented by the C++ version of QgisInterface are made available through the Python wrapper.

What this means is that when you access the Python QGIS APIs, you aren’t accessing the API directly. Instead, the wrapper connects your code to the underlying C++ objects and methods, as follows :

 Fortunately, in most cases, the QGIS Python wrappers simply hide away the complexity of the underlying C++ code, so the PyQGIS libraries work as you would expect them to. There are some gotchas, however, and we will cover these as they come up.

Deciphering the C++ documentation

As QGIS is implemented in C++, the documentation for the QGIS APIs is all based on C++. This can make it difficult for Python developers to understand and work with the QGIS APIs. For example, the API documentation for the QgsInterface.zoomToActiveLayer() method:

 If you’re not familiar with C++, this can be quite confusing. Fortunately, as a Python programmer, you can skip over much of this complexity because it doesn’t apply to you. In particular:

• The virtual keyword is an implementation detail you don’t need to worry about
• The word void indicates that the method doesn’t return a value
• The double colons in QgisInterface::zoomToActiveLayer are simply a C++ convention for separating the class name from the method name

Just like in Python, the parentheses show that the method doesn’t take any parameters. So if you have an instance of QgisInterface (for example, as the standard iface variable available in the Python Console), you can call this method simply by typing the following:

iface.zoomToActiveLayer()

Now, let’s take a look at a slightly more complex example: the C++ documentation for the QgisInterface.addVectorLayer() method looks like the following:

 Notice how the virtual keyword is followed by QgsVectorLayer* instead of void. This is the return value for this method; it returns a QgsVector object.

Technically speaking, * means that the method returns a pointer to an object of type QgsVectorLayer. Fortunately, Python wrappers automatically handle pointers, so you don’t need to worry about this. 

Notice the brief description at the bottom of the documentation for this method; while many of the C++ methods have very little, if any, additional information, other methods have quite extensive descriptions. Obviously, you should read these descriptions carefully as they tell you more about what the method does.

Even without any description, the C++ documentation is still useful as it tells you what the method is called, what parameters it accepts, and what type of data is being returned.

In the preceding method, you can see that there are three parameters listed in between the parentheses. As C++ is a strongly typed language, you have to define the type of each parameter when you define a function. This is helpful for Python programmers as it tells you what type of value to supply. Apart from QGIS objects, you might also encounter the following data types in the C++ documentation:

Data type	Description
int	A standard Python integer value
long	A standard Python long integer value
float	A standard Python floating point (real) number
bool	A Boolean value (true or false)
QString	A string value. Note that the QGIS Python wrappers automatically convert Python strings to C++ strings, so you don’t need to deal with QString objects directly
QList	This object is used to encapsulate a list of other objects. For example, QList<QString*> represents a list of strings

Just as in Python, a method can take default values for each parameter. For example, the QgisInterface.newProject() method looks like the following:

 In this case, the thePromptToSaveFlag parameter has a default value, and this default value will be used if no value is supplied.

In Python, classes are initialized using the __init__ method. In C++, this is called a constructor. For example, the constructor for the QgsLabel class looks like the following:

 Just as in Python, C++ classes inherit the methods defined in their superclass. Fortunately, QGIS doesn’t have an extensive class hierarchy, so most of the classes don’t have a superclass. However, don’t forget to check for a superclass if you can’t find the method you’re looking for in the documentation for the class itself.

Finally, be aware that C++ supports the concept of method overloading. A single method can be defined more than once, where each version accepts a different set of parameters. For example, take a look at the constructor for the QgsRectangle class—you will see that there are four different versions of this method.

The first version accepts the four coordinates as floating point numbers:

 The second version constructs a rectangle using two QgsPoint objects:

 The third version copies the coordinates from QRectF (which is a Qt data type) into a QgsRectangle object:

 The final version copies the coordinates from another QgsRectangle object:

 The C++ compiler chooses the correct method to use based on the parameters that have been supplied. Python has no concept of method overloading; just choose the version of the method that accepts the parameters you want to supply, and the QGIS Python wrappers will automatically choose the correct method for you.

If you keep these guidelines in mind, deciphering the C++ documentation for QGIS isn’t all that hard. It just looks more complicated than it really is, thanks to all the complexity specific to C++. However, it won’t take long for your brain to start filtering out C++ and use the QGIS reference documentation almost as easily as if it was written for Python rather than C++.

Organization of the QGIS Python libraries

Now that we can understand the C++-oriented documentation, let’s see how the PyQGIS libraries are structured. All of the PyQGIS libraries are organized under a package named qgis. You wouldn’t normally import qgis directly, however, as all the interesting libraries are subpackages within this main package; here are the five packages that make up the PyQGIS library:

qgis.core	This provides access to the core GIS functionality used throughout QGIS.
qgis.gui	This defines a range of GUI widgets that you can include in your own programs.
qgis.analysis	This provides spatial analysis tools to analyze vector and raster format data.
qgis.networkanalysis	This provides tools to build and analyze topologies.
qgis.utils	This implements miscellaneous functions that allow you to work with the QGIS application using Python.

 The first two packages (qgis.core and qgis.gui) implement the most important parts of the PyQGIS library, and it’s worth spending some time to become more familiar with the concepts and classes they define. Now let’s take a closer look at these two packages.

The qgis.core package

The qgis.core package defines fundamental classes used throughout the QGIS system. A large part of this package is dedicated to working with vector and raster format geospatial data, and displaying these types of data within a map. Let’s take
a look at how this is done.

Maps and map layers

A map consists of multiple layers drawn one on top of the other:

 There are three types of map layers supported by QGIS:

• Vector layer: This layer draws geospatial features such as points, lines, and polygons
• Raster layer: This layer draws raster (bitmapped) data onto a map
• Plugin layer: This layer allows a plugin to draw directly onto a map

Each of these types of map layers has a corresponding class within the qgis.core library. For example, a vector map layer will be represented by an object of type qgis.core.QgsVectorLayer.

We will take a closer look at vector and raster map layers shortly. Before we do this, though, we need to learn how geospatial data (both vector and raster data)
is positioned on a map.

Coordinate reference systems

Since the Earth is a three-dimensional object, maps will only represent the Earth’s surface as a two-dimensional plane, so there has to be a way of translating points on the Earth’s surface into (x,y) coordinates within a map. This is done using a Coordinate Reference System (CRS):

 Globe image courtesy Wikimedia (http://commons.wikimedia.org/wiki/File:Rotating_globe.gif)

A CRS has two parts: an ellipsoid, which is a mathematical model of the Earth’s surface, and a projection, which is a formula that converts points on the surface of the spheroid into (x,y) coordinates on a map.

Generally you won’t need to worry about all these details. You can simply select the appropriate CRS that matches the CRS of the data you are using. However, as there are many different coordinate reference systems that have been devised over the years, it is vital that you use the correct CRS when plotting your geospatial data. If you don’t do this, your features will be displayed in the wrong place or have the wrong shape.

The majority of geospatial data available today uses the EPSG 4326 coordinate reference system (sometimes also referred to as WGS84). This CRS defines coordinates as latitude and longitude values. This is the default CRS used for new data imported into QGIS. However, if your data uses a different coordinate reference system, you will need to create and use a different CRS for your map layer.

The qgis.core.QgsCoordinateReferenceSystem class represents a CRS. Once you create your coordinate reference system, you can tell your map layer to use that CRS when accessing the underlying data. For example:

crs = QgsCoordinateReferenceSystem(4326,
           QgsCoordinateReferenceSystem.EpsgCrsId))
layer.setCrs(crs)

Note that different map layers can use different coordinate reference systems. Each layer will use its CRS when drawing the contents of the layer onto the map.

Vector layers

A vector layer draws geospatial data onto a map in the form of points, lines, polygons, and so on. Vector-format geospatial data is typically loaded from a vector data source such as a shapefile or database. Other vector data sources can hold vector data in memory, or load data from a web service across the internet.

A vector-format data source has a number of features, where each feature represents a single record within the data source. The qgis.core.QgsFeature class represents a feature within a data source. Each feature has the following principles:

• ID: This is the feature’s unique identifier within the data source.
• Geometry: This is the underlying point, line, polygon, and so on, which represents the feature on the map. For example, a data source representing cities would have one feature for each city, and the geometry would typically be either a point that represents the center of the city, or a polygon (or a multipolygon) that represents the city’s outline.
• Attributes: These are key/value pairs that provide additional information about the feature. For example, a city data source might have attributes such as total_area, population, elevation, and so on. Attribute values can be strings, integers, or floating point numbers.

In QGIS, a data provider allows the vector layer to access the features within the data source. The data provider, an instance of qgis.core.QgsVectorDataProvider, includes:

• A geometry type that is stored in the data source
• A list of fields that provide information about the attributes stored for each feature
• The ability to search through the features within the data source, using the getFeatures() method and the QgsFeatureRequest class

You can access the various vector (and also raster) data providers by using the qgis.core.QgsProviderRegistry class.

The vector layer itself is represented by a qgis.core.QgsVectorLayer object. Each vector layer includes:

• Data provider: This is the connection to the underlying file or database
  that holds the geospatial information to be displayed
• Coordinate reference system: This indicates which CRS the geospatial
  data uses
• Renderer: This chooses how the features are to be displayed

Let’s take a closer look at the concept of a renderer and how features are displayed within a vector map layer.

Displaying vector data

The features within a vector map layer are displayed using a combination of renderer and symbol objects. The renderer chooses the symbol to use for a given feature, and the symbol that does the actual drawing.

There are three basic types of symbols defined by QGIS:

• Marker symbol: This displays a point as a filled circle
• Line symbol: This draws a line using a given line width and color
• Fill symbol: This draws the interior of a polygon with a given color

These three types of symbols are implemented as subclasses of the qgis.core.QgsSymbolV2 class:

• qgis.core.QgsMarkerSymbolV2
• qgis.core.QgsLineSymbolV2
• qgis.core.QgsFillSymbolV2

Internally, symbols are rather complex, as “symbol layers” allow multiple elements to be drawn on top of each other. In most cases, however, you can make use of the “simple” version of the symbol. This makes it easier to create a new symbol without having to deal with the internal complexity of symbol layers. For example:

symbol = QgsMarkerSymbolV2.createSimple({'width' : 1.0, 'color' : "255,0,0"})

While symbols draw the features onto the map, a renderer is used to choose which symbol to use to draw a particular feature. In the simplest case, the same symbol is used for every feature within a layer. This is called a single symbol renderer, and is represented by the qgis.core.QgsSingleSymbolRenderV2 class. Other possibilities include:

• Categorized symbol renderer (qgis.core.QgsCategorizedSymbolRendererV2): This renderer chooses a symbol based on the value of an attribute. The categorized symbol renderer has a mapping from attribute values to symbols.
• Graduated symbol renderer (qgis.core.QgsGraduatedSymbolRendererV2): This type of renderer has a range of attribute, values, and maps each range to an appropriate symbol.

Using a single symbol renderer is very straightforward:

symbol = ...
renderer = QgsSingleSymbolRendererV2(symbol)
layer.setRendererV2(renderer)

To use a categorized symbol renderer, you first define a list of qgis.core. QgsRendererCategoryV2 objects, and then use that to create the renderer.

For example:

symbol_male = ...
symbol_female = ...
 
categories = []
categories.append(QgsRendererCategoryV2("M", symbol_male, "Male"))
categories.append(QgsRendererCategoryV2("F", symbol_female, "Female"))
renderer = QgsCategorizedSymbolRendererV2("", categories)
renderer.setClassAttribute("GENDER")
layer.setRendererV2(renderer)

Notice that the QgsRendererCategoryV2 constructor takes three parameters: the desired value, the symbol used, and a label used to describe that category.

Finally, to use a graduated symbol renderer, you define a list of qgis.core.QgsRendererRangeV2 objects and then use that to create your renderer. For example:

symbol1 = ...
symbol2 = ...
 
ranges = []
ranges.append(QgsRendererRangeV2(0, 10, symbol1, "Range 1"))
ranges.append(QgsRendererRange(11, 20, symbol2, "Range 2"))
 
renderer = QgsGraduatedSymbolRendererV2("", ranges)
renderer.setClassAttribute("FIELD")
layer.setRendererV2(renderer)

Accessing vector data

In addition to displaying the contents of a vector layer within a map, you can use Python to directly access the underlying data. This can be done using the data provider’s getFeatures() method. For example, to iterate over all the features within the layer, you can do the following:

provider = layer.dataProvider()
for feature in provider.getFeatures(QgsFeatureRequest()):
  ...

If you want to search for features based on some criteria, you can use the QgsFeatureRequest object’s setFilterExpression() method, as follows:

provider = layer.dataProvider()
request = QgsFeatureRequest()
request.setFilterExpression('"GENDER" = "M"')
for feature in provider.getFeatures(QgsFeatureRequest()):
  ...

Once you have the features, it’s easy to get access to the feature’s geometry, ID,
and attributes. For example:

geometry = feature.geometry()
  id = feature.id()
  name = feature.attribute("NAME")

The object returned by the feature.geometry() call, which will be an instance of qgis.core.QgsGeometry, represents the feature’s geometry. This object has a number of methods you can use to extract the underlying data and perform various geospatial calculations.

Spatial indexes

In the previous section, we searched for features based on their attribute values. There are times, though, when you might want to find features based on their position in space. For example, you might want to find all features that lie within a certain distance of a given point. To do this, you can use a spatial index, which indexes features according to their location and extent. Spatial indexes are represented in QGIS by the QgsSpatialIndex class.

For performance reasons, a spatial index is not created automatically for each vector layer. However, it’s easy to create one when you need it:

provider = layer.dataProvider()
index = QgsSpatialIndex()
for feature in provider.getFeatures(QgsFeatureRequest()):
  index.insertFeature(feature)

Don’t forget that you can use the QgsFeatureRequest.setFilterExpression() method to limit the set of features that get added to the index.

Once you have the spatial index, you can use it to perform queries based on the position of the features. In particular:

• You can find one or more features that are closest to a given point using the nearestNeighbor() method. For example:

features = index.nearestNeighbor(QgsPoint(long, lat), 5)

Note that this method takes two parameters: the desired point as a QgsPoint object and the number of features to return.

• You can find all features that intersect with a given rectangular area by using the intersects() method, as follows:

features = index.intersects(QgsRectangle(left, bottom, right, top))

Raster layers

Raster-format geospatial data is essentially a bitmapped image, where each pixel or “cell” in the image corresponds to a particular part of the Earth’s surface. Raster data is often organized into bands, where each band represents a different piece of information. A common use for bands is to store the red, green, and blue component of the pixel’s color in a separate band. Bands might also represent other types of information, such as moisture level, elevation, or soil type.

There are many ways in which raster information can be displayed. For example:

• If the raster data only has one band, the pixel value can be used as an index into a palette. The palette maps each pixel value maps to a particular color.
• If the raster data has only one band but no palette is provided. The pixel values can be used directly as a grayscale value; that is, larger numbers are lighter and smaller numbers are darker. Alternatively, the pixel values can be passed through a pseudocolor algorithm to calculate the color to be displayed.
• If the raster data has multiple bands, then typically, the bands would be combined to generate the desired color. For example, one band might represent the red component of the color, another band might represent the green component, and yet another band might represent the blue component.
• Alternatively, a multiband raster data source might be drawn using a palette, or as a grayscale or a pseudocolor image, by selecting a particular band to use for the color calculation.

Let’s take a closer look at how raster data can be drawn onto the map.

Displaying raster data

The drawing style associated with the raster band controls how the raster data will be displayed. The following drawing styles are currently supported:

Drawing style	Description
PalettedColor	For a single band raster data source, a palette maps each raster value to a color.
SingleBandGray	For a single band raster data source, the raster value is used directly as a grayscale value.
SingleBandPseudoColor	For a single band raster data source, the raster value is used to calculate a pseudocolor.
PalettedSingleBandGray	For a single band raster data source that has a palette, this drawing style tells QGIS to ignore the palette and use the raster value directly as a grayscale value.
PalettedSingleBandPseudoColor	For a single band raster data source that has a palette, this drawing style tells QGIS to ignore the palette and use the raster value to calculate a pseudocolor.
MultiBandColor	For multiband raster data sources, use a separate band for each of the red, green, and blue color components. For this drawing style, the setRedBand(), setGreenBand(), and setBlueBand() methods can be used to choose which band to use for each color component.
MultiBandSingleBandGray	For multiband raster data sources, choose a single band to use as the grayscale color value. For this drawing style, use the setGrayBand() method to specify the band to use.
MultiBandSingleBandPseudoColor	For multiband raster data sources, choose a single band to use to calculate a pseudocolor. For this drawing style, use the setGrayBand() method to specify the band to use.

 To set the drawing style, use the layer.setDrawingStyle() method, passing in a string that contains the name of the desired drawing style. You will also need to call the various setXXXBand() methods, as described in the preceding table, to tell the raster layer which bands contain the value(s) to use to draw each pixel.

Note that QGIS doesn’t automatically update the map when you call the preceding functions to change the way the raster data is displayed. To have your changes displayed right away, you’ll need to do the following:

1. Turn off raster image caching. This can be done by calling layer.setImageCache(None).
2. Tell the raster layer to redraw itself, by calling layer.triggerRepaint().

Accessing raster data

As with vector-format data, you can access the underlying raster data via the data provider’s identify() method. The easiest way to do this is to pass in a single coordinate and retrieve the value or values of that coordinate. For example:

provider = layer.dataProvider()
values = provider.identify(QgsPoint(x, y),
              QgsRaster.IdentifyFormatValue)
if values.isValid():
  for band,value in values.results().items():
    ...

As you can see, you need to check whether the given coordinate exists within the raster data (using the isValid() call). The values.results() method returns a dictionary that maps band numbers to values.

Using this technique, you can extract all the underlying data associated with a given coordinate within the raster layer.

You can also use the provider.block() method to retrieve the band data for a large number of coordinates all at once. We will look at how to do this later in this article.

Other useful qgis.core classes

Apart from all the classes and functionality involved in working with data sources and map layers, the qgis.core library also defines a number of other classes that you might find useful:

Class	Description
QgsProject	This represents the current QGIS project. Note that this is a singleton object, as only one project can be open at a time. The QgsProject class is responsible for loading and storing properties, which can be useful for plugins.
QGis	This class defines various constants, data types, and functions used throughout the QGIS system.
QgsPoint	This is a generic class that stores the coordinates for a point within a two-dimensional plane.
QgsRectangle	This is a generic class that stores the coordinates for a rectangular area within a two-dimensional plane.
QgsRasterInterface	This is the base class to use for processing raster data. This can be used, to reproject a set of raster data into a new coordinate system, to apply filters to change the brightness or color of your raster data, to resample the raster data, and to generate new raster data by rendering the existing data in various ways.
QgsDistanceArea	This class can be used to calculate distances and areas for a given geometry, automatically converting from the source coordinate reference system into meters.
QgsMapLayerRegistry	This class provides access to all the registered map layers in the current project.
QgsMessageLog	This class provides general logging features within a QGIS program. This lets you send debugging messages, warnings, and errors to the QGIS “Log Messages” panel.

 The qgis.gui package

The qgis.gui package defines a number of user-interface widgets that you can include in your programs. Let’s start by looking at the most important qgis.gui classes, and follow this up with a brief look at some of the other classes that you might find useful.

The QgisInterface class

QgisInterface represents the QGIS system’s user interface. It allows programmatic access to the map canvas, the menu bar, and other parts of the QGIS application. When running Python code within a script or a plugin, or directly from the QGIS Python console, a reference to QgisInterface is typically available through the iface global variable.

The QgisInterface object is only available when running the QGIS application itself. If you are running an external application and import the PyQGIS library into your application, QgisInterface won’t be available.

Some of the more important things you can do with the QgisInterface object are:

• Get a reference to the list of layers within the current QGIS project via the legendInterface() method.
• Get a reference to the map canvas displayed within the main application window, using the mapCanvas() method.
• Retrieve the currently active layer within the project, using the activeLayer() method, and set the currently active layer by using the setActiveLayer() method.
• Get a reference to the application’s main window by calling the mainWindow() method. This can be useful if you want to create additional
  Qt windows or dialogs that use the main window as their parent.
• Get a reference to the QGIS system’s message bar by calling the messageBar() method. This allows you to display messages to the user directly within the QGIS main window.

The QgsMapCanvas class

The map canvas is responsible for drawing the various map layers into a window. The QgsMapCanvas class represents a map canvas. This class includes:

• A list of the currently shown map layers. This can be accessed using the layers() method.

  Note that there is a subtle difference between the list of map layers available within the map canvas and the list of map layers included in the QgisInterface.legendInterface() method. The map canvas’s list of layers only includes the list of layers currently visible, while QgisInterface.legendInterface() returns all the map layers, including those that are currently hidden. 

• The map units used by this map (meters, feet, degrees, and so on). The map’s units can be retrieved by calling the mapUnits() method.
• An extent,which is the area of the map that is currently shown within the canvas. The map’s extent will change as the user zooms in and out, and pans across the map. The current map extent can be obtained by calling the extent() method.
• A current map tool that controls the user’s interaction with the contents of the map canvas. The current map tool can be set using the setMapTool() method, and you can retrieve the current map tool (if any) by calling the mapTool() method.
• A background color used to draw the background behind all the map layers. You can change the map’s background color by calling the canvasColor() method.
• A coordinate transform that converts from map coordinates (that is, coordinates in the data source’s coordinate reference system) to pixels within the window. You can retrieve the current coordinate transform by calling the getCoordinateTransform() method.

The QgsMapCanvasItem class

A map canvas item is an item drawn on top of the map canvas. The map canvas item will appear in front of the map layers. While you can create your own subclass of QgsMapCanvasItem if you want to draw custom items on top of the map canvas, it would be more useful for you to make use of an existing subclass that will do the work for you. There are currently three subclasses of QgsMapCanvasItem that you might find useful:

• QgsVertexMarker: This draws an icon (an “X”, a “+”, or a small box) centered around a given point on the map.
• QgsRubberBand: This draws an arbitrary polygon or polyline onto the map. It is intended to provide visual feedback as the user draws a polygon onto the map.
• QgsAnnotationItem: This is used to display additional information about a feature, in the form of a balloon that is connected to the feature. The QgsAnnotationItem class has various subclasses that allow you to customize the way the information is displayed.

The QgsMapTool class

A map tool allows the user to interact with and manipulate the map canvas, capturing mouse events and responding appropriately. A number of QgsMapTool subclasses provide standard map interaction behavior such as clicking to zoom in, dragging to pan the map, and clicking on a feature to identify it. You can also create your own custom map tools by subclassing QgsMapTool and implementing the various methods that respond to user-interface events such as pressing down the mouse button, dragging the canvas, and so on.

Once you have created a map tool, you can allow the user to activate it by associating the map tool with a toolbar button. Alternatively, you can activate it from within your Python code by calling the mapCanvas.setMapTool(…) method.

We will look at the process of creating your own custom map tool in the section Using the PyQGIS library in the following table:

Other useful qgis.gui classes

While the qgis.gui package defines a large number of classes, the ones you are most likely to find useful are given in the following table:

Class	Description
QgsLegendInterface	This provides access to the map legend, that is, the list of map layers within the current project. Note that map layers can be grouped, hidden, and shown within the map legend.
QgsMapTip	This displays a tip on a map canvas when the user holds the mouse over a feature. The map tip will show the display field for the feature; you can set this by calling layer.setDisplayField(“FIELD”).
QgsColorDialog	This is a dialog box that allows the user to select a color.
QgsDialog	This is a generic dialog with a vertical box layout and a button box, making it easy to add content and standard buttons to your dialog.
QgsMessageBar	This is a user interface widget for displaying non-blocking messages to the user. We looked at the message bar class in the previous article.
QgsMessageViewer	This is a generic class that displays long messages to the user within a modal dialog.
QgsBlendModeComboBox QgsBrushStyleComboBox QgsColorRampComboBox QgsPenCapStyleComboBox QgsPenJoinStyleComboBox QgsScaleComboBox	These QComboBox user-interface widgets allow you to prompt the user for various drawing options. With the exception of the QgsScaleComboBox, which lets the user choose a map scale, all the other QComboBox subclasses let the user choose various Qt drawing options.

 Using the PyQGIS library

In the previous section, we looked at a number of classes provided by the PyQGIS library. Let’s make use of these classes to perform some real-world geospatial development tasks.

Analyzing raster data

We’re going to start by writing a program to load in some raster-format data and analyze its contents. To make this more interesting, we’ll use a Digital Elevation Model (DEM) file, which is a raster format data file that contains elevation data.

The Global Land One-Kilometer Base Elevation Project (GLOBE) provides free DEM data for the world, where each pixel represents one square kilometer of the Earth’s surface. GLOBE data can be downloaded from http://www.ngdc.noaa.gov/mgg/topo/gltiles.html. Download the E tile, which includes the western half of the USA. The resulting file, which is named e10g, contains the height information you need. You’ll also need to download the e10g.hdr header file so that QGIS can read the file—you can download this from http://www.ngdc.noaa.gov/mgg/topo/elev/esri/hdr. Once you’ve downloaded these two files, put them together into a convenient directory.

You can now load the DEM data into QGIS using the following code:

registry = QgsProviderRegistry.instance()
provider = registry.provider("gdal", "/path/to/e10g")

Unfortunately, there is a slight complexity here. Since QGIS doesn’t know which coordinate reference system is used for the data, it displays a dialog box that asks you to choose the CRS. Since the GLOBE DEM data is in the WGS84 CRS, which QGIS uses by default, this dialog box is redundant. To disable it, you need to add the following to the top of your program:

from PyQt4.QtCore import QSettings
QSettings().setValue("/Projections/defaultBehaviour", "useGlobal")

Now that we’ve loaded our raster DEM data into QGIS, we can analyze. There are lots of things we can do with DEM data, so let’s calculate how often each unique elevation value occurs within the data.

Notice that we’re loading the DEM data directly using QgsRasterDataProvider. We don’t want to display this information on a map, so we don’t want (or need) to load it into QgsRasterLayer. 

Since the DEM data is in a raster format, you need to iterate over the individual pixels or cells to get each height value. The provider.xSize() and provider.ySize() methods tell us how many cells are in the DEM, while the provider.extent() method gives us the area of the Earth’s surface covered by the DEM.

Using this information, we can extract the individual elevation values from the contents of the DEM in the following way:

raster_extent = provider.extent()
raster_width = provider.xSize()
raster_height = provider.ySize()
block = provider.block(1, raster_extent, raster_width, raster_height)

The returned block variable is an object of type QgsRasterBlock, which is essentially a two-dimensional array of values. Let’s iterate over the raster and extract the individual elevation values:

for x in range(raster_width):
  for y in range(raster_height):
    elevation = block.value(x, y)
    ....

Now that we’ve loaded the individual elevation values, it’s easy to build a histogram out of those values. Here is the entire program to load the DEM data into memory, and calculate, and display the histogram:

from PyQt4.QtCore import QSettings
QSettings().setValue("/Projections/defaultBehaviour", "useGlobal")

registry = QgsProviderRegistry.instance()
provider = registry.provider("gdal", "/path/to/e10g")
 
raster_extent = provider.extent()
raster_width = provider.xSize()
raster_height = provider.ySize()
no_data_value = provider.srcNoDataValue(1)
 
histogram = {} # Maps elevation to number of occurrences.
 
block = provider.block(1, raster_extent, raster_width,
            raster_height)
if block.isValid():
  for x in range(raster_width):
    for y in range(raster_height):
      elevation = block.value(x, y)
      if elevation != no_data_value:
        try:
          histogram[elevation] += 1
        except KeyError:
          histogram[elevation] = 1
 
for height in sorted(histogram.keys()):
  print height, histogram[height]

Note that we’ve added a no data value check to the code. Raster data often includes pixels that have no value associated with them. In the case of a DEM, elevation data is only provided for areas of land; pixels over the sea have no elevation, and we have to exclude them, or our histogram will be inaccurate.

Manipulating vector data and saving it to a shapefile

Let’s create a program that takes two vector data sources, subtracts one set of vectors from the other, and saves the resulting geometries into a new shapefile. Along the way, we’ll learn a few important things about the PyQGIS library.

We’ll be making use of the QgsGeometry.difference() function. This function performs a geometrical subtraction of one geometry from another, similar to this:

 Let’s start by asking the user to select the first shapefile and open up a vector data provider for that file:

filename_1 = QFileDialog.getOpenFileName(iface.mainWindow(),
                     "First Shapefile",
                     "~", "*.shp")
if not filename_1:
  return
 
registry = QgsProviderRegistry.instance()
provider_1 = registry.provider("ogr", filename_1)

We can then read the geometries from that file into memory:

geometries_1 = []
for feature in provider_1.getFeatures(QgsFeatureRequest()):
  geometries_1.append(QgsGeometry(feature.geometry()))

This last line of code does something very important that may not be obvious at first. Notice that we use the following:

QgsGeometry(feature.geometry())

We use the preceding line instead of the following:

feature.geometry()

This creates a new instance of the QgsGeometry object, copying the geometry into a new object, rather than just adding the existing geometry object to the list. We have to do this because of a limitation of the way the QGIS Python wrappers work: the feature.geometry() method returns a reference to the geometry, but the C++ code doesn’t know that you are storing this reference away in your Python code. So, when the feature is no longer needed, the memory used by the feature’s geometry is also released. If you then try to access that geometry later on, the entire QGIS system will crash. To get around this, we make a copy of the geometry so that we can refer to it even after the feature’s memory has been released.

Now that we’ve loaded our first set of geometries into memory, let’s do the same for the second shapefile:

filename_2 = QFileDialog.getOpenFileName(iface.mainWindow(),
                     "Second Shapefile",
                     "~", "*.shp")
if not filename_2:
  return
 
provider_2 = registry.provider("ogr", filename_2)
 
geometries_2 = []
for feature in provider_2.getFeatures(QgsFeatureRequest()):
  geometries_2.append(QgsGeometry(feature.geometry()))

With the two sets of geometries loaded into memory, we’re ready to start subtracting one from the other. However, to make this process more efficient, we will combine the geometries from the second shapefile into one large geometry, which we can then subtract all at once, rather than subtracting one at a time. This will make the subtraction process much faster:

combined_geometry = None
for geometry in geometries_2:
  if combined_geometry == None:
    combined_geometry = geometry
  else:
    combined_geometry = combined_geometry.combine(geometry)

We can now calculate the new set of geometries by subtracting one from the other:

dst_geometries = []
for geometry in geometries_1:
  dst_geometry = geometry.difference(combined_geometry)
  if not dst_geometry.isGeosValid(): continue
  if dst_geometry.isGeosEmpty(): continue
  dst_geometries.append(dst_geometry)

Notice that we check to ensure that the destination geometry is mathematically
valid and isn’t empty.

Invalid geometries are a common problem when manipulating complex shapes. There are options for fixing them, such as splitting apart multi-geometries and performing a buffer operation. 

Our last task is to save the resulting geometries into a new shapefile. We’ll first ask the user for the name of the destination shapefile:

dst_filename = QFileDialog.getSaveFileName(iface.mainWindow(),
                      "Save results to:",
                      "~", "*.shp")
if not dst_filename:
  return

We’ll make use of a vector file writer to save the geometries into a shapefile. Let’s start by initializing the file writer object:

fields = QgsFields()
writer = QgsVectorFileWriter(dst_filename, "ASCII", fields,
               dst_geometries[0].wkbType(),
               None, "ESRI Shapefile")
if writer.hasError() != QgsVectorFileWriter.NoError:
  print "Error!"
  return

We don’t have any attributes in our shapefile, so the fields list is empty. Now that the writer has been set up, we can save the geometries into the file:

for geometry in dst_geometries:
  feature = QgsFeature()
  feature.setGeometry(geometry)
  writer.addFeature(feature)

Now that all the data has been written to the disk, let’s display a message box that informs the user that we’ve finished:

QMessageBox.information(iface.mainWindow(), "",
            "Subtracted features saved to disk.")

As you can see, creating a new shapefile is very straightforward in PyQGIS, and it’s easy to manipulate geometries using Python—just so long as you copy QgsGeometry you want to keep around. If your Python code starts to crash while manipulating geometries, this is probably the first thing you should look for.

Using different symbols for different features within a map

Let’s use World Borders Dataset that you downloaded in the previous article to draw a world map, using different symbols for different continents. This is a good example of using a categorized symbol renderer, though we’ll combine it into a script that loads the shapefile into a map layer, and sets up the symbols and map renderer to display the map exactly as you want it. We’ll then save this map as an image.

Let’s start by creating a map layer to display the contents of the World Borders Dataset shapefile:

layer = iface.addVectorLayer("/path/to/TM_WORLD_BORDERS-0.3.shp",
               "continents", "ogr")

Each unique region code in the World Borders Dataset shapefile corresponds to a continent. We want to define the name and color to use for each of these regions, and use this information to set up the various categories to use when displaying the map:

from PyQt4.QtGui import QColor
categories = []
for value,color,label in [(0,   "#660000", "Antarctica"),
                          (2,   "#006600", "Africa"),
                          (9,   "#000066", "Oceania"),
                          (19,  "#660066", "The Americas"),
                          (142, "#666600", "Asia"),
                          (150, "#006666", "Europe")]:
  symbol = QgsSymbolV2.defaultSymbol(layer.geometryType())
  symbol.setColor(QColor(color))
  categories.append(QgsRendererCategoryV2(value, symbol, label))

With these categories set up, we simply update the map layer to use a categorized renderer based on the value of the region attribute, and then redraw the map:

layer.setRendererV2(QgsCategorizedSymbolRendererV2("region",
                          categories))
layer.triggerRepaint()

There’s only one more thing to do; since this is a script that can be run multiple times, let’s have our script automatically remove the existing continents layer, if it exists, before adding a new one. To do this, we can add the following to the start of our script:

layer_registry = QgsMapLayerRegistry.instance()

for layer in layer_registry.mapLayersByName(“continents”):

  layer_registry.removeMapLayer(layer.id())

When our script is running, it will create one (and only one) layer that shows the various continents in different colors. These will appear as different shades of gray
in the printed article, but the colors will be visible on the computer screen:

Now, let’s use the same data set to color each country based on its relative population. We’ll start by removing the existing population layer, if it exists:

layer_registry = QgsMapLayerRegistry.instance()
for layer in layer_registry.mapLayersByName("population"):
  layer_registry.removeMapLayer(layer.id())

Next, we open the World Borders Dataset into a new layer called “population”:

layer = iface.addVectorLayer("/path/to/TM_WORLD_BORDERS-0.3.shp",
               "population", "ogr")

We then need to set up our various population ranges:

from PyQt4.QtGui import QColor
ranges = []
for min_pop,max_pop,color in [(0,        99999,     "#332828"),
                              (100000,   999999,    "#4c3535"),
                              (1000000,  4999999,   "#663d3d"),
                              (5000000,  9999999,   "#804040"),
                              (10000000, 19999999,  "#993d3d"),
                              (20000000, 49999999,  "#b33535"),
                              (50000000, 999999999, "#cc2828")]:
  symbol = QgsSymbolV2.defaultSymbol(layer.geometryType())
  symbol.setColor(QColor(color))
  ranges.append(QgsRendererRangeV2(min_pop, max_pop,
                   symbol, ""))

Now that we have our population ranges and their associated colors, we simply set up a graduated symbol renderer to choose a symbol based on the value of the pop2005 attribute, and tell the map to redraw itself:

layer.setRendererV2(QgsGraduatedSymbolRendererV2("pop2005",
                         ranges))
layer.triggerRepaint()

The result will be a map layer that shades each country according to its population:

 Calculating the distance between two user-defined points 

In our final example of using the PyQGIS library, we’ll write some code that, when run, starts listening for mouse events from the user. If the user clicks on a point, drags the mouse, and then releases the mouse button again, we will display the distance between those two points. This is a good example of how to add your 

own map interaction logic to QGIS, using the QgsMapTool class.

This is the basic structure for our QgsMapTool subclass:

class DistanceCalculator(QgsMapTool):
  def __init__(self, iface):
    QgsMapTool.__init__(self, iface.mapCanvas())
    self.iface = iface
 
  def canvasPressEvent(self, event):
    ...
 
  def canvasReleaseEvent(self, event):
    ...

To make this map tool active, we’ll create a new instance of it and pass it to the mapCanvas.setMapTool() method. Once this is done, our canvasPressEvent() and canvasReleaseEvent() methods will be called whenever the user clicks or releases the mouse button over the map canvas.

Let’s start with the code that handles the user clicking on the canvas. In this method, we’re going to convert from the pixel coordinates that the user clicked on to the map coordinates (that is, a latitude and longitude value). We’ll then remember these coordinates so that we can refer to them later. Here is the necessary code:

def canvasPressEvent(self, event):
  transform = self.iface.mapCanvas().getCoordinateTransform()
  self._startPt = transform.toMapCoordinates(event.pos().x(),
                        event.pos().y())

When the canvasReleaseEvent() method is called, we’ll want to do the same with the point at which the user released the mouse button:

def canvasReleaseEvent(self, event):
  transform = self.iface.mapCanvas().getCoordinateTransform()
  endPt = transform.toMapCoordinates(event.pos().x(),
                    event.pos().y())

Now that we have the two desired coordinates, we’ll want to calculate the distance between them. We can do this using a QgsDistanceArea object:

     crs = self.iface.mapCanvas().mapRenderer().destinationCrs()
  distance_calc = QgsDistanceArea()
  distance_calc.setSourceCrs(crs)
  distance_calc.setEllipsoid(crs.ellipsoidAcronym())
  distance_calc.setEllipsoidalMode(crs.geographicFlag())
  distance = distance_calc.measureLine([self._startPt,
                     endPt]) / 1000

Notice that we divide the resulting value by 1000. This is because the QgsDistanceArea object returns the distance in meters, and we want to display the distance in kilometers.

Finally, we’ll display the calculated distance in the QGIS message bar:

  messageBar = self.iface.messageBar()
  messageBar.pushMessage("Distance = %d km" % distance,
              level=QgsMessageBar.INFO,
              duration=2)

Now that we’ve created our map tool, we need to activate it. We can do this by adding the following to the end of our script:

calculator = DistanceCalculator(iface)
iface.mapCanvas().setMapTool(calculator)

With the map tool activated, the user can click and drag on the map. When the mouse button is released, the distance (in kilometers) between the two points will be displayed in the message bar:

Summary

In this article, we took an in-depth look at the PyQGIS libraries and how you can use them in your own programs. We learned that the QGIS Python libraries are implemented as wrappers around the QGIS APIs implemented in C++. We saw how Python programmers can understand and work with the QGIS reference documentation, even though it is written for C++ developers. We also looked at the way the PyQGIS libraries are organized into different packages, and learned about the most important classes defined in the qgis.core and qgis.gui packages.

We then saw how a coordinate reference systems (CRS) is used to translate from points on the three-dimensional surface of the Earth to coordinates within a two-dimensional map plane.

We learned that vector format data is made up of features, where each feature has an ID, a geometry, and a set of attributes, and that symbols are used to draw vector geometries onto a map layer, while renderers are used to choose which symbol to use for a given feature.

We learned how a spatial index can be used to speed up access to vector features.

Next, we saw how raster format data is organized into bands that represent information such as color, elevation, and so on, and looked at the various ways in which a raster data source can be displayed within a map layer. Along the way, we learned how to access the contents of a raster data source.

Finally, we looked at various techniques for performing useful tasks using the PyQGIS library.

In the next article, we will learn more about QGIS Python plugins, and then go on to use the plugin architecture as a way of implementing a useful feature within a mapping application.

Resources for Article:

 

---

Further resources on this subject:

• QGIS Feature Selection Tools [article]
• Server Logs [article]

---