Chapter 5. Requêtes spatiales

Table of Contents

La raison d'être des bases de données spatiales est de réaliser à l'intérieur de la base de données des requêtes qui normalement nécessiteraient des fonctionnalités d'un SIG Desktop. Utiliser PostGIS requiert en effet la connaissance de quelles fonctions spatiales sont disponibles, comment les utiliser dans une requête, et de s'assurer de la mise en place adéquate des index, pour une bonne performance.

5.1. Déterminer les relations spatiales

Les relations spatiales indiquent comment deux géométries interagissent l'une avec l'autre. Elles constituent une fonctionnalité fondamentale pour effectuer des requêtes sur des géométries.

5.1.1. Modèle à 9 intersections dimensionnellement étendu

Selon le OpenGIS Simple Features Implementation Specification for SQL, "l'approche de base pour comparer deux géométries consiste à effectuer des tests par paire des intersections entre les intérieurs, les limites et les extérieurs des deux géométries et à classer la relation entre les deux géométries sur la base des entrées de la matrice "intersection" résultante."

Dans la théorie de la topologie des ensembles de points, les points d'une géométrie intégrée dans un espace à deux dimensions sont classés en trois ensembles :

Limites

La limite d'une géométrie est l'ensemble des géométries de la dimension immédiatement inférieure. Pour les POINT, qui ont une dimension de 0, la limite est l'ensemble vide. La limite d'une LINESTRING est constituée par ses deux extrémités. Pour les POLYGON, la limite est le réseau de lignes des anneaux extérieur et intérieur.

Intérieur

L'intérieur d'une géométrie est constitué des points de la géométrie qui ne se trouvent pas dans la limite. Pour les POINT, l'intérieur est le point lui-même. L'intérieur d'une LINESTRING est l'ensemble des points situés entre les extrémités. Pour les POLYGON, l'intérieur est la surface aréale à l'intérieur du polygone.

Extérieur

L'extérieur d'une géométrie est le reste de l'espace dans lequel la géométrie est intégrée ; en d'autres termes, tous les points qui ne se trouvent pas à l'intérieur ou sur la limite de la géométrie. Il s'agit d'une surface non fermée à 2 dimensions.

Le Dimensionally Extended 9-Intersection Model (DE-9IM) décrit la relation spatiale entre deux géométries en spécifiant les dimensions des 9 intersections entre les ensembles ci-dessus pour chaque géométrie. Les dimensions des intersections peuvent être représentées formellement dans une matrice d'intersections de 3x3.

Pour une géométrie g, les Limites, Intérieures, et Extérieures sont désignés par la notation I(g), B(g), et E(g). En outre, dim(s) désigne la dimension d'un ensemble s dont le domaine est {0,1,2,F} :

  • 0 => point

  • 1 => ligne

  • 2 => area

  • F => ensemble vide

En utilisant cette notation, la matrice d'intersection pour deux géométries a et b est :

  Intérieur Limite Extérieur
Intérieur dim( I(a) ∩ I(b) ) dim( I(a) ∩ B(b) ) dim( I(a) ∩ E(b) )
Limite dim( B(a) ∩ I(b) ) dim( B(a) ∩ B(b) ) dim( B(a) ∩ E(b) )
Extérieur dim( E(a) ∩ I(b) ) dim( E(a) ∩ B(b) ) dim( E(a) ∩ E(b) )

Visuellement, pour deux géométries polygonales qui se chevauchent, cela ressemble à ce qui suit :

 
  Intérieur Limite Extérieur
Intérieur

dim( I(a) ∩ I(b) ) = 2

