SpatiaLite logo

Recipe #17
Railways vs Populated Places

2011 January 28

Previous Slide Table of Contents Next Slide

The problem

We'll use once again the railways dataset.
But this really is an hot spiced recipe: be prepared to taste very strong flavors.
As you can now easily image by yourself, computing distances between a railway line and Populated Places isn't so difficult.
So this problem introduces a further degree of complexity (just to escape from boredom and to keep your mind active and interested).

Image that for any good reason the following classification exists:
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
The problem you are faced to resolve is:
  • identify any Populated Place laying within a distance radius of 20 Km from a Railway.
  • identify the corresponding distance Class for each one of such 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)))));
Yes, this one really looks like a complex and intimidating query.
Anyway, complexity is much more apparent than real.
You already know the trick: you simply have to break down this statement into several smallest chunks. And then you'll soon discover that there isn't nothing really difficult and complex.

Let us examine the main framework:
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 (...);
...
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)))
...
All right, now the main framework of the complex query is absolutely clear:
  • the first JOIN will include into the result-set any Populated Place falling within 20 Km from the railway line.
  • any other LEFT JOIN will then test decreasing distances, accordingly to the imposed Class boundaries.
  • and each LEFT JOIN carefully checks if the Populated Place ID is the same of the previous successfully identified Class, as in: pp_d.id = pp_c.id
  • each time one such LEFT JOIN will fail, then corresponding NULL-values will be inserted into the result-set.


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]",
...
Just a latest element to be shortly explained:
You can now play by yourself, performing further tests on this query.
i.e you can add some smart ORDER BY or WHERE clause and so on: that's really easy now, isn't ?

Previous Slide Table of Contents Next Slide

CC-BY-SA logo Author: Alessandro Furieri a.furieri@lqt.it
This work is licensed under the Attribution-ShareAlike 3.0 Unported (CC BY-SA 3.0) license.

GNU logo 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.