PostGIS pour débuter : charger, indexer, joindre

Le GeoPackage est parfait pour un fichier qu'on transporte ; dès que plusieurs personnes travaillent sur les mêmes données, que les tables dépassent quelques millions d'objets ou qu'une application web doit les interroger, on passe à un vrai serveur de base de données : PostgreSQL et son extension spatiale PostGIS. Ce tutoriel couvre les premiers pas : installer, charger des données, les indexer et faire une jointure spatiale.

GeoPackage ou PostGIS ?

Installer

• Windows : installez PostgreSQL avec l'installeur d'EDB, puis ajoutez PostGIS avec l'outil Stack Builder proposé à la fin.
• Linux : sudo apt install postgis (Debian, Ubuntu) installe PostgreSQL et PostGIS.
• Docker, le plus rapide pour essayer :

docker run -d --name postgis -p 5432:5432 -e POSTGRES_PASSWORD=motdepasse postgis/postgis

Créez ensuite une base et activez l'extension (avec psql, pgAdmin ou DBeaver) :

CREATE DATABASE sig;
-- puis, connecté à la base « sig » :
CREATE EXTENSION postgis;
SELECT postgis_full_version();

Charger des données

La méthode universelle est ogr2ogr (voir les recettes ogr2ogr) : il lit tous les formats et crée la table, ses colonnes et son index spatial.

# Un GeoPackage de communes vers la table « communes », reprojeté en Lambert 93
ogr2ogr -f PostgreSQL PG:"host=localhost dbname=sig user=postgres password=motdepasse" communes.gpkg \
  -nln communes -t_srs EPSG:2154 -nlt PROMOTE_TO_MULTI -lco GEOMETRY_NAME=geom -lco FID=id

# Un CSV de points (colonnes lon et lat)
ogr2ogr -f PostgreSQL PG:"dbname=sig user=postgres" arrets.csv -nln arrets \
  -oo X_POSSIBLE_NAMES=lon -oo Y_POSSIBLE_NAMES=lat -a_srs EPSG:4326 -t_srs EPSG:2154 -lco GEOMETRY_NAME=geom

• -nlt PROMOTE_TO_MULTI évite l'erreur classique d'un fichier qui mélange Polygon et MultiPolygon.
• Sans -lco GEOMETRY_NAME=geom, la colonne s'appellera wkb_geometry.
• Autres voies : shp2pgsql fourni avec PostGIS, ou dans QGIS le Gestionnaire BD (bouton Importer) et l'algorithme « Exporter vers PostgreSQL ».

Le type geometry

Une couche PostGIS est une table ordinaire avec une colonne de type geometry, typée et rattachée à un système de coordonnées :

CREATE TABLE parcelles (
    id     serial PRIMARY KEY,
    idu    text,
    geom   geometry(MultiPolygon, 2154)
);
-- Les tables spatiales de la base et leur système
SELECT f_table_name, f_geometry_column, type, srid FROM geometry_columns;

Les calculs se font dans l'unité du système : en Lambert 93, ST_Area rend des m² et ST_Distance des mètres. Pour des données en WGS84 (degrés), utilisez ST_Transform(geom, 2154), ou le type geography, qui calcule sur l'ellipsoïde et rend des mètres.

Indexer : le réflexe à prendre

Un index spatial GiST range les objets par leur rectangle englobant. Une requête spatiale commence alors par trouver les candidats dont le rectangle touche la zone, ce qui est très rapide, avant le calcul exact sur ces seuls candidats. ogr2ogr crée cet index par défaut ; pour une table créée à la main :

CREATE INDEX parcelles_geom_idx ON parcelles USING GIST (geom);
ANALYZE parcelles;   -- met à jour les statistiques utilisées par le planificateur

EXPLAIN ANALYZE devant une requête montre si l'index est utilisé (Index Scan ou Bitmap Index Scan sur …_geom_idx) ou si la table est parcourue en entier (Seq Scan).

Requêtes spatiales

-- Surface de chaque commune en km²
SELECT nom, round((ST_Area(geom) / 1e6)::numeric, 2) AS km2 FROM communes ORDER BY km2 DESC;

-- Jointure spatiale : nombre d'arrêts par commune
SELECT c.nom, count(a.id) AS nb_arrets
FROM communes c
LEFT JOIN arrets a ON ST_Intersects(c.geom, a.geom)
GROUP BY c.nom
ORDER BY nb_arrets DESC;

-- Écoles à moins de 300 m d'une route principale
SELECT DISTINCT e.nom
FROM ecoles e
JOIN routes r ON ST_DWithin(e.geom, r.geom, 300)
WHERE r.importance <= 2;

-- Rattacher à chaque point la commune qui le contient (nouvelle table)
CREATE TABLE arrets_communes AS
SELECT a.*, c.nom AS commune
FROM arrets a JOIN communes c ON ST_Within(a.geom, c.geom);

Réparer et vérifier

-- Géométries invalides et raison
SELECT id, ST_IsValidReason(geom) FROM parcelles WHERE NOT ST_IsValid(geom);
-- Réparation (voir le tutoriel Topologie et géométries)
UPDATE parcelles SET geom = ST_Multi(ST_CollectionExtract(ST_MakeValid(geom), 3)) WHERE NOT ST_IsValid(geom);

Les causes et les réparations sont détaillées dans Topologie et géométries valides.

Utiliser PostGIS dans QGIS

1. Dans l'Explorateur, clic droit sur PostgreSQL › Nouvelle connexion : hôte, port 5432, base, identifiants (de préférence dans une configuration d'authentification plutôt qu'en clair dans le projet).
2. Les tables apparaissent ; un glisser-déposer les ajoute à la carte, éditables comme un fichier.
3. Le Gestionnaire BD exécute du SQL et charge le résultat comme une couche, à condition de renvoyer une colonne d'identifiant unique et une géométrie.
4. Une vue (CREATE VIEW) se charge comme une couche et reste toujours à jour.

Démonstration : pourquoi l'index change tout

Une jointure spatiale « points dans communes » calculée dans votre navigateur, sur les communes de la Métropole de Lyon et des points tirés au hasard. Sans index, chaque point est testé contre chaque commune. Avec l'index, on ne teste exactement que les communes dont le rectangle englobant contient le point : c'est le principe d'un index GiST.

À essayer : lancez la jointure sans index, puis avec, et comparez le nombre de tests exacts « point dans polygone ». Augmentez le nombre de points : l'écart grandit. Cochez « Rectangles englobants » pour voir ce que l'index a réellement en mémoire.