Bases de données spatiales et PostGIS

Bases de données spatiales et PostGIS#

Hide code cell source

 
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import matplotlib.colors as mcolors
import numpy as np
import pandas as pd
import seaborn as sns
from shapely.geometry import Point, LineString, Polygon, MultiPolygon
from shapely.ops import unary_union
import shapely

sns.set_theme(style="whitegrid", palette="muted", font_scale=1.1)

Note : Les requêtes SQL PostGIS de ce chapitre supposent une connexion à un serveur PostgreSQL avec l’extension PostGIS installée. Ils sont présentés comme blocs illustratifs. Les cellules Python exécutables utilisent la bibliothèque shapely pour les calculs et visualisations géométriques.

Données géographiques : les types fondamentaux#

Les données spatiales représentent des entités du monde réel par leur forme géométrique et leur position dans l’espace.

Définition 111

Les trois types géométriques fondamentaux en SQL spatial sont :

  • Point : une paire de coordonnées (x, y) ou (longitude, latitude). Représente un lieu précis : une ville, un capteur, une adresse.

  • LineString : une séquence ordonnée de points reliés par des segments. Représente une route, une rivière, un tracé de frontière.

  • Polygon : une région fermée délimitée par un anneau extérieur et éventuellement des anneaux intérieurs (trous). Représente un pays, un quartier, un bâtiment.

Ces types de base sont combinés en types multiples : MultiPoint, MultiLineString, MultiPolygon, et GeometryCollection.

WKT et WKB#

Définition 112

Deux formats standards d’échange des géométries :

  • WKT (Well-Known Text) : représentation textuelle lisible par un humain. Exemples : POINT(2.3522 48.8566), LINESTRING(0 0, 1 1, 2 0), POLYGON((0 0, 4 0, 4 4, 0 4, 0 0)).

  • WKB (Well-Known Binary) : représentation binaire compacte, utilisée pour le stockage et les échanges efficaces entre systèmes. Non lisible directement, mais interopérable entre tous les SIG.

Systèmes de coordonnées : SRID et projections#

Définition 113

Un système de référence spatiale (SRS) définit comment les coordonnées d’une géométrie correspondent à des positions sur la Terre. Il est identifié par un SRID (Spatial Reference ID), un entier issu du registre EPSG.

Les deux SRID les plus courants :

  • EPSG:4326 (WGS 84) : système géographique utilisé par le GPS. Les coordonnées sont en degrés décimaux (longitude, latitude). C’est le standard pour stocker des données géographiques mondiales.

  • EPSG:3857 (Web Mercator) : projection utilisée par Google Maps, OpenStreetMap. Les coordonnées sont en mètres. Pratique pour l’affichage, mais distord les surfaces aux hautes latitudes.

Remarque 71

Le choix de la projection est crucial pour les calculs :

  • Calculer une distance en degrés (EPSG:4326) sans reprojection donne un résultat incorrect car les degrés n’ont pas la même longueur selon la latitude.

  • Pour des calculs de distances métriques exacts, il faut soit utiliser le type GEOGRAPHY de PostGIS (calculs sur sphéroïde), soit reprojecter les données dans un système métrique local.

PostGIS : installation et types#

PostGIS est l’extension officielle de PostgreSQL pour les données spatiales. Elle ajoute des types, des fonctions et des index spéciaux.

Définition 114

PostGIS ajoute deux types principaux à PostgreSQL :

  • GEOMETRY : géométrie plane (coordonnées cartésiennes). Les calculs sont euclidiens. Plus rapide mais moins précis pour les longues distances sur la Terre.

  • GEOGRAPHY : géométrie géodésique (coordonnées sphériques, défaut EPSG:4326). Les calculs (distance, aire) tiennent compte de la courbure de la Terre. Plus précis pour les données mondiales.

Installation et initialisation (blocs illustratifs — requiert PostGIS) :

-- Activer l'extension
CREATE EXTENSION IF NOT EXISTS postgis;

-- Vérifier la version
SELECT PostGIS_Version();

