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, we need to add the following to the top of our 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 it. While there are lots of things we can do with DEM data, let's calculate how often each unique elevation value occurs within the data.

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 there 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 then 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, like 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 includes an important feature. Notice that we use the following:

QgsGeometry(feature.geometry())

We use the preceding line instead of the following:

feature.geometry()

This is to get the geometry object to add to the list. In other words, we had to create a new geometry object based on the feature's existing geometry object. This is 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 see whether the destination geometry is mathematically valid and is not empty.

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 the QgsGeometry objects 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 chapter 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. We'll then save the resulting 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())

Now when our script is run, 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 book, 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 responds when the user clicks on the canvas. In this method, we're going to convert from the pixel coordinates that the user clicked on to the corresponding 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: