Tag: python

Adding a layer list widget to a PyQGIS application

14.02.2011 09:31 ·  GIS  ·  qgis, python, howto

QGIS is not only a ready-to-use desktop GIS, but also a set of libraries that can be used to create custom GIS applications. Unfortunately, when creating such applications, developers have to implement some GUI elements from scratch. An example of such a widget is the list of loaded layers (sometimes incorrectly called a legend).

In this post I will show the process of embedding a layer list widget into a PyQGIS application. The widget was developed by Germán Carrillo and you can get the code and read how to use it in German’s post “Layer list widget for PyQGIS applications” on the GeoTux blog.

Read more ››

Extracting vertices from polygon geometry

12.10.2010 18:01 ·  GIS  ·  gdal, python, tips

Imagine you have a vector layer containing polygons, both single- and multipart. Your task is to get vertices of each geometry, including all rings and parts. Here is how it can be done using GDAL and Python.

Let’s start with simplest case — singlepart polygons. In that case everything is straightforward — we need to find out the number of vertices, iterate over them and get coordinated. For demostration purposes the code will look like this:

def dump_geometry(geom):
    print "Exterior ring"
    ring = geom.GetGeometryRef(0)
    for i in range(ring.GetPointCount()):
        print i, ring.GetX(i), ring.GetY(i)

A polygon with “holes” is considered as a set of “rings”: one outer ring defining the boundary and several inner rings describing the “holes”. The outer ring always comes first. Thus, the algorithm will transform as follows: loop over all rings; if it is the first iteration - this is an outer ring, otherwise - inner ring; extract vertices from the ring and output their coordinates.

def dump_geometry(geom):
    nRings = geom.GetGeometryCount()
    for i in range(nRings):
        ring = geom.GetGeometryRef(i)
        if i == 0:
            print "Exterior ring"
        else:
            print "Interior ring", i

        for j in range(ring.GetPointCount()):
            print j, ring.GetX(j), ring.GetY(j)

Handling multipart polygons is more complicated: they can be either singlepart polygons or contain “holes”, “islands”, or even “holes” and “islands” at the same time. Let’s try to figure out how to deal with them. A polygon consisting of several “islands” can be considered as a set of independent polygons. In turn, each of those polygons can have one or more rings. So what we have now:

  1. if the polygon is a singlepart polygon and contains only one ring, we extract its vertices immediately
  2. if the number of rings in polygon is greater than 1 — we loop through the rings and extract vertices from them
  3. if the polygon is a multipart polygon — loop through the parts and for each for part perform checks starting from the step 1

So, this is a classical recursion and the final function will look like this:

def dump_geometry(geom, file):
    geomType = geom.GetGeometryType() # geometry type
    # singlepart polygon
    if geomType == ogr.wkbPolygon or geomType == ogr.wkbPolygon25D:
        nRings = geom.GetGeometryCount() # number of rings
        print "Rings", nRings
        if nRings == 1: # there are no "holes" in the polygons
            print "Exterior ring"
            ring = geom.GetGeometryRef(0)
            for i in range(ring.GetPointCount()):
                print i, ring.GetX(i), ring.GetY(i)
        else: # polygon with "holes"
            for i in range(nRings):
                ring = geom.GetGeometryRef(i)
                if i == 0:
                    print "Exterior ring"
                else:
                    print "Interior ring", i

                for j in range(ring.GetPointCount()):
                    print j, ring.GetX(j), ring.GetY(j)
    # multipart polygon
    elif geomType == ogr.wkbMultiPolygon or geomType == ogr.wkbMultiPolygon25D:
        for i in range(geom.GetGeometryCount()):
            subGeom = geom.GetGeometryRef(i) # for each polygon part
            dump_geometry(subGeom) # recursively call dump_geometry() function

The code is a bit bloated, but my aim was to demonstrate a possible way of approaching the task, not to come up with an ideal solution.

Geographic coordinates to raster (row, column) coordinates

21.07.2010 15:37 ·  GIS  ·  qgis, python, howto

I have noticed that recently many people have asked how to get image coordinates (row, column) from real world coordinates (latitude/longitude). The following code shows how to do this in the QGIS Python console:

import math # needed for floor() function
iface = qgis.utils.iface
# let's assume that raster is a current layer
layer = iface.mapCanvas().currentLayer()
# geographic coordinates
x = 75.0791666
y = 57.2541666
extent = layer.extent()
width = layer.width()
height = layer.height()
# raster resolution (pixel size)
xres = (extent.xMaximum() - extent.xMinimum()) / width
yres = (extent.yMaximum() - extent.yMinimum()) / height
# calculate raster coordinates (row and column)
col = int(math.floor((x - extent.xMinimum()) / xres))
row = int(math.floor((extent.yMaximum() - y) / yres))

Nothing complicated or fancy. The same code can be used in plugin with minimal changes.

Working with vectors using GDAL and Python

15.04.2010 09:55 ·  GIS  ·  gdal, python, howto

