QGIS Planet

Getting multipolygon vertexes using PostGIS

EN | PT

Today I needed to create a view in PostGIS that returned the vertexes of a multi-polygon layer. Besides, I needed that they were numerically ordered starting in 1, and with the respective XY coordinates.

Screenshot from 2015-11-05 23:58:19

It seemed to be a trivial task – all I would need was to use the ST_DumpPoints() function to get all vertexes – if it wasn’t for the fact that PostGIS polygons have a duplicate vertex (the last vertex must be equal to the first one) that I have no interess in showing.

After some try-and-fail, I came up with the following query:

CREATE OR REPLACE VIEW public.my_polygons_vertexes AS
WITH t AS -- Transfor polygons in sets of points
    (SELECT id_polygon,
            st_dumppoints(geom) AS dump
     FROM public.my_polygons),
f AS -- Get the geometry and the indexes from the sets of points 
    (SELECT t.id_polygon,
           (t.dump).path[1] AS part,
           (t.dump).path[3] AS vertex,
           (t.dump).geom AS geom
     FROM t)
-- Get all points filtering the last point for each geometry part
SELECT row_number() OVER () AS gid, -- Creating a unique id
       f.id_polygon,
       f.part,
       f.vertex,
       ST_X(f.geom) as x, -- Get point's X coordinate
       ST_Y(f.geom) as y, -- Get point's Y coordinate
       f.geom::geometry('POINT',4326) as geom -- make sure of the resulting geometry type
FROM f 
WHERE (f.id_polygon, f.part, f.vertex) NOT IN
      (SELECT f.id_polygon,
              f.part,
              max(f.vertex) AS max
       FROM f
       GROUP BY f.id_polygon,
                f.part);

The interesting part occurs in the WHERE clause, basically, from the list of all vertexes, only the ones not included in the list of vertexes with the maximum index by polygon part are showed, that is, the last vertex of each polygon part.

Here’s the result:

Screenshot from 2015-11-05 23:58:40

The advantage of this approach (using PostGIS) instead of using “Polygons to Lines” and “Lines to points” processing tools is that we just need to change the polygons layer, and save it, to see our vertexes get updated automatically. It’s because of this kind of stuff that I love PostGIS.

Learn More

Creating a Tissot’s Indicatrix in QGIS

The task of projecting, or unfolding the spherical Earth onto a flat map, is an age old problem in cartography. Projection almost always introduces distortion, most projections cannot preserve angles, areas and distances at the same time, they may be conformal (angle-preserving), equal-area (area-preserving) or equidistant (distance preserving) but not all at once. The only exception is a Globe, which preserves angles, areas and distances perfectly. Thus a projection is a compromise.

The choice of projection depends on a map’s use, scale and audience. Conformal projections, for example, are preferred for nautical charts or small scale maps because they locally preserve angles necessary for navigation and survey drawings. Equal-area projections are best suited for maps of broad continental region as they preserve the relative sizes of countries, seas and oceans and allow comparison between regions. Finally, there are hybrid projections that minimise the distortion by merging conformal and equal-area projections, these can be used to create visually pleasing maps of the entire Earth (for a guide to selecting a map projection see Fig. 9 in Jenny (2011), link below).

But how does one measure the degree and type of distortion in a map projection?

One elegant method was developed in the 1880s century by the French cartographer Monsieur Nicolas Auguste Tissot, the Tissot’s Indicatrix (or Tissot’s Ellipse). This mathematical contrivance consists of a grid of infinitely small circles that measures the degree and type of distortion caused by projection. While Monsieur Tissot’s approach is mathematical, involving infinity small circles, his technique can be approximated overlaying a regular grid of large circles and crosses to a map.

The Indicatrix Mapper plugin for QGIS by Ervin Wirth and Peter Kun creates a Tissot’s Indicatrix by adding a vector layer of circles and crosses in a gridded pattern on a map. The degree and type of distortion of the Tissot’s Indicatrix reveals the class of map projection as follows: –

  • If a projection is conformal, the area of circles and sizes of the crosses will change while the shapes of circles remain the same and intersection angle of the crosses will always meet at 90 degrees e.g. Mercator projection
  • If a projection is equal-area, the area of the circles will remain the same while the shape of the circles change and intersection angle of the crosses will not always meet at 90 degrees e.g. Mollweide projection or Hammer projection
  • If a projection preserves neither property, the area of the circles and their shape will change, and the intersection angle of the crosses will not always meet at 90 degrees e.g. Robinson

After adding the Indicatrix Mapper plugin to QGIS (menu Plugins – Manage and Install Plugins) first add a basemap using the OpenLayers plugin e.g. Bing Aerial layer, then click the Indicatrix Mapper icon and run the plugin using default settings. You can then select different projections (lower right in world icon QGIS) to see the effects of various protections on the Tissot Indicatrix. If the circles appear as squares after selecting a different projection, right click the Circles layer in the layers panel, then select the Rendering tab and deselect the Simplify geometry check box. Also, turn off the basemap layer when using different projections, unfortunately the OpenLayers plugin only supports Google Mercator projection (EPSG 3857). To create the basemap below, that can be displayed using different projections, I styled vector data downloaded from Natural Earth and OpenStreetMap.

