Blog GIS & Territories

We’ll take up the example covered in the article Eliminating overlays and gaps between polygons in a layer (with QGis and Postgis) to see the same type of solution when you don’t have a PostgreSQL/Postgis database. Here, we’ll be using only the possibilities offered by QGis, working on a Geopackage data format.

If you have the opportunity to work with French cadastral data corresponding to localities, you’ll notice that spatial consistency (topology) is not respected:

In fact, the boundary of each locality is the edge of the adjoining road, leaving empty spaces where the roads are. While this may not be a problem in some cases, it becomes one when you add the “Communes” layer (from the same cadastre) to your map:

At this scale it doesn’t look very clean, but if you zoom in it’s even worse.

And if we look at the locality layer itself, we find numerous topological anomalies:

Here’s how to solve the first problem: How to correct the topology of the lieux-dits layer by removing polygon overlaps and empty spaces between them.

For the second problem: How to recreate a common layer consistent with the new lieux-dits layer and superimposed exactly on the lieux-dits boundaries, please refer to the article cited above (Eliminating superimpositions and gaps between polygons in a layer(with QGis and Postgis), as this part simply uses QGis and Geopackage.

The following example is based on data from the official cadastre of the 71 département, downloaded from https://cadastre.data.gouv.fr/data/etalab-cadastre/2023-01-01/shp/departements/71/

How to correct the topology of the lieux-dits layer

The solution proposed here does not require a PostgreSQL/Postgis database. However, given the limitations of SQLite3 and a bug in QGis 3.30’s DB Manager, some additional manipulations will be required. Long processing times are also to be expected, which is why it is not suitable for processing large layers.

1-Create a layer with a buffer around localities

To manage the spaces between polygons, even if certain snapping tools can be used to deal with certain problems, the simplest and most radical solution is to build a buffer around the polygons, so as to transform empty spaces into superimpositions. In the case of localities, empty spaces are of the order of 10m. We build a 6m buffer around the localities, giving an overlay of the order of 2-3m.

Choose a permanent output file in Geopackage format, which will allow you to follow the procedure in the next step directly without modification.

The result of the buffer layer is as follows:

We no longer have to manage spaces between polygons, just overlays.

2- Open QGis DataBase Manager

You will see the DB Manager window

If your geopackage file doesn’t appear in the list, right-click on Geopackage->New connection and load it.

Click on the layer buffered in the previous step, then on the Database menu ->SQL window.

3- Execute the SQL script lines

All the lines in the following script correspond to the following tasks

1. Create a view with two collections of geometries. Initially, these are the polygons of the localities in duplicate.
2. The first collection contains those parts of the polygons that do not overlap with others
3. The overlapping parts are kept in the second collection.
4. Finally, we add the superimposed part to the smaller of the two polygons concerned by the superimposition.

Unlike the SQL window in pgAdmin and PostgreSQL/Postgis, in the DB Manager SQL window you need to execute query by query.

To adapt the script to your situation, you need to replace the following elements with a word processor:

1. all tampon occurrences with the name of your table containing the buffered localities,
2. all occurrences of table_corrigee with the name of your desired result table,
3. and check that the identifier of your buffered table is id. By default, it can be named fid. In this case, replace all occurrences of id with fid.

Here are the script lines:

First query (if you’ve already run the script)

DROP TABLE IF EXISTS table_corrigee ; 

Second request (if you’ve already run the script)

DROP VIEW auto_jointure;

Third request

CREATE OR REPLACE VIEW auto_jointure AS SELECT

    tampon.id as id,
    tampon.nom as "nom",
    tampon.commune as "commune",
-- retrieve all polygons from the buffer table in two collections: one for complete geometries and another for overlapping parts. For the moment, both collections contain the complete geometries.
    ST_Union(tampon.geom) AS full_geom,
    ST_Union(tampon_bis.geom) AS shared_geom
    FROM  tampon ,tampon AS tampon_bis
    WHERE
      --check polygon validity as a precaution
      ST_IsValid(tampon.geom) AND ST_IsValid(tampon_bis.geom)
      --filter to retain intersecting polygons
      AND ST_intersects(tampon.geom,tampon_bis.geom)   
      --eliminate self intersecting polygons
      AND tampon.id <> tampon_bis.id          
      --for an intersection of 2 polygons, keep only the smaller one
      AND ST_Area(tampon_bis.geom) < ST_Area(tampon.geom)
      --since we're making "unions", we need to perform a grouping on the other attributes
    GROUP BY tampon.id,tampon."nom" , tampon.commune ;

Fourth request

/*
We're going to create a table that will contain the final result: the first step is to remove all overlapping areas. The second step is to add the overlapping areas to the smaller of the two polygons. The holes created by the first step are filled by the second step. The final polygons are seamless, with no overlaps or gaps.
*/

    CREATE TABLE table_corrigee AS SELECT
    id,
    "nom",
    commune,
/* Intersections are subtracted from the layer containing all polygons, thus removing conflicting partss*/
    ST_Multi(ST_Difference(full_geom,shared_geom)) as geom,
    ST_Area(full_geom) as area
    FROM auto_jointure
    WHERE ST_IsValid(full_geom) AND ST_IsValid(shared_geom)
/*The intersections just subtracted must then be added to fill the spaces created*/
    UNION 
      SELECT
      id,
      "nom",
      commune,
      geom,
      ST_Area(geom)
      FROM tampon
      WHERE id NOT IN (SELECT id FROM auto_jointure);

DB Manager troubleshooting

After running your queries, you may be surprised to find that neither the stable nor the views created during scripting appear in the list of Geopackage tables.

And yet, if you type a query referring to the script tables, they work. The tables are somewhere, but DB Manager doesn’t see them… It’s a bug, known and recognized, and…accepted! As the QGis developers are preparing a replacement for DB Manager, Data Source Manager, they have no intention of solving this problem.

In short, if we can’t solve the problem, we can work around it.

With DB Manager, we’ll create a table with the standard option

And we create a table with the same fields as the buffered layer

This one appears in the list of geopackage table.

Now all we have to do is load the result of our script into this “official” table with the following SQL query:

insert into recup
select fid,geom,nom,commune from lieux_dits_corr

We load this new layer into QGis:

We’ve solved the problem of topological consistency for localities.

For the purists out there, a final explanation of the table display problem in DB Manager.

The problem doesn’t come from QGIS, but from the geopackage standard.
OGR (and qgis) only look for tables in the gpkg_contents table. To make them visible, you need to add a new line to these tables:

But you’d expect QGIS DB Manager to execute these queries automatically.

Si cet article vous a intéressé et que vous pensez qu'il pourrait bénéficier à d'autres personnes, n'hésitez pas à le partager sur vos réseaux sociaux en utilisant les boutons ci-dessous. Votre partage est apprécié !

Atilio Francois