Table of Contents for
QGIS: Becoming a GIS Power User

Version ebook / Retour

Cover image for bash Cookbook, 2nd Edition QGIS: Becoming a GIS Power User by Alexander Bruy Published by Packt Publishing, 2017
  1. Cover
  2. Table of Contents
  3. QGIS: Becoming a GIS Power User
  4. QGIS: Becoming a GIS Power User
  5. QGIS: Becoming a GIS Power User
  6. Credits
  7. Preface
  8. What you need for this learning path
  9. Who this learning path is for
  10. Reader feedback
  11. Customer support
  12. 1. Module 1
  13. 1. Getting Started with QGIS
  14. Running QGIS for the first time
  15. Introducing the QGIS user interface
  16. Finding help and reporting issues
  17. Summary
  18. 2. Viewing Spatial Data
  19. Dealing with coordinate reference systems
  20. Loading raster files
  21. Loading data from databases
  22. Loading data from OGC web services
  23. Styling raster layers
  24. Styling vector layers
  25. Loading background maps
  26. Dealing with project files
  27. Summary
  28. 3. Data Creation and Editing
  29. Working with feature selection tools
  30. Editing vector geometries
  31. Using measuring tools
  32. Editing attributes
  33. Reprojecting and converting vector and raster data
  34. Joining tabular data
  35. Using temporary scratch layers
  36. Checking for topological errors and fixing them
  37. Adding data to spatial databases
  38. Summary
  39. 4. Spatial Analysis
  40. Combining raster and vector data
  41. Vector and raster analysis with Processing
  42. Leveraging the power of spatial databases
  43. Summary
  44. 5. Creating Great Maps
  45. Labeling
  46. Designing print maps
  47. Presenting your maps online
  48. Summary
  49. 6. Extending QGIS with Python
  50. Getting to know the Python Console
  51. Creating custom geoprocessing scripts using Python
  52. Developing your first plugin
  53. Summary
  54. 2. Module 2
  55. 1. Exploring Places – from Concept to Interface
  56. Acquiring data for geospatial applications
  57. Visualizing GIS data
  58. The basemap
  59. Summary
  60. 2. Identifying the Best Places
  61. Raster analysis
  62. Publishing the results as a web application
  63. Summary
  64. 3. Discovering Physical Relationships
  65. Spatial join for a performant operational layer interaction
  66. The CartoDB platform
  67. Leaflet and an external API: CartoDB SQL
  68. Summary
  69. 4. Finding the Best Way to Get There
  70. OpenStreetMap data for topology
  71. Database importing and topological relationships
  72. Creating the travel time isochron polygons
  73. Generating the shortest paths for all students
  74. Web applications – creating safe corridors
  75. Summary
  76. 5. Demonstrating Change
  77. TopoJSON
  78. The D3 data visualization library
  79. Summary
  80. 6. Estimating Unknown Values
  81. Interpolated model values
  82. A dynamic web application – OpenLayers AJAX with Python and SpatiaLite
  83. Summary
  84. 7. Mapping for Enterprises and Communities
  85. The cartographic rendering of geospatial data – MBTiles and UTFGrid
  86. Interacting with Mapbox services
  87. Putting it all together
  88. Going further – local MBTiles hosting with TileStream
  89. Summary
  90. 3. Module 3
  91. 1. Data Input and Output
  92. Finding geospatial data on your computer
  93. Describing data sources
  94. Importing data from text files
  95. Importing KML/KMZ files
  96. Importing DXF/DWG files
  97. Opening a NetCDF file
  98. Saving a vector layer
  99. Saving a raster layer
  100. Reprojecting a layer
  101. Batch format conversion
  102. Batch reprojection
  103. Loading vector layers into SpatiaLite
  104. Loading vector layers into PostGIS
  105. 2. Data Management
  106. Joining layer data
  107. Cleaning up the attribute table
  108. Configuring relations
  109. Joining tables in databases
  110. Creating views in SpatiaLite
  111. Creating views in PostGIS
  112. Creating spatial indexes
  113. Georeferencing rasters
  114. Georeferencing vector layers
  115. Creating raster overviews (pyramids)
  116. Building virtual rasters (catalogs)
  117. 3. Common Data Preprocessing Steps
  118. Converting points to lines to polygons and back – QGIS
  119. Converting points to lines to polygons and back – SpatiaLite
  120. Converting points to lines to polygons and back – PostGIS
  121. Cropping rasters
  122. Clipping vectors
  123. Extracting vectors
  124. Converting rasters to vectors
  125. Converting vectors to rasters
  126. Building DateTime strings
  127. Geotagging photos
  128. 4. Data Exploration
  129. Listing unique values in a column
  130. Exploring numeric value distribution in a column
  131. Exploring spatiotemporal vector data using Time Manager
  132. Creating animations using Time Manager
  133. Designing time-dependent styles
  134. Loading BaseMaps with the QuickMapServices plugin
  135. Loading BaseMaps with the OpenLayers plugin
  136. Viewing geotagged photos
  137. 5. Classic Vector Analysis
  138. Selecting optimum sites
  139. Dasymetric mapping
  140. Calculating regional statistics
  141. Estimating density heatmaps
  142. Estimating values based on samples
  143. 6. Network Analysis
  144. Creating a simple routing network
  145. Calculating the shortest paths using the Road graph plugin
  146. Routing with one-way streets in the Road graph plugin
  147. Calculating the shortest paths with the QGIS network analysis library
  148. Routing point sequences
  149. Automating multiple route computation using batch processing
  150. Matching points to the nearest line
  151. Creating a routing network for pgRouting
  152. Visualizing the pgRouting results in QGIS
  153. Using the pgRoutingLayer plugin for convenience
  154. Getting network data from the OSM
  155. 7. Raster Analysis I
  156. Using the raster calculator
  157. Preparing elevation data
  158. Calculating a slope
  159. Calculating a hillshade layer
  160. Analyzing hydrology
  161. Calculating a topographic index
  162. Automating analysis tasks using the graphical modeler
  163. 8. Raster Analysis II
  164. Calculating NDVI
  165. Handling null values
  166. Setting extents with masks
  167. Sampling a raster layer
  168. Visualizing multispectral layers
  169. Modifying and reclassifying values in raster layers
  170. Performing supervised classification of raster layers
  171. 9. QGIS and the Web
  172. Using web services
  173. Using WFS and WFS-T
  174. Searching CSW
  175. Using WMS and WMS Tiles
  176. Using WCS
  177. Using GDAL
  178. Serving web maps with the QGIS server
  179. Scale-dependent rendering
  180. Hooking up web clients
  181. Managing GeoServer from QGIS
  182. 10. Cartography Tips
  183. Using Rule Based Rendering
  184. Handling transparencies
  185. Understanding the feature and layer blending modes
  186. Saving and loading styles
  187. Configuring data-defined labels
  188. Creating custom SVG graphics
  189. Making pretty graticules in any projection
  190. Making useful graticules in printed maps
  191. Creating a map series using Atlas
  192. 11. Extending QGIS
  193. Defining custom projections
  194. Working near the dateline
  195. Working offline
  196. Using the QspatiaLite plugin
  197. Adding plugins with Python dependencies
  198. Using the Python console
  199. Writing Processing algorithms
  200. Writing QGIS plugins
  201. Using external tools
  202. 12. Up and Coming
  203. Preparing LiDAR data
  204. Opening File Geodatabases with the OpenFileGDB driver
  205. Using Geopackages
  206. The PostGIS Topology Editor plugin
  207. The Topology Checker plugin
  208. GRASS Topology tools
  209. Hunting for bugs
  210. Reporting bugs
  211. Bibliography
  212. Index

