IX. Analyse spatiale

IX.5 Un exemple d'application : dites-le avec du SQL !

• Création d'une nouvelle base SpatiaLite ou PostGIS
  • Création d'une base SpatiaLite
  • Création d'une base PostGIS
• Import de données dans SpatiaLite ou PostGIS
• Lancer une requête simple
• Création d'une grille
  • Créer une grille avec SpatiaLite
  • Créer une grille avec PostGIS
• Union et agrégation
• Évolution temporelle : soustraction de 2 maillages
  • Relancer l'opération sur les données Corine Land Cover 2012
  • Soustraire les 2 maillages 2012 et 2000

Comment faire pour automatiser les opérations réalisées au chapitre précédent, afin de pouvoir rendre le processus plus reproductible ?

Une première solution serait d'utiliser un modèle. Nous ne verrons pas ici le pas à pas pour créer le modèle en question, mais vous pouvez essayer vous-même en vous référant ici !

Une autre solution pour automatiser ce traitement est d'utiliser des requêtes SQL. Il s'agit d'une partie « pour aller plus loin », d'un niveau plus avancé, et vous pouvez très bien décider de vous arrêter ici. Nous nous appuierons sur cette partie.

Il est possible d'utiliser 2 logiciels différents à partir de QGIS pour lancer des requêtes SQL :

• SQLite et son module spatial SpatiaLite ne nécessitent pas d'installation, mais les fonctions disponibles sont plus limitées que celles de Postgresql/PostGIS
• Postgresql et son module spatial PostGIS doivent être installés en plus de QGIS, mais une fois cette opération réalisée beaucoup de possibilités s'offriront à vous !

Il n'y a pas une solution meilleure qu'une autre mais elles répondent à des besoins différents.

Pour cet exercice les 2 logiciels peuvent être utilisés. Si vous choisissez Postgresql/PostGIS, l'installation ne sera pas détaillée ici mais on trouve de nombreuses ressources sur internet. La syntaxe étant légèrement différente d'un logiciel à l'autre, les requêtes seront proposées en 2 versions.

Pour ce chapitre, vous pouvez soit télécharger les données Corine Land Cover : Données Métropole 2000 et Données Métropole 2012 puis les filtrer pour ne garder que les vignes, comme détaillé dans le tutoriel, ou bien utiliser les données en téléchargement déjà filtrées (pour un téléchargement moins lourd).

Création d'une nouvelle base SpatiaLite ou PostGIS

La première étape consiste à créer une base de données vides dans le logiciel choisi, et à y importer les données de départ, à savoir la couche Corine Land Cover.

Si vous choisissez d'utiliser Postgresql/PostGIS, assurez-vous d'avoir installé ces logiciels avant de poursuivre !

Création d'une base SpatiaLite

Nous allons tout d'abord créer une base SpatiaLite vide.

Si vous naviguez dans l'explorateur de fichiers de votre ordinateur jusqu'à votre base de données, vous verrez qu'un fichier SQLite a été créé. Par rapport à un format comme le shapefile, une base de données est comme une boîte qui peut contenir plusieurs jeux de données.

Une base SpatiaLite est une base SQLite qui peut gérer des données spatiales grâce au module SpatiaLite. Elle est constituée d'une seule fichier qu'on peut copier d'un dossier à l'autre.

Ça n'est pas le cas des logiciels de bases de données « traditionnels » tels que PostgreSQL, qui fonctionnent selon une logique client-serveur : la base est installé sur un serveur, et des utilisateurs peuvent s'y connecter à partir d'autres ordinateurs. Bien sûr il est possible de se servir de son propre ordinateur comme de serveur, et la base ne sera accessible que depuis celui-ci : c'est ce qu'on appelle une base « en local ».

Création d'une base PostGIS

...Et c'est ce que nous allons faire ici ! Si vous avez décidé de travailler avec SpatiaLite, passez directement à l'import des données.

Je pars du principe ici que vous avez déjà installé PostgreSQL et PostGIS, et que vous êtes connecté à un serveur PostgreSQL, local ou distant.

QGIS peut uniquement se connecter à une base déjà existante, la création d'une nouvelle base PostgreSQL/PostGIS ne peut se faire dans QGIS.

Il existe plusieurs méthodes pour cela, qui sortent un peu de ce tutoriel : cette partie consistera surtout en des renvois vers des tutoriels existants.

PostgreSQL en soi ne possède pas d'interface graphique mais fonctionne en lignes de commande. Si ce mot vous fait peur, pas de panique ! Il existe plusieurs logiciels qui peuvent servir d'interface graphique à PostgreSQL et vous permettre d'arriver à vos fins en quelques clics.

Le plus connu est peut-être pgAdmin, qui est peut-être déjà installé sur votre ordinateur selon la manière dont vous avez installé PostgreSQL. Vous pouvez également utiliser par exemple DBeaver.

Pour créer une nouvelle base de données PostgreSQL, en la nommant maillage_clc :

• en ligne de commande ou ici
• avec pgAdmin en déroulant pour aller à la partie « 3) Creating a new database using pgAdmin »
• avec DBeaver, le principe est le même qu'avec pgAdmin : clic droit sur le serveur postgres → Créer → Database

Il faut ensuite donner les 2 extensions postgis et postgis_topology à la base pour pouvoir y mettre des données géographiques. Encore une fois, vous pouvez procéder au choix :

• en ligne de commande avec ces 2 commandes : CREATE EXTENSION postgis; et CREATE EXTENSION postgis_topology;
• avec pgAdmin, en faisant un clic droit sur la base → Créer → Extension puis en choisissant postgis et postgis_topology
• avec DBeaver, sensiblement de la même manière qu'avec pgAdmin

Votre base PostGIS est créée, il ne reste plus qu'à s'y connecter dans QGIS !

Une dernière étape consiste à créer un schéma nommé tutoqgis dans cette base. Un schéma correspond à un sous-dossier dans la base de données. Le schéma par défaut se nomme public mais c'est une bonne pratique de créer vos propres schémas.

Vous remarquerez que cette étape est plus longue et complexe avec PostgreSQL. Néanmoins, une fois ce logiciel configuré et avec un peu d'habitude, il ne vous faudra plus que quelque secondes pour créer une nouvelle base !

Import de données dans SpatiaLite ou PostGIS

Cette étape est presque équivalente pour les 2 logiciels, en passant par le gestionnaire de base de données de QGIS.

Avant d'importer les données, nous allons sélectionner dans la couche Corine Land Cover les entités correspondant au vignoble, pour n'importer que celles-ci. Cette étape pourrait tout à fait être effectuée une fois les données importées dans la base, mais comme cela le temps d'import dans SpatiaLite ou PostGIS sera plus court.

Il ne reste plus qu'à importer les entités sélectionnes, dans SpatiaLite ou PostGIS selon votre choix.

Bravo ! Vous disposez maintenant d'une base de données avec une couche de données.

Nous allons pouvoir rentrer dans le vif du sujet !

Lancer une requête simple

A ce stade, si vous n'avez aucune notion de SQL, je vous conseille de d'abord suivre cette partie sur les requêtes SQL.

Nous allons ici lancer une requête simple dans SpatiaLite ou PostGIS, en guise d'introduction.

Remarquez que cet onglet porte le nom de votre base : il ne sera pas possible d'y lancer une requête concernant une autre base (mais on peut très bien ouvrir plusieurs onglets de requête).

L'onglet comporte 2 parties :

1. en haut, vous pouvez taper une requête
2. en bas, vous pourrez visualiser le résultat de cette requête.

Pour tester cet outil, nous allons sélectionner les polygones ayant une surface inférieure à 25 hectares.

Les différents mots clés (select, from, where...) peuvent être écrits indifféremment en minuscules ou majuscules. De même, la requête peut être écrite en une seule ligne ou bien avec des retours à la ligne.

Si l'on détaille cette requête ligne par ligne :

SELECT *

signifie que nous allons sélectionner (select) toutes (la mention *) les colonnes de la table attributaire, ainsi que la géométrie, qui est considérée comme une colonne nommée geom, comme vous pouvez le vérifier dans l'onglet Info.

FROM clc00_vignes

signifie que nous allons sélectionner les colonnes de la couche clc00_vignes.

WHERE area_ha < 25

applique un critère à la requête : seules seront sélectionnées les lignes répondant à ce critère, c'est-à-dire dont la valeur pour le champ area_ha est inférieure à 25.

Ici, la requête ne crée pas de nouvelles couches mais renvoie les lignes sélectionnées. Comment faire pour créer une nouvelle couche à partir de cette sélection ?

Pour SpatiaLite, les choses sont se compliquent un petit peu : la même requête renvoie bien une table avec les mêmes colonnes, mais la colonne geom n'est pas reconnue comme colonne de géométrie. Il s'agit d'une simple table sans composante géographique.

Vous pouvez d'ailleurs le voir dans l'illustration ci-dessus : l'icône de inf25ha dans la base SpatiaLite est celle d'une table, alors que c'est un polygone dans la base PostGIS.

Pour expliquer cette requête : RecoverGeometryColumn est la fonction permettant de transformer une colonne ordinaire en colonne de géométrie (sous réserve bien sûr que son contenu corresponde bien à des géométries). Cette fonction prend plusieurs arguments :

• inf25ha est le nom de la couche
• geom est le nom de la colonne de géométrie
• 2154 est le code EPSG du SCR de la couche
• 'MULTIPOLYGON' est le type de géométries de la couche (ici, c'est le même type que la couche clc00_vignes qui a servi a créer inf25ha)
• 'XY' correspond à la dimension de la couche

Une autre différence entre PostGIS et SpatiaLite : dans PostGIS, il est possible de lancer plusieurs requêtes à la suite les unes des autres, pourvu que chaque requête se termine par un point-virgule. Dans SpatiaLite, il faut lancer les requêtes une à une, en cliquant sur le bouton Exécuter entre chaque.

Nous allons maintenant utiliser des requêtes plus complexes pour créer une grille !

Création d'une grille

Notre première étape consiste à créer une grille ayant la même étendue que notre couche clc00_vignes, avec une maille de 50km. C'est l'équivalent de cette étape réalisée au chapitre précédent.

Il existe dans SpatiaLite une fonction spécifique pour créer une grille ; la fonction équivalente n'est accessible dans PostGIS qu'à partir de la version 3.1. A moins de disposer de cette version, il faudra donc utiliser une fonction « fait maison ».

Créer une grille avec SpatiaLite

Dans SpatiaLite, nous allons pouvoir utiliser la fonction ST_SquareGrid. 4 étapes seront nécessaires :

1. Création d'une table vide
2. Ajout d'une colonne de géométrie à cette table
3. Mise à jour de cette colonne pour créer une grille avec une seule entité multiparties
4. Séparation de l'entité multiparties pour que 1 case = 1 polygone

Créer une grille avec PostGIS

Malheureusement, la fonction ST_SquareGrid permettant la génération d'une grille avec PostGIS n'est accessible qu'à partir de la version 3.1. A moins de disposer de cette version, il faudra donc utiliser notre propre fonction !

Une fonction est un bout de code pouvant être « appelé ». C'est en quelque sorte un raccourci qui permet d'éviter de taper une série d'instructions, en tapant seulement le nom de cette série d'instructions.

Les fonctions peuvent prendre des arguments en entrée : par exemple, une couche, une taille de maille... Et peuvent renvoyer un résultat un sortie, par exemple une grille.

Nous allons ici utiliser cette fonction créée par Alexander Palamarchuk pour générer une grille.

La fonction makegrid_2d est maintenant accessible dans PostGIS : vous n'aurez plus besoin de retaper ce code.

Il ne reste plus qu'à appeler cette fonction avec en entrée :

• l'étendue de la grille, c'est-à-dire l'étendue de clc00_vignes
• la taille de maille, soit 50 000 mètres

Quelques explications, cette fois-ci nous prendrons cette requête non plus ligne par ligne, mais fonction par fonction...

Une petite astuce : en cliquant sur une parenthèse ouvrante ou fermante dans l'onglet de requête du gestionnaire de bases de données, cette parenthèse et son alter ego sont surlignées en vert, ce qui permet de mieux comprendre l'emboîtement des fonctions.

create table tutoqgis.grid00 as

Cette première ligne veut simplement dire que cette requête va créer une nouvelle table nommée grid00 dans le schéma tutoqgis.

SELECT row_number() over () as gid, geom FROM (SELECT...) AS q_grid

Cette nouvelle table contiendra 2 colonnes : une colonne gid d'identifiant unique, créée avec la fonction row_number(), et une colonne de géométrie geom. Comme la table est créée à partir d'un deuxième SELECT, il faut donner un nom (alias) à cette sous-requête, ici q_grid.

Vous pouvez essayer de relancer la requête en omettant la partie AS q_grid, vous obtiendrez un message d'erreur vous indiquant que la sous-requête doit avoir un alias : ERROR: subquery in FROM must have an alias.

(SELECT (ST_Dump(...)).geom AS geom

La sous-requête utilise la fonction ST_Dump, qui permet de créer des entités à une seule partie. Pour récupérer une géométrie en retour avec cette fonction, on a ajouté .geom, et pour nommer cette géométrie geom AS geom.

makegrid_2d(..., ...)

La fonction ST_Dump prend un seul paramètre en entrée correspondant à une géométrie. Ici, cette géométrie sera celle renvoyée en sortie par la fonction st_makegrid créée précédemment, qui prend elle en entrée 2 arguments, séparés par une virgule.

SELECT ST_SetSRID(ST_Extent(geom), 2154) FROM tutoqgis.clc00_vignes

Le premier argument de la fonction makegrid_2d correspond à une étendue. On utilise pour la créer la fonction ST_Extent sur la couche clc00_vignes. Il faut également attribuer un système de coordonnées (SRID dans la jargon PostGIS) à cette étendue, ce qui est fait avec la fonction St_SetSRID, qui utilise 2 paramètres : le résultat de ST_Extent, et le code EPSG 2154 (RGF93 Lambert 93).

Si vous testez la requête sans utiliser la fonction ST_SetSRID, en remplaçant la partie ci-dessus par « select st_extent(geom) from tutoqgis.clc00_vignes », vous obtiendrez un message d'erreur vous indiquant que la géométrie en entrée n'a pas de système de coordonnées : ERROR: ST_Transform: Input geometry has unknown (0) SRID

50000

50000 est le deuxième paramètre utilisé par la fonction makegrid_2d, qui correspond à la longueur du côté d'une maille dans le SCR utilisé (ici Lambert 93, donc en mètres).

A ce stade, si vous avez suivi le chapitre précédent et créé une grille avec l'outil Créer une grille de QGIS, l'opération paraît bien plus compliquée en SQL. Avec un peu de chance la partie suivante vous donnera l'impression inverse !

Union et agrégation

Nous allons maintenant donner à chaque case de cette grille une valeur correspondant à sa surface en vignes, à partir de la couche clc00_vignes.

Cette opération regroupe les 3 parties du chapitre précédent : union, recalcul de la surface et agrégation des données par maille (je vous avais bien dit que le SQL avait des avantages).

Comme d'habitude, à vous de choisir votre logiciel préféré pour cette opération, qui nécessite donc 2 couches : une grille, et une couche de polygones.

Comment fonctionne cette requête ? Prenons ses lignes une à une, mais dans le désordre (ici il s'agit de la requête PostGIS, mais les explications sont les mêmes pour SpatiaLite) :

CREATE TABLE tutoqgis.grid00_surf AS

Pas de souci ici, il s'agit de créer une table nommée grid00_surf dans le schéma tutoqgis. On passe 2 lignes plus loin !

FROM tutoqgis.grid00 AS g, tutoqgis.clc00_vignes AS v

Comme prévu, cette requête va utiliser les 2 couches grid00 et clc00_vignes, toutes 2 dans le schéma tutoqgis. Dans le reste de la requête, on se référera à ces tables par leur alias : g pour la grille, v pour les vignes. On remonte à la ligne du dessus...

Si on veut relancer cette requête sur d'autres couches, grâce à l'alias il n'y aura besoin de modifier les noms des couches qu'à un seul endroit !

SELECT g.gid, g.geom, sum(ST_Area(ST_Intersection(v.geom, g.geom)))/10000 AS surf

La couche créée va comporter 3 colonnes :

• La colonne d'identifiant unique gid
• Une colonne de géométrie, identique à celle de la couche grid00
• Une colonne nommée surf, correspondant à la somme (fonction SUM) des surfaces du résultat de l'intersection entre la grille et les vignes (fonction ST_Intersection). On divise cette somme par 10 000 pour obtenir non plus des surfaces en mètres mais en hectares.

On descend 2 lignes plus bas...

GROUP BY g.gid, g.geom

La clause GROUP BY permet de regrouper toutes les entités ayant la même valeur de géométrie et d'identifant. Elle implique d'utiliser une fonction d'agrégation pour les colonnes autres que geom et gid. Ici, c'est la fonction d'agrégation SUM qui est utilisée dans la deuxième ligne de la requête pour créer la colonne surf.

La requête pourrait être lancée telle quelle. Pour qu'elle soit moins longue à s'exécuter, on a rajouté la ligne

WHERE ST_Intersects(g.geom, v.geom)

qui permet de ne prendre en compte que les cases de la grille qui sont superposées avec des vignes. Avec ce critère, la table créée ne contient donc que ces cases. Vous pouvez tester en supprimant cette ligne, le résultat sera un peu plus long à créer.

ORDER BY g.gid

Enfin, cette dernière ligne, optionnelle, permet de choisir l'ordre des lignes dans la table, ici un ordre croissant sur le champ gid.

Si vous utilisez SpatiaLite, il reste une requête pour que la colonne géométrie soit reconnue en tant que telle.

On peut maintenant visualiser le résultat, que vous ayiez utilisé SpatiaLite ou PostGIS.

Avec 5 ou 6 lignes de SQL, vous avez accompli l'équivalent de 3 outils QGIS et de beaucoup de clics !

Et surtout, il sera très facile de relancer toute cette opération sur d'autres données, comme nous allons le faire ci-dessous.

Évolution temporelle : soustraction de 2 maillages

Relancer l'opération sur les données Corine Land Cover 2012

Nous avons jusqu'ici travaillé sur les données Corine Land Cover 2000. Nous allons maintenant utiliser les données équivalentes pour l'année 2012, ce qui nous permettra de visualiser l'évolution entre ces 2 années.

Nous allons donc relancer l'opération précédente (union et agrégation) sur la couche CLC 2012.

Rapide, non ? Quand vous travaillez avec des requêtes SQL, une pratique peut être de copier les requêtes dans un fichier texte. Ainsi vous en gardez la mémoire, et si vous avez besoin de les relancer sur d'autres données, il vous suffira de faire un rechercher-remplacer sur les noms des couches, puis de recopier ces requêtes dans le gestionnaire de bases de données QGIS.

Il est possible en SQL d'ajouter des commentaires, non pris en compte : la ligne doit alors être précédée par 2 tirets. Ceci vous permet d'expliquer vos requêtes, ce qui est toujours utile quand on reprend un travail quelques semaines/mois/années plus tard, ou pour vos collègues.

Soustraire les 2 maillages 2012 et 2000

En guise de dernière application pour ce chapitre, nous allons voir 2 manières de visualiser l'évolution de la surface en vignes entre 2000 et 2012 :

• Avec des requêtes SQL
• En mode raster

On pourrait aussi travailler en mode raster avec des requêtes SQL ! Mais ce chapitre est déjà bien rempli.

L'idée est donc de soustraire pour chaque case de grille les données 2000 aux données 2012, afin de visualiser une évolution, négative ou positive.

Si on utilise le langage SQL, on va pouvoir faire une jointure attributaire sur le champ gid entre les 2 couches grid00_surf et gri12_surf.

La requête sera la même pour SpatiaLite et PostGIS, il faut simplement ajouter le nom du schéma tutoqgis devant le nom des tables pour PostGIS, et mettre à jour la colonne de géométrie pour SpatiaLite.

Cette requête crée une couche evol_00_12 avec 3 colonnes :

• gid, ici en récupérant le champ gid de la couche grid00 (mais on aurait très bien pu remplacer g1.gid par g2.gid, le résultat serait le même)
• une colonne diff_surf correspondant à la différence de surface en vignes entre 2012 et 2000
• la géométrie geom, idem on aurait pu remplacer g1.geom par g2.geom

La ligne

WHERE g1.gid = g2.gid

est celle qui fait la jointure entre les 2 couches, sur le champ gid.

Une autre manière de faire pour cette opération est de convertir les 2 maillages en rasters, et de soustraire ces rasters l'un à l'autre.

Bravo, vous êtes arrivés au bout de ce chapitre ! Vous pouvez vous reposer avec le suivant sur les modes de représentation et la mise en page de cartes.

chapitre précédent partie X : représentation et mise en page
haut de page