Mercator

Mercator Projection – the area of the circles and size of the crosses increase towards the poles but their shape remains the same.

Mollweide

Mollweide Projection – the area of the circles remain the same but their shapes are distorted, the crosses do not always intersect at 90 degrees.

Robinson

Robinson Projection – both the area of the circles and intersection angle of the crosses circles vary.

It is important to note that a Tissot’s Indicatrix generated in QGIS is an approximation of mathematical ideal, we are not no longer dealing with infinity small circles. As a result, here will be some minor distortion visible towards the edge of a map independent of the projection used; notice that the circles in the Mercator projection nearest the poles are not quite symmetrical or the circles at the edge of the Mollweide projection do not appear to preserve area as they should. This anomalous distortion can be minimised by reducing the size and spacing of the circles and crosses created by the Indicatrix Mapper plugin. However, despite these limitations a Tissot’s Indicatrix elegantly reveals the distortion present. This is something to important to understand when when choosing a map projection.

References:

Jenny, B., 2012. Adaptive composite map projections [PDF]. Visualization and Computer Graphics, IEEE Transactions on, 18, 2575–2582.
Learn More

Rule-based labelling in QGIS 2.12

The new QGIS 2.12 (Lyon) will be out soon!

In this release, we have revamped the labelling engine and made it more flexible in-line with the rest of vector styling.

In this release, we have:

  • Revamped the labelling engine
  • Added support for mutually exclusive layer tree groups
  • Developed raster alignment tool

For users

In previous versions of QGIS, users can select a field value or use an expression as labels for a vector layer. In QGIS 2.12, users can additionally define rules to label vector layers. Rule-based labelling works in the similar way as “Style”. A list of rules will be defined by users and they will be applied from top-to-bottom.

Field based labelling (Click to enlarge)

Rule based labelling (Click to enlarge)

To achieve the same effect in the earlier versions of QGIS, users should add a complex expression.The example below shows the expression used in earlier versions of QGIS:

CASE WHEN  length( "htmlname" ) > 13 AND strpos("htmlname",' ') > 6  THEN  replace("htmlname",' ','  ') WHEN  length( "htmlname" ) > 20 AND "htmlname"  LIKE '%Golf Course' THEN  regexp_replace("htmlname",'Golf Course',' Golf Course') WHEN  length( "htmlname" ) > 20 AND "htmlname"  LIKE '%Nature Reserve' THEN  regexp_replace("htmlname",'Nature Reserve',' Nature Reserve') WHEN  length( "htmlname" ) > 20 AND "htmlname"  LIKE '%Church Of England%' THEN  regexp_replace("htmlname",'Church Of England',' Church Of England ')  WHEN  length( "htmlname" ) > 13 AND "htmlname"  LIKE '% Of The %' THEN  regexp_replace("htmlname",'Of The','Of The ') WHEN  length( "htmlname" ) > 13 AND "htmlname"  LIKE '% of %' AND  "fontcolour" <> 2 AND  "fontcolour" <>  4 THEN  regexp_replace("htmlname",' of ',' of  ')  WHEN "htmlname" LIKE '%/%' THEN regexp_replace("htmlname",'/','/  ') WHEN  length( "htmlname" ) > 30 THEN  replace("htmlname",' ','  ')  ELSE  "htmlname"  END

As you can see, without rule-based labelling users have to define each case and the labelling text. Additionally, with the cumbersome task of defining each case, users are limited to using specific labelling formats (e.g. font color, font size, visibility range, etc). With the new labelling engine, users can define their labels similarly to styles - in fact they can copy the rules from styling tabs and use it within the new labelling section!

Examples of label rules (Click to enlarge)

For developers

In order to facilitate addition of rule-based labelling, some internal changes were made to the QGIS labelling engine interface. The labelling is now driven by the new class QgsLabelingEngineV2 which may have several label providers associated with it. The label providers are objects derived from QgsAbstractLabelProvider and they are responsible for:

  1. providing “label features” that define properties of each label and geometry of the feature they represent; and

  2. drawing of labels at the positions determined by the engine.

Currently there are label provider implementations for diagrams, simple labelling and rule-based labelling.

The existing labelling engine (QgsPalLabeling class) is now built on top of the new labelling engine and works as a wrapper for it so that existing code that uses QgsPalLabeling still works.

As of now, the API for the new labelling engine is not considered as complete and therefore not available in Python. The idea is to make it easier to use, more polished and better prepared for the future use cases. It will be likely finished during the 2.14 release cycle and some changes may need to be postponed to QGIS 3.0 where backwards incompatible API changes will be allowed (as of now, QGIS 3.0 is being intensively discussed on QGIS developer mailing list).

Funding

This rule-based labeling has been funded by Tuscany Region (Italy). Special thanks also to the QGIS developers for their help with bug fixing.

You may also like...

Mergin Maps, a field data collection app based on QGIS. Mergin Maps makes field work easy with its simple interface and cloud-based sync. Available on Android, iOS and Windows. Screenshots of the Mergin Maps mobile app for Field Data Collection
Get it on Google Play Get it on Apple store
Learn More

Presentations at FOSS4G 2015 in Seoul

Slides from our presentations at FOSS4G 2015 in Seoul: Keynote: The QGIS project and its evolution from a desktop GIS to a GIS platform - Slides New QGIS functions for power users - Slides QGIS Plugins - From Must-Haves to insider tips - Slides Building an OpenLayers 3 map viewer with React - Slides Thanks to the organizers of this great conference! It was a pleasure to get in contact with so many users from around the world.
Learn More

Final GSoC report

What do I have completed this week?

  • Added a dedicated threadPool for the Processing Toolbox to run the algorithms.

  • Bugfix on the cancel option. Avoid connecting  the signals every time the algorithm runs. This was making the slots to be called multiple times when running and cancelling the algorithm several times.

Is there any blocking issue?

  • I wasn’t able to fix the issue with the cancel option when you have more than one algorithm running. The issue turned out to be more complicated than I first expected.

Learn More

Parallel execution of QGIS algorithms

This is just an example of the advantage provided by the multithreading on the Toolbox. As each algorithm is running in a separate thread, this allows us to run more than one algorithm at the same time.

 

I’m polishing the implementation and trying to find some major problems in the current implementation.  It has some minor bugs that can be easily solved by adding a parameter to the  cancel signal (e.g. when you have two algorithms running and try to cancel one of them, the signal emitted will cancel two algorithms instead of just one)

Learn More

Twelfth GSoC report – Multithreading on Processing

What do I have completed this week?

  • I have finished the cancel option in order to stop the execution of non-QGIS algorithms. It’s important to remember that as well as the QGIS algorithms the cancel option do not cover all the non-QGIS algorithms and this implementation only covers GDAL and SAGA algorithms. However this option can be easily replicated to the othersthird party algorithms.

  • I have changed the python console to use the QThreadPool.

  • I have been testing and debugging the multithreading implementation looking for unexpected behaviours or incorrect results.

  • Core refactoring.

What am I going to achieve for the next week?

  • Write documentation for the multithreading implementation.

  • Code cleaning and refactoring.

  • Search for problems in the multithreading implementation and correct them.

Is there any blocking issue?

  • No.

Learn More

Eleventh GSoC report – Multithreading on Processing

What do I have completed this week?

  • During this week I have changed the multithreading implementation to use the QThreadPool to allow thread recycling and avoid thread creation costs every time we want to run an algorithm. This new implementation adds some flexibility to the multithread support and avoid unexpected behaviours in different machines. The number of threads in the QThreadPool is initialised accordingly to the number of cores in the PC and may happen that the QT is not able to detect the number of cores. In this case we will only have one thread in the thread pool, which can be occupied with something else. I’m trying to detect when there is no available threads on the ThreadPool and increase the number of threads if necessary.

What am I going to achieve for the next week?

  • Start working on a test suite for the multithreading to ensure the correct behaviour of the multithreading implementation. The test suite must check if the worker thread is starting, check if signals are triggering the functions when it is supposed to and verify if the cancel option is working properly.

  • Figure out how to stop the execution of a non-QGIS algorithm.

Is there any blocking issue?

  • I have spotted some issues when testing the code in a different machine (with Linux Mint 17) that may have to do with thread unavailability.  Therefore, I decided to postpone the current tasks and focus on changing the multithreading implementation in order to correct this kind of unexpected behaviours.

Learn More

Tenth GSoC report – Multithreading on Processing

What do I have completed this week?

  • During this week I worked on a mechanism to stop the algorithm execution on QGIS algorithms. The approach used to cancel the execution explores python exceptions and QT signals to stop the algorithm from the main thread without having to wait in order to stop the execution and terminate the worker thread. I have used signals to trigger an exception inside the algorithm and stop the execution from the main thread. Since there is a significant amount of algorithms, the mechanism had to be implemented in just two of the QGIS algorithms. This implementation can be further replicated  towards the remaining  QGIS algorithms.

What am I going to achieve for the next week?

  • Start working on a test suite for the multithreading to ensure the correct behaviour of the multithreading implementation. The test suite must check if the worker thread is starting, check if signals are triggering the functions when it is supposed to and verify if the cancel option is working properly.

  • Figure out how to stop the execution of a non-QGIS algorithm.

Is there any blocking issue?

  • No blocking issues
Learn More

Ninth GSoC Report – Multithreading on Processing

What do I have completed this week?

  • Solved an issue that makes QGIS crash when the main thread doesn’t wait for the worker thread to start.
  • I spent most of the time looking up some tutorials about unit testing on QGIS and trying to configure python paths to run the processing tests.

What am I going to achieve for the next week?

  • This week I’ll change back to the multithreading implementation in order to add the cancel option to stop the algorithm execution. I will be working on a mechanism to stop the execution on QGIS algorithm.

Is there any blocking issue?

  • This week my mentor proposed to start working on a suite test and I changed the plan to start working on a unit test to ensure the correct behaviour of the multithreading implementation. At the end of the week we discussed about asking the community to test the new implementation and report the issues if there is any.
Learn More