SpatiaLite Cookbook

	Recette #17 Railways vs Populated Places
Février 2011

Table des matières

Le problème

ici, on utilisera également la table railways.
Il s'agit par contre d'une recette très épicée: nous passons à la vitesse supérieure.
Comme vous vous en doutez, calculer des distances entre un rail et des Populated Places n'est pas si difficile.
On va donc introduire un peu de complexité au problème.

Imaginons la classification suivante:

Class	Min. distance	Max. distance
A-class	0 Km	1 Km
B-class	1 Km	2.5 Km
C-class	2.5 Km	5 Km
D-class	5 Km	10 Km
E-class	10 Km	20 Km

Le problème que vous devez résoudre est le suivant:

identifier chaque Populated Place situé à moins de 20 Km d'un Railway.
identifier la classe de distance pour chacune de ces Populated Places.

Railway	PopulatedPlace	A class [< 1Km]	B class [< 2.5Km]	C class [< 5Km]	D class [< 10Km]	E class [< 20Km]
Ferrovia Adriatica	Zapponeta	NULL	NULL	NULL	NULL	1
Ferrovia Adriatica	Villamagna	NULL	NULL	NULL	NULL	1
Ferrovia Adriatica	Villalfonsina	NULL	NULL	NULL	1	0
Ferrovia Adriatica	Vasto	1	0	0	0	0
...	...	...	...	...	...	...

SELECT rw.name AS Railway,
  pp_e.name AS PopulatedPlace,
  (ST_Distance(rw.geometry,
    Transform(pp_a.geometry, 23032)) <= 1000.0)
      AS "A class [< 1Km]",
  (ST_Distance(rw.geometry,
    Transform(pp_b.geometry, 23032)) > 1000.0)
      AS "B class [< 2.5Km]",
  (ST_Distance(rw.geometry,
    Transform(pp_c.geometry, 23032)) > 2500.0)
      AS "C class [< 5Km]",
  (ST_Distance(rw.geometry,
    Transform(pp_d.geometry, 23032)) > 5000.0)
      AS "D class [< 10Km]",
  (ST_Distance(rw.geometry,
    Transform(pp_e.geometry, 23032)) > 10000.0)
      AS "E class [< 20Km]"
FROM railways AS rw
JOIN populated_places AS pp_e ON (
  ST_Distance(rw.geometry,
    Transform(pp_e.geometry, 23032)) <= 20000.0
  AND pp_e.id IN (
    SELECT pkid
    FROM idx_populated_places_geometry
    WHERE pkid MATCH RTreeIntersects(
      MbrMinX(
        Transform(
          ST_Envelope(rw.geometry), 4326)),
      MbrMinY(
        Transform(
          ST_Envelope(rw.geometry), 4326)),
      MbrMaxX(
        Transform(
          ST_Envelope(rw.geometry), 4236)),
      MbrMaxY(
        Transform(
          ST_Envelope(rw.geometry), 4326)))))
LEFT JOIN populated_places AS pp_d ON (
  pp_e.id = pp_d.id
  AND ST_Distance(rw.geometry,
    Transform(pp_d.geometry, 23032)) <= 10000.0
  AND pp_d.id IN (
    SELECT pkid
    FROM idx_populated_places_geometry
    WHERE pkid MATCH RTreeIntersects(
      MbrMinX(
        Transform(
          ST_Envelope(rw.geometry), 4326)),
      MbrMinY(
        Transform(
          ST_Envelope(rw.geometry), 4326)),
      MbrMaxX(
        Transform(
          ST_Envelope(rw.geometry), 4236)),
      MbrMaxY(
        Transform(
          ST_Envelope(rw.geometry), 4326)))))
LEFT JOIN populated_places AS pp_c ON (
  pp_d.id = pp_c.id
  AND ST_Distance(rw.geometry,
    Transform(pp_c.geometry, 23032)) <= 5000.0
  AND pp_c.id IN (
    SELECT pkid
    FROM idx_populated_places_geometry
    WHERE pkid MATCH RTreeIntersects(
      MbrMinX(
        Transform(
          ST_Envelope(rw.geometry), 4326)),
      MbrMinY(
        Transform(
          ST_Envelope(rw.geometry), 4326)),
      MbrMaxX(
        Transform(
          ST_Envelope(rw.geometry), 4236)),
      MbrMaxY(
        Transform(
          ST_Envelope(rw.geometry), 4326)))))
LEFT JOIN populated_places AS pp_b ON (
  pp_c.id = pp_b.id
  AND ST_Distance(rw.geometry,
    Transform(pp_b.geometry, 23032)) <= 2500.0
  AND pp_b.id IN (
    SELECT pkid
    FROM idx_populated_places_geometry
    WHERE pkid MATCH RTreeIntersects(
      MbrMinX(
        Transform(
          ST_Envelope(rw.geometry), 4326)),
      MbrMinY(
        Transform(
          ST_Envelope(rw.geometry), 4326)),
      MbrMaxX(
        Transform(
          ST_Envelope(rw.geometry), 4236)),
      MbrMaxY(
        Transform(
          ST_Envelope(rw.geometry), 4326)))))