dim( I(a) ∩ B(b) = 1

dim( I(a) ∩ E(b) ) = 2

Limite

dim( B(a) ∩ I(b) ) = 1

dim( B(a) ∩ B(b) ) = 0

dim( B(a) ∩ E(b) ) = 1

Extérieur

dim( E(a) ∩ I(b) ) = 2

dim( E(a) ∩ B(b) ) = 1

dim( E(a) ∩ E(b) = 2

En lisant de gauche à droite et de haut en bas, la matrice d'intersection est représentée par la chaîne de texte "212101212".

Pour plus d'informations, voir :

5.1.2. Relations spatiales nommées

Pour faciliter la détermination des relations spatiales communes, l'OGC SFS définit un ensemble de prédicats de relations spatiales. PostGIS fournit ces prédicats sous la forme des fonctions ST_Contains, ST_Crosses, ST_Disjoint, ST_Equals, ST_Intersects, ST_Overlaps, ST_Touches, ST_Within. Il définit également les prédicats de relation non standard ST_Covers, ST_CoveredBy, et ST_ContainsProperly.

Les prédicats spatiaux sont généralement utilisés comme conditions dans les clauses SQL WHERE ou JOIN. Les prédicats spatiaux nommés utilisent automatiquement un index spatial s'il existe, de sorte qu'il n'est pas nécessaire d'utiliser également l'opérateur de boîte de délimitation &&. Par exemple :

SELECT city.name, state.name, city.geom
FROM city JOIN state ON ST_Intersects(city.geom, state.geom);

Pour plus de détails et d'illustrations, voir l' Atelier PostGIS.

5.1.3. Relations spatiales générales

Dans certains cas, les relations spatiales nommées sont insuffisantes pour fournir une condition de filtrage spatial souhaitée.

Prenons l'exemple d'un ensemble de données linéaires représentant un réseau routier. Il peut être nécessaire d'identifier tous les segments de route qui se croisent, non pas en un point, mais en une ligne (peut-être pour valider une règle commerciale). Dans ce cas, ST_Crosses ne fournit pas le filtre spatial nécessaire, puisque pour les caractéristiques linéaires, il renvoie true uniquement lorsqu'elles se croisent en un point.

Une solution en deux étapes consisterait à calculer d'abord l'intersection réelle (ST_Intersection) des paires de lignes routières qui se croisent dans l'espace (ST_Intersects), puis vérifier si le ST_GeometryType de l'intersection est "LINESTRING" (en traitant correctement les cas qui renvoient des GEOMETRYCOLLECTIONs de [MULTI]POINTs, [MULTI]LINESTRINGs, etc. ).

Il est évident qu'une solution plus simple et plus rapide est souhaitable.

Un deuxième exemple est la localisation des quais qui coupent les limites d'un lac sur une ligne et dont l'une des extrémités se trouve sur le rivage. En d'autres termes, lorsqu'un quai est situé à l'intérieur d'un lac mais n'est pas complètement contenu par celui-ci, qu'il coupe la limite d'un lac sur une ligne et que l'une des extrémités du quai se trouve à l'intérieur ou sur la limite du lac. Il est possible d'utiliser une combinaison de prédicats spatiaux pour trouver les caractéristiques requises :

Ces exigences peuvent être satisfaites en calculant la matrice d'intersection DE-9IM complète. PostGIS fournit la fonction ST_Relate pour ce faire :

SELECT ST_Relate( 'LINESTRING (1 1, 5 5)',
                  'POLYGON ((3 3, 3 7, 7 7, 7 3, 3 3))' );
st_relate
-----------
1010F0212

Pour tester une relation spatiale particulière, un modèle de matrice d'intersection est utilisé. Il s'agit de la représentation matricielle augmentée des symboles supplémentaires {T,*} :

  • T => la dimension d'intersection est non vide ; c'est-à-dire qu'elle se trouve dans {0,1,2}

  • * => indifférent

Les modèles de matrice d'intersection permettent d'évaluer des relations spatiales spécifiques de manière plus succincte. Les fonctions ST_Relate et ST_RelateMatch peuvent être utilisées pour tester les modèles de matrice d'intersection. Pour le premier exemple ci-dessus, le modèle de matrice d'intersection spécifiant deux lignes se croisant dans une ligne est '1*1***1**' :

-- Find road segments that intersect in a line
SELECT a.id
FROM roads a, roads b
WHERE a.id != b.id
      AND a.geom && b.geom
      AND ST_Relate(a.geom, b.geom, '1*1***1**');

Pour le deuxième exemple, le modèle de matrice d'intersection spécifiant une ligne située en partie à l'intérieur et en partie à l'extérieur d'un polygone est "102101FF2" :

-- Find wharves partly on a lake's shoreline
SELECT a.lake_id, b.wharf_id
FROM lakes a, wharfs b
WHERE a.geom && b.geom
      AND ST_Relate(a.geom, b.geom, '102101FF2');

5.2. Utilisation des index spatiaux

Lors de la construction de requêtes utilisant des conditions spatiales, il est important, pour de meilleures performances, de s'assurer qu'un index spatial est utilisé, s'il existe (voir Section 4.9, “Index spatiaux”). Pour ce faire, un opérateur spatial ou une fonction sensible à l'index doit être utilisé dans une clause WHERE ou ON de la requête.

Les opérateurs spatiaux comprennent les opérateurs de boîte de délimitation (dont le plus couramment utilisé est && ; voir Section 7.10.1, “Opérateurs de Bounding Box” pour la liste complète) et les opérateurs de distance utilisés dans les requêtes sur les plus proches voisins (le plus courant étant <-> ; voir Section 7.10.2, “Opérateurs de distance” pour la liste complète.)

Les fonctions tenant compte de l'index ajoutent automatiquement un opérateur de boîte de délimitation à la condition spatiale. Les fonctions tenant compte de l'index comprennent les prédicats de relation spatiale nommés ST_Contains, ST_ContainsProperly, ST_CoveredBy, ST_Covers, ST_Crosses, ST_Intersects, ST_Overlaps, ST_Touches, ST_Within, ST_Within, et ST_3DIntersects, et les prédicats de distance ST_DWithin, ST_DFullyWithin, ST_3DDFullyWithin, et ST_3DDWithin . )

Les fonctions telles que ST_Distance n'utilisent pas les index pour optimiser leur fonctionnement. Par exemple, la requête suivante serait assez lente sur une grande table :

SELECT geom
FROM geom_table
WHERE ST_Distance( geom, 'SRID=312;POINT(100000 200000)' ) < 100

Cette requête sélectionne toutes les géométries de la geom_table qui se trouvent à moins de 100 unités du point (100000, 200000). Elle sera lente car elle calcule la distance entre chaque point du tableau et le point spécifié, c'est-à-dire qu'un calcul ST_Distance() est effectué pour chaque ligne du tableau.

Le nombre de lignes traitées peut être considérablement réduit en utilisant la fonction d'indexation ST_DWithin :

SELECT geom
FROM geom_table
WHERE ST_DWithin( geom, 'SRID=312;POINT(100000 200000)', 100 )

Cette requête sélectionne les mêmes géométries, mais de manière plus efficace. Ceci est possible en ST_DWithin() utilisant l'opérateur && en interne sur une boîte de délimitation élargie de la géométrie de la requête. S'il existe un index spatial sur geom, le planificateur de requêtes reconnaîtra qu'il peut utiliser l'index pour réduire le nombre de lignes analysées avant de calculer la distance. L'index spatial permet d'extraire uniquement les enregistrements dont les boîtes de délimitation chevauchent l'étendue élargie et qui, par conséquent, pourraient se trouver à l'intérieur de la distance requise. La distance réelle est ensuite calculée pour confirmer l'inclusion de l'enregistrement dans l'ensemble des résultats.

Pour plus d'informations et d'exemples, voir l' atelier de travail PostGIS.

5.3. Exemples de SQL spatial

Les exemples de cette section utilisent une table de routes linéaires et une table de limites de communes polygonales. La définition de la table bc_roads est la suivante :

Column    | Type              | Description
----------+-------------------+-------------------
gid       | integer           | Unique ID
name      | character varying | Road Name
geom      | geometry          | Location Geometry (Linestring)

La définition de la table bc_municipality est la suivante :

Column   | Type              | Description
---------+-------------------+-------------------
gid      | integer           | Unique ID
code     | integer           | Unique ID
name     | character varying | City / Town Name
geom     | geometry          | Location Geometry (Polygon)

5.3.1.

Quelle est la longueur totale de toutes les routes, exprimée en kilomètres ?

Vous pouvez répondre à cette question à l'aide d'un code SQL très simple :

SELECT sum(ST_Length(geom))/1000 AS km_roads FROM bc_roads;

km_roads
------------------
70842.1243039643

5.3.2.

Quelle est la superficie de la ville de Prince George, en hectares ?

Cette requête combine une condition d'attribut (sur le nom de la municipalité) avec un calcul spatial (de la surface du polygone) :

SELECT
  ST_Area(geom)/10000 AS hectares
FROM bc_municipality
WHERE name = 'PRINCE GEORGE';

hectares
------------------
32657.9103824927

5.3.3.

Quelle est la plus grande municipalité de la province, en termes de superficie ?

Cette requête utilise une mesure spatiale comme valeur de référence. Il existe plusieurs façons d'aborder ce problème, mais la plus efficace est la suivante :

SELECT
  name,
  ST_Area(geom)/10000 AS hectares
FROM bc_municipality
ORDER BY hectares DESC
LIMIT 1;

name           | hectares
---------------+-----------------
TUMBLER RIDGE  | 155020.02556131

Notez que pour répondre à cette question, nous devons calculer la surface de chaque polygone. Si nous faisions cela souvent, il serait logique d'ajouter une colonne de surface à la table qui pourrait être indexée pour des raisons de performance. En ordonnant les résultats dans une direction décroissante et en utilisant la commande "LIMIT" de PostgreSQL, nous pouvons facilement sélectionner la plus grande valeur sans utiliser une fonction d'agrégation comme MAX().

5.3.4.

Quelle est la longueur des routes entièrement comprises dans chaque municipalité ?

Il s'agit d'un exemple de "jointure spatiale", qui rassemble les données de deux tables (avec une jointure) en utilisant une interaction spatiale ("contenu") comme condition de jointure (plutôt que l'approche relationnelle habituelle de jointure sur une clé commune) :

SELECT
  m.name,
  sum(ST_Length(r.geom))/1000 as roads_km
FROM bc_roads AS r
JOIN bc_municipality AS m
  ON ST_Contains(m.geom, r.geom)
GROUP BY m.name
ORDER BY roads_km;

name                        | roads_km
----------------------------+------------------
SURREY                      | 1539.47553551242
VANCOUVER                   | 1450.33093486576
LANGLEY DISTRICT            | 833.793392535662
BURNABY                     | 773.769091404338
PRINCE GEORGE               | 694.37554369147
...

Cette requête prend un certain temps, car chaque route du tableau est résumée dans le résultat final (environ 250 000 routes pour le tableau de l'exemple). Pour des ensembles de données plus petits (plusieurs milliers d'enregistrements sur plusieurs centaines), la réponse peut être très rapide.

5.3.5.

Créez un nouveau tableau avec toutes les routes de la ville de Prince George.

Il s'agit d'un exemple de "superposition", qui prend en compte deux tableaux et produit un nouveau tableau composé de résultats coupés dans l'espace. Contrairement à la "jointure spatiale" présentée ci-dessus, cette requête crée de nouvelles géométries. Une superposition est comme une jointure spatiale turbocompressée, et est utile pour des travaux d'analyse plus précis :

CREATE TABLE pg_roads as
SELECT
  ST_Intersection(r.geom, m.geom) AS intersection_geom,
  ST_Length(r.geom) AS rd_orig_length,
  r.*
FROM bc_roads AS r
JOIN bc_municipality AS m
  ON ST_Intersects(r.geom, m.geom)
WHERE
  m.name = 'PRINCE GEORGE';

5.3.6.

Quelle est la longueur en kilomètres de "Douglas St" à Victoria ?

SELECT
  sum(ST_Length(r.geom))/1000 AS kilometers
FROM bc_roads r
JOIN bc_municipality m
  ON ST_Intersects(m.geom, r.geom
WHERE
  r.name = 'Douglas St'
  AND m.name = 'VICTORIA';

kilometers
------------------
4.89151904172838

5.3.7.

Quel est le plus grand polygone municipal comportant un trou ?

SELECT gid, name, ST_Area(geom) AS area
FROM bc_municipality
WHERE ST_NRings(geom) 
> 1
ORDER BY area DESC LIMIT 1;

gid  | name         | area
-----+--------------+------------------
12   | SPALLUMCHEEN | 257374619.430216