Working with geospatial data in the console

So far, we have used the QGIS Console as a glorified Python interpreter, running standard Python programs and manipulating the QGIS user interface. But QGIS is a Geographical Information System (GIS), and one of the main uses of a GIS is to manipulate and query geospatial data. So, let's write some Python code to work with geospatial data directly within the QGIS Console.

In the previous chapter, we loaded three shapefiles into a QGIS project using Python. Here is a typical instruction we used to load a shapefile into a QGIS map layer:

layer = iface.addVectorLayer("/path/to/shapefile.shp", "layer_name", "ogr")

While this is useful if you want to create a QGIS project programmatically, you may just want to load a shapefile so you can analyze its contents, without putting the data into a map layer. To do this, we have to get an appropriate data provider and ask it to open the shapefile, like this:

registry = QgsProviderRegistry.instance()
provider = registry.provider("ogr","/path/to/shapefile.shp")
if not provider.isValid():
    print "Invalid shapefile."
    return

The isValid() method will return False if the shapefile cannot be loaded; this allows us to fail gracefully if there is an error.

Once we have the data provider, we can ask it for the list of fields used to hold the attribute values for each of the shapefile's features:

for field in provider.fields():
      print field.name(), field.typeName()

We can also scan through the features within the shapefile using a QgsFeatureRequest object. For example:

for feature in provider.getFeatures(QgsFeatureRequest()):
    print feature.attribute("name")

Of course, this is just a taste of what can be done using the QGIS libraries to query and manipulate geospatial data. However, let's use what we've learned to build a simple program that calculates and displays information about the contents of a shapefile. Shapefiles hold geospatial features such as polygons, lines and points, and each feature can have any number of attributes associated with it. We'll write a program that opens and scans through a shapefile, identifying the features and calculating the length of each line feature and the area of each polygon feature. We'll also calculate the total length and area across all the features.

One of the challenges we'll have to deal with is the fact that the shapefile can be in any map projection. This means that our calculation of the area and length has to take the map projection into account; if, for example, we simply calculated the linear length of a feature in a shapefile that uses the EPSG 4326 projection (that is, lat/long coordinates), then the calculated length will be in degrees of latitude and longitude—which is a completely meaningless figure. We'll want to calculate the feature lengths in kilometers, and the areas in square kilometers. This is possible but requires us to do a bit more work.

Let's get started with our program. Start by creating a new Python script and enter the following:

from PyQt4.QtGui import *

To make the program easier to use, we're going to define a function and place all our program logic inside this function, like this:

def analyze_shapefile():
    ...

analyze_shapefile()

Now, let's start writing the contents of the analyze_shapefile() function. So far, we've been hardwiring the name of the shapefile, but this time, let's use QGIS's graphical interface to prompt the user to select a shapefile:

def analyze_shapefile():
    filename = QFileDialog.getOpenFileName(iface.mainWindow(),
                                           "Select Shapefile",
                                           "~", '*.shp')
    if not filename:
        print "Cancelled."
        return

We can then open the selected shapefile:

    registry = QgsProviderRegistry.instance()
    provider = registry.provider("ogr",filename)
    if not provider.isValid():
        print "Invalid shapefile."
        return

In order to identify a feature, we need to display a meaningful label for the feature. To do this, we'll look for an attribute with a likely-looking name. If there is no suitable attribute, we'll have to use the feature's ID instead.

Let's start by building a list of the various attributes stored in this shapefile:

    attr_names = []
    for field in provider.fields():
        attr_names.append(field.name())

We're now ready to start scanning through the shapefile's features. Before we do this, though, let's initialize a couple of variables to hold the totals we need to calculate:

    tot_length = 0
    tot_area = 0

We also need to set up a QgsDistanceArea object to do the distance and area calculations for us.

    crs = provider.crs()
    calculator = QgsDistanceArea()
    calculator.setSourceCrs(crs)
    calculator.setEllipsoid(crs.ellipsoidAcronym())
    calculator.setEllipsoidalMode(crs.geographicFlag())

We'll use this object to calculate the true length and area of the shapefile's features in meters and square meters respectively.

We're now ready to scan through the contents of the shapefile, processing each feature in turn:

    for feature in provider.getFeatures(QgsFeatureRequest()):
        ...

For each feature, we want to calculate a label that identifies that feature. We'll do this by looking for an attribute called "name", "NAME", or "Name", and using that attribute's value as the feature label. If there is no attribute with one of these field names, we'll fall back to using the feature's ID instead. Here is the relevant code:

        if "name" in attr_names:
            feature_label = feature.attribute("name")
        elif "Name" in attr_names:
            feature_label = feature.attribute("Name")
        elif "NAME" in attr_names:
            feature_label = feature.attribute("NAME")
        else:
            feature_label = str(feature.id())

Next, we need to obtain the geometry object associated with the feature. The geometry object represents a polygon, line, or point. Getting a reference to the feature's underlying geometry object is simple:

        geometry = feature.geometry()

We can now use the QgsDistanceArea calculator we initialized earlier to calculate the length of a line feature and the area of a polygon feature. To do this, we'll first have to identify the type of feature we are dealing with:

        if geometry.type() == QGis.Line:
            ...
        elif geometry.type() == QGis.Polygon:
            ...
        else:
            ...

For line geometries, we'll calculate the length of the line and update the total length:

        if geometry.type() == QGis.Line:
            length = int(calculator.measure (geometry) / 1000)
            tot_length = tot_length + length
            feature_info = "line of length %d kilometers" % length

For polygon geometries, we'll calculate the area of the polygon and update the total area:

        elif geometry.type() == QGis.Polygon:
            area = int(calculator.measure (geometry) / 1000000)
            tot_area = tot_area + area
            feature_info = "polygon of area %d square kilometers" % area

Finally, for the other types of geometries, we'll simply display the geometry's type:

        else:
            geom_type = qgis.vectorGeometryType(geometry.type())
            feature_info = "geometry of type %s" % geom_type

Now that we've done these calculations, we can display the feature's label together with the information we calculated about this feature:

        print "%s: %s" % (feature_label, feature_info)

Finally, when we've finished iterating over the features, we can display the total line length and polygon area for all the features in the shapefile:

    print "Total length of all line features: %d" % tot_length
    print "Total area of all polygon features: %d" % tot_area

This completes our program for analyzing the contents of a shapefile. The full source for this program is available in the code samples provided with this book. To test our program, type or copy and paste it into the console's script editor, save the file, and click on the Run Script button (or press Ctrl + Shift + E). Here's an example of what the program's output looks like:

Antigua and Barbuda: polygon of area 549 square kilometers
Algeria: polygon of area 2334789 square kilometers
Azerbaijan: polygon of area 86109 square kilometers
Albania: polygon of area 28728 square kilometers
Armenia: polygon of area 29732 square kilometers
...
Jersey: polygon of area 124 square kilometers
South Georgia South Sandwich Islands: polygon of area 3876 square kilometers
Taiwan: polygon of area 36697 square kilometers
Total length of all line features: 0
Total area of all polygon features: 147363163

As you can see, it is quite possible to create Python programs that read and analyze geospatial data, and you can run these programs directly from within the QGIS Console. It is also possible to create and manipulate geospatial data sources using the PyQGIS libraries.