Tag: tips

Correct scale for tile maps in QGIS

26.12.2014 19:01 ·  GIS  ·  qgis, tips

When viewing TMS layers in QGIS, as well as in any other GIS, the background map may be blurred if the current map scale does not match the scale of the tiles. The reason for this is a feature of TMS, namely the use of a fixed set of so-called “zoom levels”: tiles are generated only for certain scales defined by the data provider.

So to get a sharp image, you should only use the scales that correspond to the zoom levels of the selected TMS service. For OpenStreetMap, the formula for calculating the scale for each zoom level can be found on this page, and a similar approach can be used for other services.

Some time ago I added a special widget to QGIS to select a map scale from a given set. Later, it became possible to edit this list and define a set of scales on a global level as well as on a project level. So if your project uses TMS layers, you can create your own list of scales and switch between them. At the same time, you still have the option of using any intermediate scale values.

Another option is to install the Tile Map Scale plugin. This plugin allows you to easily connect popular TMS layers (remember about ToS!). It also monitors map scale changes and automatically sets the nearest correct TMS scale. Note that you will not be able to use intermediate scale values in this case.

Finally, QGIS has a built-in widget for changing the scale according to the layer’s scale list. It can be found in “View → Panels → Tile scale” or from the toolbar context menu.

The choice of method is up to you.

Automatic stretching of the raster histogram in QGIS

01.03.2011 15:01 ·  GIS  ·  qgis, tips

QGIS allows you to save a default raster band combination, standard deviation and contrast enhancement algorithm for rasters.

But it seems that not many people use this feature. At least no one noticed that although the band combination and standard deviation can be saved, they are of no use because the saved values are not applied when the raster is loaded. There were also other problems (for example, the default band combination is 5-4-3, but the raster properties show 1-2-3) that no one reported.

All this has been around for a long time, somewhere since r13582 or even earlier. And only in r15256 were all these problems finally solved.

I’ll try to explain why you might need it and how to use it.

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.

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.