-- Créer une table spatiale
CREATE TABLE villes (
    id       SERIAL PRIMARY KEY,
    nom      TEXT NOT NULL,
    pays     TEXT,
    geom     GEOMETRY(POINT, 4326)   -- Point en WGS 84
);

-- Insertion via WKT et ST_GeomFromText
INSERT INTO villes (nom, pays, geom) VALUES
  ('Paris',    'France',   ST_GeomFromText('POINT(2.3522 48.8566)',  4326)),
  ('Lyon',     'France',   ST_GeomFromText('POINT(4.8357 45.7640)',  4326)),
  ('Marseille','France',   ST_GeomFromText('POINT(5.3698 43.2965)',  4326)),
  ('Genève',   'Suisse',   ST_GeomFromText('POINT(6.1432 46.2044)',  4326)),
  ('Madrid',   'Espagne',  ST_GeomFromText('POINT(-3.7038 40.4168)', 4326));

-- Vue de la géométrie en WKT
SELECT nom, ST_AsText(geom) FROM villes;

Fonctions spatiales essentielles#

Définition 115

Les fonctions spatiales PostGIS les plus utilisées :

Fonction

Description

ST_Distance(a, b)

Distance entre deux géométries

ST_Contains(a, b)

a contient-il complètement b ?

ST_Intersects(a, b)

Les géométries se croisent-elles ?

ST_Area(g)

Aire d’un polygone

ST_Buffer(g, r)

Tampon de rayon r autour de g

ST_Transform(g, srid)

Reprojeter dans un autre système

ST_Union(a, b)

Union de deux géométries

ST_Centroid(g)

Centre géométrique

ST_Length(g)

Longueur d’une ligne

ST_Within(a, b)

a est-il dans b ?

Requêtes spatiales illustratives (PostGIS)#

-- Distance entre Paris et Lyon (en degrés, GEOMETRY)
SELECT ST_Distance(
    ST_GeomFromText('POINT(2.3522 48.8566)', 4326),
    ST_GeomFromText('POINT(4.8357 45.7640)', 4326)
) AS distance_degres;

-- Distance réelle en mètres avec le type GEOGRAPHY
SELECT ST_Distance(
    ST_GeographyFromText('POINT(2.3522 48.8566)'),
    ST_GeographyFromText('POINT(4.8357 45.7640)')
) AS distance_metres;

-- Villes dans un rayon de 500 km autour de Paris
SELECT v.nom,
       ROUND(ST_Distance(
           v.geom::geography,
           ST_GeographyFromText('POINT(2.3522 48.8566)')
       ) / 1000) AS distance_km
FROM villes v
WHERE ST_DWithin(
    v.geom::geography,
    ST_GeographyFromText('POINT(2.3522 48.8566)'),
    500000  -- 500 km en mètres
)
ORDER BY distance_km;

-- Tampon de 50 km autour de Paris (en projection métrique)
SELECT ST_Buffer(
    ST_Transform(ST_GeomFromText('POINT(2.3522 48.8566)', 4326), 2154),
    50000   -- 50 km en mètres (EPSG:2154 = Lambert 93, projection française)
) AS zone_tampon;

Index spatial GiST#

Définition 116

PostGIS utilise un index GiST (Generalized Search Tree) basé sur un R-tree pour accélérer les requêtes spatiales. Un R-tree décompose l’espace en rectangles minimaux englobants (MBR — Minimum Bounding Rectangle). Sans index, toute requête spatiale nécessite un parcours complet de la table (seq scan).

-- Créer un index GiST sur la colonne géométrique
CREATE INDEX idx_villes_geom ON villes USING GIST (geom);

-- L'index est utilisé automatiquement par les opérateurs &&, @>, etc.
-- Les fonctions ST_DWithin, ST_Intersects, ST_Contains bénéficient de l'index.
EXPLAIN ANALYZE
SELECT nom FROM villes
WHERE  ST_DWithin(geom::geography,
                  ST_GeographyFromText('POINT(2.3522 48.8566)'),
                  300000);

Remarque 72

Le R-tree supporte aussi la recherche du plus proche voisin (k-NN search) efficacement :