Getting to know the Python Console

The most direct way to interact with the QGIS API (short for Application Programming Interface) is through the Python Console, which can be opened by going to Plugins | Python Console. As you can see in the following screenshot, the Python Console is displayed within a new panel below the map:

Getting to know the Python Console

Our access point for interaction with the application, project, and data is the iface object. To get a list of all the functions available for iface, type help(iface). Alternatively, this information is available online in the API documentation at http://qgis.org/api/classQgisInterface.html.

Loading and exploring datasets

One of the first things we will want to do is to load some data. For example, to load a vector layer, we use the addVectorLayer() function of iface:

v_layer = iface.addVectorLayer('C:/Users/anita/Documents/Geodata/qgis_sample_data/shapefiles/airports.shp','airports','ogr')

When we execute this command, airports.shp will be loaded using the ogr driver and added to the map under the layer name of airports. Additionally, this function returns the created layer object. Using this layer object—which we stored in v_layer—we can access vector layer functions, such as name(), which returns the layer name and is displayed in the Layers list:

v_layer.name()

This is the output:

u'airports'

Note

The u in front of the airports layer name shows that the name is returned as a Unicode string.

Of course, the next logical step is to look at the layer's features. The number of features can be accessed using featureCount():

v_layer.featureCount()

Here is the output:

76L

This shows us that the airport layer contains 76 features. The L in the end shows that it's a numerical value of the long type. In our next step, we will access these features. This is possible using the getFeatures() function, which returns a QgsFeatureIterator object. With a simple for loop, we can then print the attributes() of all features in our layer:

my_features = v_layer.getFeatures()
for feature in my_features:
    print feature.attributes()

This is the output:

[1, u'US00157', 78.0, u'Airport/Airfield', u'PA', u'NOATAK' ...
[2, u'US00229', 264.0, u'Airport/Airfield', u'PA', u'AMBLER'...
[3, u'US00186', 585.0, u'Airport/Airfield', u'PABT', u'BETTL...
...

Tip

When using the preceding code snippet, it is worth noting that the Python syntax requires proper indentation. This means that, for example, the content of the for loop has to be indented, as shown in the preceding code. If Python encounters such errors, it will raise an Indentation Error.

You might have noticed that attributes() shows us the attribute values, but we don't know the field names yet. To get the field names, we use this code:

for field in v_layer.fields():
    print field.name()

The output is as follows:

ID
fk_region
ELEV
NAME
USE

Once we know the field names, we can access specific feature attributes, for example, NAME:

for feature in v_layer.getFeatures():
    print feature.attribute('NAME')

This is the output:

NOATAK
AMBLER
BETTLES
...

A quick solution to, for example, sum up the elevation values is as follows:

sum([feature.attribute('ELEV') for feature in v_layer.getFeatures()])

Here is the output:

22758.0

Note

In the previous example, we took advantage of the fact that Python allows us to create a list by writing a for loop inside square brackets. This is called list comprehension, and you can read more about it at https://docs.python.org/2/tutorial/datastructures.html#list-comprehensions.

Loading raster data is very similar to loading vector data and is done using addRasterLayer():

r_layer = iface.addRasterLayer('C:/Users/anita/Documents/Geodata/qgis_sample_data/raster/SR_50M_alaska_nad.tif','hillshade')
r_layer.name()

The following is the output:

u'hillshade'

To get the raster layer's size in pixels we can use the width() and height() functions, like this:

r_layer.width(), r_layer.height()

Here is the output:

(1754, 1394)

If we want to know more about the raster values, we use the layer's data provider object, which provides access to the raster band statistics. It's worth noting that we have to use bandStatistics(1) instead of bandStatistics(0) to access the statistics of a single-band raster, such as our hillshade layer (for example, for the maximum value):

r_layer.dataProvider().bandStatistics(1).maximumValue

The output is as follows:

251.0

Other values that can be accessed like this are minimumValue, range, stdDev, and sum. For a full list, use this line:

help(r_layer.dataProvider().bandStatistics(1))

Styling layers

Since we now know how to load data, we can continue to style the layers. The simplest option is to load a premade style (a .qml file):

v_layer.loadNamedStyle('C:/temp/planes.qml')
v_layer.triggerRepaint()

Make sure that you call triggerRepaint() to ensure that the map is redrawn to reflect your changes.

Note

You can create planes.qml by saving the airport style you created in Chapter 2, Viewing Spatial Data (by going to Layer Properties | Style | Save Style | QGIS Layer Style File), or use any other style you like.

Of course, we can also create a style in code. Let's take a look at a basic single symbol renderer. We create a simple symbol with one layer, for example, a yellow diamond:

from PyQt4.QtGui import QColor
symbol = QgsMarkerSymbolV2()
symbol.symbolLayer(0).setName('diamond')
symbol.symbolLayer(0).setSize(10)
symbol.symbolLayer(0).setColor(QColor('#ffff00'))
v_layer.rendererV2().setSymbol(symbol)
v_layer.triggerRepaint()

A much more advanced approach is to create a rule-based renderer. We discussed the basics of rule-based renderers in Chapter 5, Creating Great Maps. The following example creates two rules: one for civil-use airports and one for all other airports. Due to the length of this script, I recommend that you use the Python Console editor, which you can open by clicking on the Show editor button, as shown in the following screenshot:

Styling layers

Each rule in this example has a name, a filter expression, and a symbol color. Note how the rules are appended to the renderer's root rule:

from PyQt4.QtGui import QColor
rules = [['Civil','USE LIKE \'%Civil%\'','green'], ['Other','USE NOT LIKE \'%Civil%\'','red']]
symbol = QgsSymbolV2.defaultSymbol(v_layer.geometryType())
renderer = QgsRuleBasedRendererV2(symbol)
root_rule = renderer.rootRule()
for label, expression, color_name in rules:
    rule = root_rule.children()[0].clone()
    rule.setLabel(label)
    rule.setFilterExpression(expression)
    rule.symbol().setColor(QColor(color_name))
    root_rule.appendChild(rule)
root_rule.removeChildAt(0)
v_layer.setRendererV2(renderer)
v_layer.triggerRepaint()

To run the script, click on the Run script button at the bottom of the editor toolbar.

Tip

If you are interested in reading more about styling vector layers, I recommend Joshua Arnott's post at http://snorf.net/blog/2014/03/04/symbology-of-vector-layers-in-qgis-python-plugins/.

Filtering data

To filter vector layer features programmatically, we can specify a subset string. This is the same as defining a Feature subset query in in the Layer Properties | General section. For example, we can choose to display airports only if their names start with an A:

v_layer.setSubsetString("NAME LIKE 'A%'")

To remove the filter, just set an empty subset string:

v_layer.setSubsetString("")

Creating a memory layer

A great way to create a temporary vector layer is by using so-called memory layers. Memory layers are a good option for temporary analysis output or visualizations. They are the scripting equivalent of temporary scratch layers, which we used in Chapter 3, Data Creation and Editing. Like temporary scratch layers, memory layers exist within a QGIS session and are destroyed when QGIS is closed. In the following example, we create a memory layer and add a polygon feature to it.

Basically, a memory layer is a QgsVectorLayer like any other. However, the provider (the third parameter) is not 'ogr' as in the previous example of loading a file, but 'memory'. Instead of a file path, the first parameter is a definition string that specifies the geometry type, the CRS, and the attribute table fields (in this case, one integer field called MYNUM and one string field called MYTXT):

mem_layer = QgsVectorLayer("Polygon?crs=epsg:4326&field=MYNUM:integer&field=MYTXT:string", "temp_layer", "memory")
if not mem_layer.isValid():
    raise Exception("Failed to create memory layer")

Once we have created the QgsVectorLayer object, we can start adding features to its data provider:

mem_layer_provider = mem_layer.dataProvider()
my_polygon = QgsFeature()
my_polygon.setGeometry(QgsGeometry.fromRect(QgsRectangle(16,48,17,49)))
my_polygon.setAttributes([10,"hello world"])
mem_layer_provider.addFeatures([my_polygon])
QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

Note

Note how we first create a blank QgsFeature, to which we then add geometry and attributes using setGeometry() and setAttributes(), respectively. When we add the layer to QgsMapLayerRegistry, the layer is rendered on the map.

Exporting map images

The simplest option for saving the current map is by using the scripting equivalent of Save as Image (under Project). This will export the current map to an image file in the same resolution as the map area in the QGIS application window:

iface.mapCanvas().saveAsImage('C:/temp/simple_export.png')

If we want more control over the size and resolution of the exported image, we need a few more lines of code. The following example shows how we can create our own QgsMapRendererCustomPainterJob object and configure to our own liking using custom QgsMapSettings for size (width and height), resolution (dpi), map extent, and map layers:

from PyQt4.QtGui import QImage, QPainter
from PyQt4.QtCore import QSize
# configure the output image
width = 800
height = 600
dpi = 92
img = QImage(QSize(width, height), QImage.Format_RGB32)
img.setDotsPerMeterX(dpi / 25.4 * 1000)
img.setDotsPerMeterY(dpi / 25.4 * 1000)
# get the map layers and extent
layers = [ layer.id() for layer in iface.legendInterface().layers() ]
extent = iface.mapCanvas().extent()
# configure map settings for export
mapSettings = QgsMapSettings()
mapSettings.setMapUnits(0)
mapSettings.setExtent(extent)
mapSettings.setOutputDpi(dpi)
mapSettings.setOutputSize(QSize(width, height))
mapSettings.setLayers(layers)
mapSettings.setFlags(QgsMapSettings.Antialiasing | QgsMapSettings.UseAdvancedEffects | QgsMapSettings.ForceVectorOutput | QgsMapSettings.DrawLabeling)
# configure and run painter
p = QPainter()
p.begin(img)
mapRenderer = QgsMapRendererCustomPainterJob(mapSettings, p)
mapRenderer.start()
mapRenderer.waitForFinished()
p.end()
# save the result
img.save("C:/temp/custom_export.png","png")