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:
defdump_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.
defdump_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:
if the polygon is a singlepart polygon and contains only one ring, we extract its vertices immediately
if the number of rings in polygon is greater than 1 — we loop through the rings and extract vertices from them
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:
defdump_geometry(geom, file):
geomType = geom.GetGeometryType() # geometry type# singlepart polygonif 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 polygonelif 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.
The collaborative project to translate the QGIS 1.5 User Guide into Russian is nearing completion. I’m busy with the final editing and formatting.
Recently, I had to look into the OGC CSW (Catalogue Service Web) standard. It’s an interesting and useful thing, but most of the public servers I managed to find implement various aspects of it in a strange way :-(. So I have to add a monstrous construction to the code to catch exceptions and give the user a more or less understandable error message.
I am also slowly writing Pascal wrappers for the GDAL/OGR and Proj libraries. I couldn’t find any Pascal bindings for GDAL/OGR, so I’m making my own. I already have a more or less working module for OGR and some examples of its use. A module for OSR is also ready, but has not been tested yet. Next in line is a module for GDAL.
The situation with the Proj wrapper was a bit more interesting: first, I found a Pascal module for Proj 4.4.3 (the current Proj version is 4.7.0), and this wrapper required changes to the source code of the Proj library. I didn’t like this approach and started to make my own module, partly based on the one I found. Later, when updating the FPC compiler, I accidentally noticed that the available packages included modules for Proj 4.6.1. Now I’m not sure what to do: should I continue developing my own wrapper or not?
I plan to push the results of my work upstream, and if that doesn’t work, I’ll just publish them somewhere.
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.
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:
support for GetCapabilities, GetMap, GetStyle, GetFeatureInfo requests and customisable styles based on the Styled Layer Descriptor specification (supported standards include WMS 1.3.0, WMS 1.1.1 and SLD 1.0.0).
SOAP over HTTP POST. Compatibility with ORCHESTRA and SANY Service Oriented Architecture
SLD extensions (diagrams and custom symbols in SVG format)
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.
I am leaving for Crimea tomorrow morning at 6 a.m. and will be back late in the afternoon on the 5th.
I will not have internet access there, but you can write to me :-). I promise to answer all emails, comments, bug reports and feature requests when I’m back.
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:
add a vector layer to a project
add a table to a project (e.g. by going to the “Layer → Add vector layer” menu and selecting a *.csv or *.dbf file)
open the vector layer properties dialog and select a “Join” section. Then click on a “+” button to add a new joinJoining a table to a vector layer
select the fields that will be used to create the join and press “OK” button
new fields from the joined table should appear in the layer’s attribute table, as well as in the “Attributes” tab of the layer properties dialogue and in the “Identify Results” dialogue
As this feature is still in the early stages of development, there are some issues:
joined fields are read-only, their values cannot be changed
poor performance, especially when performing classification or searching on joined fields
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:
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.