-- Les 3 villes les plus proches de Paris
SELECT nom,
       ST_Distance(geom::geography,
                   ST_GeographyFromText('POINT(2.3522 48.8566)')) / 1000
          AS distance_km
FROM  villes
ORDER BY geom::geography <->
         ST_GeographyFromText('POINT(2.3522 48.8566)')
LIMIT 3;

L’opérateur <-> calcule la distance en utilisant l’index GiST — c’est beaucoup plus rapide que ORDER BY ST_Distance(...) sans index.

Requêtes spatiales avancées#

Exemple 39

Trouver les points dans un polygone — cas d’usage typique : quels sites sont dans une zone administrative ?

-- Villes dans un polygone définissant approximativement la France métropolitaine
SELECT nom
FROM   villes
WHERE  ST_Contains(
    ST_GeomFromText('POLYGON((
        -5 42, 9 42, 9 51, -5 51, -5 42
    ))', 4326),
    geom
);

Exemple 40

Trouver les routes qui intersectent une zone :

CREATE TABLE routes (
    id    SERIAL PRIMARY KEY,
    nom   TEXT,
    geom  GEOMETRY(LINESTRING, 4326)
);

-- Routes qui traversent le département de l'Hérault
SELECT r.nom
FROM   routes r
JOIN   departements d ON ST_Intersects(r.geom, d.geom)
WHERE  d.code = '34';

Calculs avec shapely (Python exécutable)#

# Création de géométries avec shapely
paris     = Point(2.3522, 48.8566)
lyon      = Point(4.8357, 45.7640)
marseille = Point(5.3698, 43.2965)
geneve    = Point(6.1432, 46.2044)
madrid    = Point(-3.7038, 40.4168)

villes = {
    "Paris":     paris,
    "Lyon":      lyon,
    "Marseille": marseille,
    "Genève":    geneve,
    "Madrid":    madrid,
}

# Distance euclidienne (en degrés) — illustrative
print("Distances depuis Paris (degrés, euclidien) :")
for nom, pt in villes.items():
    if nom != "Paris":
        d = paris.distance(pt)
        print(f"  Paris → {nom:<12} : {d:.4f}°")

Distances depuis Paris (degrés, euclidien) :
  Paris → Lyon         : 3.9664°
  Paris → Marseille    : 6.3262°
  Paris → Genève       : 4.6266°
  Paris → Madrid       : 10.3878°

# Conversion approx. degrés → km (1° ≈ 111 km)
print("Distances approx. depuis Paris (km) :")
for nom, pt in villes.items():
    if nom != "Paris":
        d_deg = paris.distance(pt)
        d_km  = d_deg * 111
        print(f"  Paris → {nom:<12} : {d_km:.0f} km")

Distances approx. depuis Paris (km) :
  Paris → Lyon         : 440 km
  Paris → Marseille    : 702 km
  Paris → Genève       : 514 km
  Paris → Madrid       : 1153 km

# Création d'une zone tampon (buffer) autour de Paris
# En coordonnées géographiques, 1° ≈ 111 km
rayon_3deg = 3.0  # ~333 km

zone_paris = paris.buffer(rayon_3deg)

# Villes dans la zone
print("Villes dans un rayon de ~333 km autour de Paris :")
for nom, pt in villes.items():
    if zone_paris.contains(pt):
        print(f"  ✓ {nom}")
    else:
        print(f"  ✗ {nom}")

Villes dans un rayon de ~333 km autour de Paris :
  ✓ Paris
  ✗ Lyon
  ✗ Marseille
  ✗ Genève
  ✗ Madrid

# Polygone représentant approximativement la France métropolitaine
france_approx = Polygon([
    (-5, 42), (9, 42), (9, 51), (-5, 51), (-5, 42)
])

print("Villes dans l'approximation de la France métropolitaine :")
for nom, pt in villes.items():
    print(f"  {'✓' if france_approx.contains(pt) else '✗'} {nom}")

Villes dans l'approximation de la France métropolitaine :
  ✓ Paris
  ✓ Lyon
  ✓ Marseille
  ✓ Genève
  ✗ Madrid