GDAL is a free library for working with raster and vector data, OGR is a part of the GDAL and is used to work with vector data. The command line utilities included in the library are widely used to perform a variety of tasks. Thanks to the developed API, you can work with OGR functions from many programming languages. This article is dedicated to using the OGR API in Python and is based on the GDAL Vector API Tutorial.

Read more ››

Working with rasters using GDAL and Python

01.04.2010 09:10 ·  GIS  ·  gdal, python, howto

GDAL is a free library for working with raster and vector data. The command line utilities included in the library are widely used to perform a variety of tasks. Thanks to the developed API, you can work with GDAL functions from many programming languages. This article is dedicated to using the GDAL API in Python and is based on the GDAL Raster API Tutorial.

Read more ››

About rasterlang and GDAL

06.01.2010 15:15 ·  GIS  ·  qgis, gdal, python

The rasterlang QGIS plugin refused to work, complaining about a self-compiled GDAL with almost everything enabled. Investigation showed that in addition to all GDAL dependencies, the NumPy package must also be installed. Otherwise you some things will not work.

NumPy is needed to build gdal_array: a small but essential module used to represent images as arrays and then do all sorts of naughty things with them.

After installing NumPy, rebuilding and updating the GDAL package, everything worked.

Again about non-ASCII characters in matplotlib

04.08.2009 08:07 ·  Notes  ·  matplotlib, python, tips

A bit more about Cyrillic and other non-ASCII characters in matplotlib. I decided to create a more comprehensive example to show how to output Cyrillic (or any other non-ASCII) characters in different parts of the plot. Actually, there is nothing complex here; you just need to find the time and read the manual, which is quite good, by the way. But as not everyone likes to read manuals… Anyway, here is the code:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from pylab import *
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm

plt.rcParams["text.usetex"] = False
fp = fm.FontProperties(fname="/home/alex/devel/mpl_cyr/CharisSILR.ttf")
plt.text(0.5, 0.5, u"довільний текст", fontproperties=fp)

# create some data to use for the plot
x = arange(0.0, 3.0, 0.01)
y1 = sin(2 * pi * x)
y2 = cos(2 * pi * x)

# the main axes is subplot(111) by default
plot(x, y1, label=u"Синусоїда")
plot(x, y2, label=u"Косинусоїда")
xlabel(u"Підпис вісі X", fontproperties=fp)
ylabel(u"Підпис вісі Y", fontproperties=fp)
title(u"Назва графіка", fontproperties=fp)
legend(prop=fp)
show()

Cyrillic and other non-ASCII characters in matplotlib

23.07.2009 17:49 ·  Notes  ·  matplotlib, python, tips

There is a wonderful Python library - matplotlib. It is a plotting library that supports a wide range of plot types and is designed to emulate MATLAB commands and behaviour. The library is easy to learn; to draw a simple plot, you literally need two commands. I have used this library in my Statist plugin for QGIS and in another GIS project. To make installation of matplotlib more convenient for inexperienced users, I recently packaged it for OSGeo4W.

Sometimes, I needed to display non-ASCII (namely Cyrillic) characters on matplotlib plots. And there was a problem: such text was drawn as empty squares. Reading manuals, googling, and asking on the mailing list led to two solutions that I would like to share with you.

Method 1: the almighty TeX

matplotlib can use LaTeX to display both plain text and mathematical symbols. Moreover, a limited subset of TeX and the corresponding parser, fonts, and renderer are built into the library, so for this subset, you do not even need to have a full TeX installation. Unfortunately, this TeX subset only contains mathematical characters and letters of the Greek alphabet. In all other cases, an external LaTeX installation is required. To use LaTeX for text rendering, we should set the option

text.usetex: True

in the rc-file. This can be done either globally, by editing the rc-file once, or as needed at runtime. Below is an example of run-time initialisation:

# -*- coding: utf-8 -*-
from matplotlib import rc

rc("font", {"family": "serif"})
rc("text", usetex=True)
rc("text.latex", unicode=True)
rc("text.latex", preamble="usepackage[utf8]{inputenc}")
rc("text.latex", preamble="usepackage[russian]{babel}")

Now we can output Cyrillic (or any other non-ASCII characters)

xlabel(u"Вісь Х: довжина, см")

The big disadvantage of this method is that the user needs to have LaTeX installed.

Method 2: unicode + fonts

matplotlib uses its own font rendering engine with full Unicode support. Therefore, we can explicitly specify a font that contains the required character sets and render the text using that font. Here is a small example:

# -*- coding: utf-8 -*-

import matplotlib.pyplot as plt
import matplotlib.font_manager as fm

plt.rcParams["text.usetex"] = False
fp = fm.FontProperties(fname="/home/alex/.fonts/academy.ttf")
plt.text(0.5, 0.5, u"кириличний текст", fontproperties=fp)
plt.show()

There is also a disadvantage - the required font may not be available on the target system or may be in a different directory, so the font has to be stored in the same directory as your program. But in my opinion, this is much better than having LaTeX as a dependency.

Personally, I chose the latter option.