---
jupytext:
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.0
kernelspec:
  name: python3
  display_name: Python 3
---

# Bases de données spatiales et PostGIS

```{code-cell} python
:tags: [hide-input]

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.

```{prf:definition}
:label: ch14-def-geometries

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

```{prf:definition}
:label: ch14-def-wkt-wkb

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

```{prf:definition}
:label: ch14-def-srid

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.
```

```{prf:remark}
:label: ch14-rem-projections

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.

```{prf:definition}
:label: ch14-def-postgis-types

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) :

```sql
-- 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

```{prf:definition}
:label: ch14-def-fonctions-spatiales

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)

```sql
-- 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

```{prf:definition}
:label: ch14-def-gist

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*).
```

```sql
-- 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);
```

```{prf:remark}
:label: ch14-rem-gist-knn

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

```sql
-- 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

```{prf:example}
:label: ch14-ex-points-dans-polygone

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

```sql
-- 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
);
```

```{prf:example}
:label: ch14-ex-routes-zone

Trouver les routes qui intersectent une zone :

```sql
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)

```{code-cell} python
# 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}°")
```

```{code-cell} python
# 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")
```

```{code-cell} python
# 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}")
```

```{code-cell} python
# 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}")
```

```{code-cell} python
# 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")
```

## Visualisation : géométries avec matplotlib

```{code-cell} python
:tags: [hide-input]

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()
```

```{code-cell} python
# 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'])
```

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

```{prf:example}
:label: ch14-ex-isochrone

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 ?".

```sql
-- 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;
```

```{prf:remark}
:label: ch14-rem-cas-usage

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é

```{prf:remark}
:label: ch14-rem-synthese

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.
```