# Route entre Paris et Marseille
route_pm = LineString([paris, lyon, marseille])
print(f"Longueur de la route Paris→Lyon→Marseille : {route_pm.length:.4f}° "
      f"(≈ {route_pm.length * 111:.0f} km)")

# Intersection route et zone tampon
intersection = route_pm.intersection(zone_paris)
print(f"Longueur de la route dans la zone Paris : {intersection.length * 111:.0f} km")

Longueur de la route Paris→Lyon→Marseille : 6.4910° (≈ 721 km)
Longueur de la route dans la zone Paris : 333 km

Visualisation : géométries avec matplotlib#

Hide code cell source

 
palette = sns.color_palette("muted", 6)
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# --- Carte simple avec les villes ---
ax = axes[0]
ax.set_title("Villes et zone tampon (shapely)", fontsize=12, fontweight='bold')

# Dessiner la zone tampon
tampon_xy = np.array(zone_paris.exterior.coords)
ax.fill(tampon_xy[:, 0], tampon_xy[:, 1], alpha=0.15, color=palette[0],
        label="Zone ~333 km autour de Paris")
ax.plot(tampon_xy[:, 0], tampon_xy[:, 1], color=palette[0], lw=1.5)

# France approximative
france_xy = np.array(france_approx.exterior.coords)
ax.plot(france_xy[:, 0], france_xy[:, 1], color=palette[3], lw=2,
        linestyle='--', label="Approx. France")

# Route
route_xy = np.array(route_pm.coords)
ax.plot(route_xy[:, 0], route_xy[:, 1], color=palette[4], lw=2.5,
        label="Route Paris→Lyon→Marseille", zorder=4)

# Points des villes
couleurs = {
    "Paris": palette[3], "Lyon": palette[2],
    "Marseille": palette[2], "Genève": palette[1], "Madrid": palette[5]
}
for nom, pt in villes.items():
    c = couleurs[nom]
    ax.scatter(pt.x, pt.y, s=80, color=c, zorder=5, edgecolors='white', lw=1.2)
    ax.annotate(nom, (pt.x, pt.y), textcoords="offset points",
                xytext=(6, 4), fontsize=9, fontweight='bold', color=c)

ax.set_xlabel("Longitude (°)")
ax.set_ylabel("Latitude (°)")
ax.set_xlim(-7, 11)
ax.set_ylim(38, 54)
ax.legend(fontsize=8, loc='upper left')
ax.set_aspect('equal')

# --- Opérations booléennes sur polygones ---
ax2 = axes[1]
ax2.set_title("Opérations géométriques", fontsize=12, fontweight='bold')

poly1 = Polygon([(0,0),(4,0),(4,3),(0,3),(0,0)])
poly2 = Polygon([(2,1),(6,1),(6,4),(2,4),(2,1)])
inter = poly1.intersection(poly2)
union_p = poly1.union(poly2)
diff  = poly1.difference(poly2)

def draw_poly(ax, poly, color, alpha, label, edgecolor='white'):
    if poly.is_empty:
        return
    coords = np.array(poly.exterior.coords) if hasattr(poly, 'exterior') else None
    if coords is not None:
        ax.fill(coords[:,0], coords[:,1], alpha=alpha, color=color, label=label)
        ax.plot(coords[:,0], coords[:,1], color=edgecolor, lw=1.2)

draw_poly(ax2, poly1, palette[0], 0.25, "Polygone A")
draw_poly(ax2, poly2, palette[2], 0.25, "Polygone B")
draw_poly(ax2, inter, palette[3], 0.60, "Intersection A∩B")

# Dessiner la différence A\B
diff_xy = np.array(diff.exterior.coords)
ax2.fill(diff_xy[:,0], diff_xy[:,1], alpha=0.45, color=palette[4],
         label="Différence A\\B")

# Centroïdes
for poly, label, c in [(poly1,"centroïde A",palette[0]),
                        (poly2,"centroïde B",palette[2])]:
    cx, cy = poly.centroid.x, poly.centroid.y
    ax2.scatter(cx, cy, s=60, color=c, zorder=5, marker='x', linewidths=2)
    ax2.annotate(label, (cx, cy), xytext=(5, 3),
                 textcoords='offset points', fontsize=8, color=c)

