Category: GIS

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 their coordinates. 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 polygon 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.

QGIS 1.5 User Guide translation

25.08.2010 16:37 ·  GIS  ·  qgis

A collaborative project to translate the QGIS 1.5 User Guide into Russian is announced. I will lead and coordinate the project. Please share information and get involved, and don’t hesitate to contact me if you have any questions.

QGIS mapserver

20.08.2010 15:17 ·  GIS  ·  qgis

Marco Hugentobler wants to merge his project, QGIS mapserver, into the QGIS source code.

QGIS mapserver is an open source WMS server that runs on Linux, Windows and Mac OS X (other Unix-like systems should also be supported, but have not been tested). The server is a FastCGI application written in C++. It uses QGIS libraries for rendering and GIS operations.

Among the main features:

QGIS 1.5 "Tethys"

06.08.2010 10:30 ·  Notes, GIS  ·  qgis, release

I returned to Zaporizhzhia yesterday, sunburned and a little refreshed. The heat in Crimea is terrible, and the proximity of the sea does not help. In Choban-Kale, at 11 o’clock in the afternoon, we had 53°C in the direct sun (they say it was almost 60°C in Simeiz), and the water temperature was ~28–29°C. However, it is no better in Zaporizhzhia.

While I was on holiday, the next version of QGIS — 1.5.0 “Tethys” — was released. The official announcement is on the QGIS blog. There are many changes: a lot of bugs have been fixed, new tools have been added, and old tools have been improved (fTools, GdalTools, annotations, new georeferencing module). The documentation is actively updated.

Table joins in QGIS

22.07.2010 14:44 ·  GIS  ·  qgis

Marco Hugentobler has implemented initial support for table joins in QGIS. If you want to try out the new functionality, you can get the source code from the new branch of the repository

svn co https://svn.osgeo.org/qgis/branches/table_join_branch table_join

and build it.

To join a table to a layer’s attribute table, you need:

As this feature is still in the early stages of development, there are some issues:

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.

GdalTools — QGIS plugin for raster data processing

06.06.2010 13:07 ·  GIS  ·  qgis, plugins, gdaltools

GdalTools (or Raster Tools) is a plugin for the open source GIS QGIS. The main purpose of the plugin is to simplify the use of the GDAL command line utilities by providing the user with a graphical interface for the most common operations.

Initially, the plugin was developed by Faunalia, later other developers, including myself, joined them. Robert Szczepanek created nice icons for most of the tools, and the development of some features was sponsored by Silvio Grosso.

At the time of the first announcement (September 2009), the plugin provided only two tools: gdalbuildvrt and gdal_contour. Now the number of integrated tools is approaching 20.

Read more ››

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 ››

GdalTools update

03.02.2010 15:36 ·  GIS  ·  qgis, plugins, gdaltools

Today we merged the experimental branch of the GdalTools plugin into the main development branch.

GdalTools (aka Raster Tools) provides users with a simple graphical interface to perform the most common raster processing tasks. Originally the plugin was created by Faunalia (Paolo Cavallini, Giuseppe Sucameli and Lorenzo Masini), the icons for the extension were created by Robert Szczepanek. About a month ago I also joined the work (and this experimental branch is my work).

This is what we ended up with: