It is no secret that one of the most common raster formats, TIFF, is able to store georeferencing information within the file, turning it into a GeoTIFF. Everything would be fine if it weren’t for one nuance: the vast majority of graphic editors suffer from “star fever”, believing that unrecognised tags have no place in the file and happily removing this valuable information when saving the file.
This feature of graphic editors is often used to “reset” georeferencing (for example, if the image is incorrectly georeferenced and needs to be georeferenced from scratch). A natural question arises: what to do if you need to edit the image in a graphics editor but do not want to lose the georeferencing? The answer is not original, it is necessary to save it in an external file (the so-called world file).
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 Germán’s post “Layer list widget for PyQGIS applications” on the GeoTux blog.
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.
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:
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.
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.
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.
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"] =Falsefp = 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 plotx = arange(0.0, 3.0, 0.01)
y1 = sin(2* pi * x)
y2 = cos(2* pi * x)
# the main axes is subplot(111) by defaultplot(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()
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:
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"] =Falsefp = 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.