ax2.set_xlim(-1, 7)
ax2.set_ylim(-1, 5.5)
ax2.legend(fontsize=8.5, loc='upper left')
ax2.set_xlabel("x")
ax2.set_ylabel("y")
ax2.set_aspect('equal')

plt.savefig("_build_postgis_geom.png", dpi=120, bbox_inches='tight')
plt.show()
_images/90f3d20f7c924aa5260deffbc9bdd9196b7cb9dff3186b8f23faed88b7584519.png 
# Tableau récapitulatif des opérations spatiales
poly_a = Polygon([(0,0),(4,0),(4,3),(0,3),(0,0)])
poly_b = Polygon([(2,1),(6,1),(6,4),(2,4),(2,1)])

resultats = {
    "ST_Area(A)":              poly_a.area,
    "ST_Area(B)":              poly_b.area,
    "ST_Area(A ∩ B)":          poly_a.intersection(poly_b).area,
    "ST_Area(A ∪ B)":          poly_a.union(poly_b).area,
    "ST_Length(A → centre B)": LineString([poly_a.centroid, poly_b.centroid]).length,
    "ST_Intersects(A, B)":     poly_a.intersects(poly_b),
    "ST_Contains(A, B)":       poly_a.contains(poly_b),
    "ST_Distance(A, B)":       poly_a.distance(poly_b),
}
pd.DataFrame.from_dict(resultats, orient='index', columns=['Valeur'])
Valeur
ST_Area(A) 12.0
ST_Area(B) 12.0
ST_Area(A ∩ B) 4.0
ST_Area(A ∪ B) 20.0
ST_Length(A → centre B) 2.236068
ST_Intersects(A, B) True
ST_Contains(A, B) False
ST_Distance(A, B) 0.0

Cas d’usage : villes dans un rayon, intersections#

Exemple 41

Un cas d’usage courant est la zone d’isochrone : délimiter la zone accessible en moins de N minutes depuis un point. Les outils comme Valhalla ou OSRM calculent ces polygones, qui sont ensuite stockés dans PostGIS pour des requêtes du type « quels clients habitent à moins de 30 minutes de notre magasin ? ».

-- Supposons que la table isochrones contienne des polygones préalculés
SELECT c.nom, c.adresse
FROM   clients c
JOIN   isochrones i ON ST_Within(c.geom, i.geom)
WHERE  i.magasin_id = 42
  AND  i.minutes    = 30;

Remarque 73

Applications typiques des bases de données spatiales :

  • Urbanisme : trouver tous les bâtiments dans une zone inondable.

  • Logistique : optimiser les tournées de livraison, calculer les zones de chalandise.

  • Environnement : superposer des données climatiques et des aires protégées.

  • Télécommunications : calculer les zones de couverture des antennes.

  • Assurance : évaluer l’exposition au risque selon la localisation.

Résumé#

Remarque 74

Ce chapitre a introduit les bases de données spatiales et PostGIS :

Fondements :

  • Trois types géométriques fondamentaux : POINT, LINESTRING, POLYGON, et leurs variantes multiples.

  • WKT (texte lisible) et WKB (binaire) sont les formats d’échange standards.

  • Le SRID identifie le système de référence spatiale (EPSG:4326 = WGS 84 pour le GPS).

PostGIS :

  • GEOMETRY pour les calculs plans, GEOGRAPHY pour les calculs géodésiques précis.

  • Les fonctions spatiales (ST_Distance, ST_Contains, ST_Intersects, ST_Buffer, ST_Transform) couvrent la quasi-totalité des besoins analytiques.

  • L’index GiST (R-tree) est indispensable pour les performances sur les grandes tables spatiales.

  • L’opérateur <-> permet la recherche efficace du plus proche voisin via l’index.

Python avec shapely :

  • shapely permet de créer et manipuler des géométries sans serveur de base de données.

  • Les opérations booléennes (intersection, union, difference), les mesures (area, length, distance) et les prédicats (contains, intersects) sont directement disponibles.

Contenu