LEFT JOIN populated_places AS pp_a ON (
  pp_b.id = pp_a.id
  AND ST_Distance(rw.geometry,
    Transform(pp_a.geometry, 23032)) <= 1000.0
  AND pp_a.id IN (
    SELECT pkid
    FROM idx_populated_places_geometry
    WHERE pkid MATCH RTreeIntersects(
      MbrMinX(
        Transform(
          ST_Envelope(rw.geometry), 4326)),
      MbrMinY(
        Transform(
          ST_Envelope(rw.geometry), 4326)),
      MbrMaxX(
        Transform(
          ST_Envelope(rw.geometry), 4236)),
      MbrMaxY(
        Transform(
          ST_Envelope(rw.geometry), 4326)))));

Oui, cette requête est complexe et intimidante.
Cependant, en y regardant de plus près, elle reste abordable.
Vous connaissez déjà l'astuce: disséquer la requête pièce par pièce.

Examinons le squelette de la requète:

SELECT rw.name AS Railway, ...
FROM railways AS rw
JOIN populated_places AS pp_e ON (...)
LEFT JOIN populated_places AS pp_d ON (...)
LEFT JOIN populated_places AS pp_c ON (...)
LEFT JOIN populated_places AS pp_b ON (...)
LEFT JOIN populated_places AS pp_a ON (...);

on réalise simplement une JOINture des tables railways AS rw et populated_places AS pp_e
puis on réalise un LEFT JOIN de la table populated_places AS pp_d (2ème instance de la table): comme vous le savez LEFT JOIN insère une ligne dans les résultats même si le terme de droite est évalué NULL.
enfin, on répète le LEFT JOIN avec pp_c, pp_b et pp_a.
On va commencer par évaluer la plus grande distance et on continue selon les distances décroissantes.

...
JOIN populated_places AS pp_e ON (
ST_Distance(rw.geometry,
Transform(pp_e.geometry, 23032)) <= 20000.0
...

Chaque JOIN (ou LEFT JOIN) évalue simplement la distance entre Railway et Populated Place, et vérifie qu'elle appartient bien à la classe correspondante
la fonction Transform() est indispensable car les SRID ne sont pas identiques pour les deux jeux de donnée.

... AND pp_e.id IN (
  SELECT pkid
  FROM idx_populated_places_geometry
  WHERE pkid MATCH RTreeIntersects(
    MbrMinX(
      Transform(
        ST_Envelope(rw.geometry), 4326)),
    MbrMinY(
      Transform(
        ST_Envelope(rw.geometry), 4326)),
    MbrMaxX(
      Transform(
        ST_Envelope(rw.geometry), 4236)),
    MbrMaxY(
      Transform(
        ST_Envelope(rw.geometry), 4326)))
...

l'utilisation de l'Index Spatial est également nécessaire afin de filtrer rapidement les géométries des Populated Places.
ici encore, Transform()est indispensable pour reprojeter la géométrie de railway.
note: transformer les coordonnées d'un POINT peut être considéré comme une opération rapide;
mais transformer (des centaines et des centaines de fois …) des LINESTRING ou POLYGON est beaucoup plus lent.
On utilisera ainsi la fonction ST_Envelope() afin de simplifier au maximum la géométrie avant la transformation des coordonnées.

Bien, la structure de la requête est maintenant claire:

le premier JOIN inclus dans les résultats les Populated Places situés à moins de 20 Km des lignes de train(railway).
les autres LEFT JOIN vont tester les autres classes de distances.
et chaque LEFT JOIN vérifie scrupuleusement si l'ID des Populated Places est identique à celui de la classe précédente, comme avec: pp_d.id = pp_c.id
chaque fois qu'un LEFT JOIN échouera,la valeur NULL correspondante sera insérée dans les résultats.

SELECT rw.name AS Railway,
  pp_e.name AS PopulatedPlace,
  (ST_Distance(rw.geometry,
    Transform(pp_a.geometry, 23032)) <= 1000.0)
      AS "A class [< 1Km]",
...

Juste un dernier point a expliquer:

considérons une distance e.g. 3.8 Km: elle est incluse dans la classe C-class. mais cette condition satisfait le LEFT JOIN pour les classes D et E également .
il est donc nécessaire d'effectuer une vérification supplémentaire.
fort heureusement SQL est un langage très intelligent. Toute colonne incluse dans les résultats peut représenter n'importe quelle expression voulue.
En SQL une expression logique est évaluée comme:
- 0 [FAUX]
- 1 [TVRAI]
- ou NULL, si un des opérateurs évalués est NULL.
et c'est tout.

A vous de jouer, et de pousser plus loin cette analyse.
i.e vous pouvez rajouter des clauses ORDER BY ou WHERE etc... : c'est facile désormais, n'est-ce pas?

Table des matières

	Author: Alessandro Furieri a.furieri@lqt.it Traduced From English By RIVIERE Romain
	This work is licensed under the Attribution-ShareAlike 3.0 Unported (CC BY-SA 3.0) license.

	Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.3 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts.