Kapitel 5. Räumliche Abfrage

Inhaltsverzeichnis

Der Sinn von räumlichen Datenbanken liegt darin Abfragen in der Datenbank ausführen zu können, die normalerweise Desktop-GIS-Funktionionalität verlangen würden. Um PostGIS effektiv verwenden zu können muss man die verfügbaren räumlichen Funktionen kennen, wissen wie sie in Abfragen verwendet werden und sicherstellen, dass für gute Performanz die passenden Indizes vorhanden sind.

5.1. Räumliche Beziehungen feststellen

Räumliche Beziehungen geben an wie zwei Geometrien miteinander interagieren. Sie sind die fundamentale Fähigkeit zum Abfragen von Geometrie.

5.1.1. Dimensionell erweitertes 9-Schnitte-Modell

Laut OpenGIS Simple Features Implementation Specification for SQL besteht der grundlegende Ansatz für den Vergleich zweier Geometrien darin, paarweise Tests der Schnittpunkte zwischen den Innen-, Rand- und Außenbereichen der beiden Geometrien durchzuführen und die Beziehung zwischen den beiden Geometrien auf der Grundlage der Einträge in der resultierenden 'Schnittpunkt'-Matrix zu klassifizieren".

In der Theorie der Punktmengentopologie werden die Punkte in einer Geometrie, die in den 2-dimensionalen Raum eingebettet ist, in drei Gruppen eingeteilt:

Grenze

Die Begrenzung einer Geometrie ist die Menge der Geometrien der nächstniedrigeren Dimension. Für POINTs, die eine Dimension von 0 haben, ist die Begrenzung die leere Menge. Die Begrenzung eines LINESTRING sind die beiden Endpunkte. Für POLYGONs ist die Begrenzung das Liniennetz der äußeren und inneren Ringe.

Innenbereich

Das Innere einer Geometrie sind die Punkte einer Geometrie, die nicht in der Begrenzung liegen. Für POINTs ist das Innere der Punkt selbst. Das Innere eines LINESTRING ist die Menge der Punkte zwischen den Endpunkten. Für POLYGONs ist das Innere die Fläche innerhalb des Polygons.

Äußeres

Das Äußere einer Geometrie ist der Rest des Raums, in den die Geometrie eingebettet ist; mit anderen Worten, alle Punkte, die sich nicht im Inneren oder auf dem Rand der Geometrie befinden. Es handelt sich um eine 2-dimensionale, nicht geschlossene Fläche.

Das Dimensionally Extended 9-Intersection Model (DE-9IM) beschreibt die räumliche Beziehung zwischen zwei Geometrien durch Angabe der Dimensionen der 9 Schnittpunkte zwischen den oben genannten Mengen für jede Geometrie. Die Schnittpunktdimensionen können formal in einer 3x3 Schnittpunktmatrix dargestellt werden.

Für eine Geometrie g werden Interior, Boundary, und Exterior mit den Bezeichnungen I(g), B(g), und E(g) bezeichnet. Außerdem bezeichnet dim(s) die Dimension einer Menge s mit dem Bereich {0,1,2,F}:

  • 0 => Punkt

  • 1 => Zeile

  • 2 => Fläche

  • F => leere Menge

Unter Verwendung dieser Notation lautet die Schnittpunktmatrix für zwei Geometrien a und b:

  Innenbereich Grenze Äußeres
Innenbereich dim( I(a) ∩ I(b) ) dim( I(a) ∩ B(b) ) dim( I(a) ∩ E(b) )
Grenze dim( B(a) ∩ I(b) ) dim( B(a) ∩ B(b) ) dim( B(a) ∩ E(b) )
Äußeres dim( E(a) ∩ I(b) ) dim( E(a) ∩ B(b) ) dim( E(a) ∩ E(b) )

The following overlapping polygons provide a concrete example. ST_Relate computes the matrix 212101212; it matches the overlap pattern T*T***T** used by ST_Overlaps.

Code
WITH example AS (
  SELECT
    'POLYGON ((140 140,140 122,135 105,126 100,118 99,110 94,100 86,97 73,102 59,98 49,87 38,70 30,55 29,40 30,28 38,20 50,14 66,10 84,6 100,4 119,4 143,6 166,10 180,18 191,28 195,40 190,55 189,67 186,86 179,112 165,124 158,133 148,140 140),(48 177,40 177,30 174,24 166,22 159,25 155,30 153,40 154,45 157,52 162,53 169,51 174,48 177))'::geometry AS a,
    'POLYGON ((75 63,79 50,87 38,95 31,108 25,124 22,140 18,154 11,166 6,176 10,184 21,188 35,190 58,190 82,193 104,190 121,185 139,178 154,166 163,154 171,139 172,124 171,112 165,96 152,92 142,92 126,86 116,79 110,75 104,72 94,73 86,75 76,75 63))'::geometry AS b
), relation AS (
  SELECT a, b, ST_Relate(a, b) AS matrix
  FROM example
)
SELECT matrix,
       ST_RelateMatch(matrix, 'T*T***T**') AS matches_overlap_pattern,
       ST_Overlaps(a, b) AS overlaps
FROM relation;
Ausgabe von Rastern
matrix   | matches_overlap_pattern | overlaps
-----------+-------------------------+----------
 212101212 | t                       | t
(1 row)
Figure
Geometry figure for visual-de9im-overlapping-polygons

The executable query and generated overlay above are the visual source for the relation. The query below renders the nine cells of the same kind of matrix as a 3-by-3 set of examples. The final exterior/exterior cell is clipped to a finite viewport, since the true exterior/exterior intersection is unbounded.

Code
WITH example AS (
  SELECT
    ST_MakeEnvelope(0, 0, 6, 4) AS a,
    ST_MakeEnvelope(3, -1, 8, 3) AS b,
    ST_MakeEnvelope(-1, -2, 9, 5) AS viewport
), parts AS (
  SELECT
    a,
    b,
    ST_Boundary(a) AS ba,
    ST_Boundary(b) AS bb,
    viewport
  FROM example
)
SELECT
  ST_Relate(a, b) AS matrix,
  ST_AsText(ST_Normalize(a)) AS input_a,
  ST_AsText(ST_Normalize(b)) AS input_b,
  ST_AsText(ST_Normalize(ST_Intersection(a, b))) AS "2 I(a) ∩ I(b)",
  ST_AsText(ST_Normalize(ST_Intersection(a, bb))) AS "1 I(a) ∩ B(b)",
  ST_AsText(ST_Normalize(ST_Difference(a, b))) AS "2 I(a) ∩ E(b)",
  ST_AsText(ST_Normalize(ST_Intersection(ba, b))) AS "1 B(a) ∩ I(b)",
  ST_AsText(ST_Normalize(ST_Intersection(ba, bb))) AS "0 B(a) ∩ B(b)",
  ST_AsText(ST_Normalize(ST_Collect(
      'LINESTRING(0 0,0 4,6 4,6 3)'::geometry,
      'LINESTRING(0 0,3 0)'::geometry
  ))) AS "1 B(a) ∩ E(b)",
  ST_AsText(ST_Normalize(ST_Difference(b, a))) AS "2 E(a) ∩ I(b)",
  ST_AsText(ST_Normalize(ST_Collect(
      'LINESTRING(3 -1,8 -1,8 3,6 3)'::geometry,
      'LINESTRING(3 -1,3 0)'::geometry
  ))) AS "1 E(a) ∩ B(b)",
  ST_AsText(ST_Normalize(ST_Difference(viewport, ST_Union(a, b))))
    AS "2 E(a) ∩ E(b) clipped"
FROM parts;
Ausgabe von Rastern
-[ RECORD 1 ]----------+-----------------------------------------------------------------------------
matrix                 | 212101212
input_a                | POLYGON((0 0,0 4,6 4,6 0,0 0))
input_b                | POLYGON((3 -1,3 3,8 3,8 -1,3 -1))
2 I(a) ∩ I(b)          | POLYGON((3 0,3 3,6 3,6 0,3 0))
1 I(a) ∩ B(b)          | LINESTRING(3 0,3 3,6 3)
2 I(a) ∩ E(b)          | POLYGON((0 0,0 4,6 4,6 3,3 3,3 0,0 0))
1 B(a) ∩ I(b)          | LINESTRING(3 0,6 0,6 3)
0 B(a) ∩ B(b)          | MULTIPOINT((6 3),(3 0))
1 B(a) ∩ E(b)          | MULTILINESTRING((0 0,0 4,6 4,6 3),(0 0,3 0))
2 E(a) ∩ I(b)          | POLYGON((3 -1,3 0,6 0,6 3,8 3,8 -1,3 -1))
1 E(a) ∩ B(b)          | MULTILINESTRING((3 -1,8 -1,8 3,6 3),(3 -1,3 0))
2 E(a) ∩ E(b) clipped  | POLYGON((-1 -2,-1 5,9 5,9 -2,-1 -2),(0 0,3 0,3 -1,8 -1,8 3,6 3,6 4,0 4,0 0))
Figure
Geometry figure for visual-de9im-matrix-cells

Von links nach rechts und von oben nach unten gelesen, wird die Kreuzungsmatrix als Textstring '212101212' dargestellt.

Weitere Informationen finden Sie unter:

5.1.2. Benannte räumliche Beziehungen

Um die Bestimmung allgemeiner räumlicher Beziehungen zu erleichtern, definiert das OGC SFS eine Reihe von genannten räumlichen Beziehungsprädikaten. PostGIS stellt diese als die Funktionen ST_Contains, ST_Crosses, ST_Disjoint, ST_Equals, ST_Intersects, ST_Overlaps, ST_Touches, ST_Within zur Verfügung. Es definiert auch die Nicht-Standard-Beziehungsprädikate ST_Covers, ST_CoveredBy und ST_ContainsProperly.

Räumliche Prädikate werden normalerweise als Bedingungen in SQL WHERE oder JOIN Klauseln verwendet. Die benannten räumlichen Prädikate verwenden automatisch einen räumlichen Index, wenn einer vorhanden ist, so dass es nicht notwendig ist, auch den Bounding-Box-Operator && zu verwenden. Zum Beispiel:

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

Weitere Einzelheiten und Abbildungen finden Sie im PostGIS Workshop.

5.1.3. Allgemeine räumliche Beziehungen

In some cases the named spatial relationships are insufficient to provide a desired spatial filter condition. These requirements can be expressed by computing the full DE-9IM intersection matrix with ST_Relate.

Um eine bestimmte räumliche Beziehung zu testen, wird ein Kreuzungsmatrixmuster verwendet. Dabei handelt es sich um die Matrixdarstellung, die um die zusätzlichen Symbole {T,*} erweitert wurde:

  • T => Schnittmenge ist nicht leer, d. h. sie liegt in {0,1,2}

  • * => ist mir egal

Using intersection matrix patterns, specific spatial relationships can be evaluated succinctly. The ST_Relate and ST_RelateMatch functions can both test these patterns.

These road segments share a line segment. ST_Crosses does not find this relationship for linear features, because it only returns true when their interiors intersect at a point. The pattern 1*1***1** tests for a one-dimensional intersection instead.

Code
WITH roads AS (
  SELECT
    'LINESTRING(10 10,40 90,70 110,140 110,170 130,190 190)'::geometry AS "road A",
    'LINESTRING(10 190,50 130,90 110,130 110,160 70,180 10)'::geometry AS "road B"
), relation AS (
  SELECT "road A", "road B", ST_Relate("road A", "road B") AS matrix
  FROM roads
)
SELECT
  ST_Intersection("road A", "road B") AS overlap,
  matrix,
  ST_Crosses("road A", "road B") AS "ST_Crosses",
  ST_RelateMatch(matrix, '1*1***1**') AS "matches 1*1***1**"
FROM relation;
Ausgabe von Rastern
overlap           |  matrix   | ST_Crosses | matches 1*1***1**
----------------------------+-----------+------------+-------------------
 LINESTRING(90 110,130 110) | 1F1FF0102 | f          | t
(1 row)
Figure
Geometry figure for visual-relate-overlapping-roads

The next example finds a wharf which runs from the lake interior onto its shoreline, with one endpoint on the shoreline. The other wharves provide nearby counterexamples in the figure. The pattern 102101FF2 expresses the complete required relationship in one test.

Code
WITH data AS (
  SELECT
    'POLYGON((-20 10,30 30,80 70,110 80,145 85,190 110,300 220,230 220,-20 220,-20 10))'::geometry AS lake,
    'MULTILINESTRING((120 180,145 85,110 80),(10 130,30 70),(90 140,110 80,94.81118881118876 75.83496503496498),(150 160,180 80))'::geometry AS wharves
), wharf AS (
  SELECT lake, (item).path[1] AS id, (item).geom
  FROM data CROSS JOIN LATERAL ST_Dump(wharves) AS item
), relation AS (
  SELECT id, geom, ST_Relate(lake, geom) AS matrix
  FROM wharf
)
SELECT
  id AS wharf,
  geom AS "matching wharf",
  matrix
FROM relation
WHERE ST_RelateMatch(matrix, '102101FF2')
ORDER BY id;
Ausgabe von Rastern
wharf |           matching wharf           |  matrix
-------+-----------------------------------+-----------
     1 | LINESTRING(120 180,145 85,110 80) | 102101FF2
(1 row)
Figure
Geometry figure for visual-relate-lake-wharves

5.2. Räumliche Indizes verwenden

Bei der Erstellung von Abfragen mit räumlichen Bedingungen ist es für eine optimale Leistung wichtig, dass ein räumlicher Index verwendet wird, sofern ein solcher existiert (siehe Abschnitt 4.9, „Räumliche Indizes“). Zu diesem Zweck muss ein räumlicher Operator oder eine indexfähige Funktion in einer WHERE oder ON Klausel der Abfrage verwendet werden.

Zu den räumlichen Operatoren gehören die Bounding-Box-Operatoren (von denen der am häufigsten verwendete && ist; eine vollständige Liste finden Sie unter Abschnitt 7.10.1, „Bounding-Box-Operatoren“ ) und die Abstandsoperatoren, die bei Abfragen nach dem nächsten Nachbarn verwendet werden (der am häufigsten verwendete ist <->; eine vollständige Liste finden Sie unter Abschnitt 7.10.2, „Operatoren“ ).

Indexgestützte Funktionen fügen der räumlichen Bedingung automatisch einen Begrenzungsrahmenoperator hinzu. Zu den indexbasierten Funktionen gehören die benannten räumlichen Beziehungsprädikate ST_Contains, ST_ContainsProperly, ST_CoveredBy, ST_Covers, ST_Crosses, ST_Intersects, ST_Overlaps, ST_Touches, ST_Within, ST_Within und ST_3DIntersects sowie die Entfernungsprädikate ST_DWithin, ST_DFullyWithin, ST_3DDFullyWithin und ST_3DDWithin ).

Funktionen wie ST_Distance und nicht verwenden Indizes, um ihren Betrieb zu optimieren. Die folgende Abfrage wäre zum Beispiel bei einer großen Tabelle ziemlich langsam:

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

Diese Abfrage wählt alle Geometrien in geom_table aus, die innerhalb von 100 Einheiten des Punktes (100000, 200000) liegen. Sie ist langsam, weil sie den Abstand zwischen jedem Punkt in der Tabelle und dem angegebenen Punkt berechnet, d. h. eine ST_Distance() Berechnung wird für jede Zeile in der Tabelle berechnet.

Die Anzahl der zu verarbeitenden Zeilen kann durch die Verwendung der indexbasierten Funktion ST_DWithin erheblich reduziert werden:

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

Diese Abfrage wählt dieselben Geometrien aus, allerdings auf effizientere Weise. Dies wird durch ST_DWithin() ermöglicht, die den && Operator intern auf einem erweiterten Begrenzungsrahmen der Abfragegeometrie verwendet. Wenn es einen räumlichen Index auf geom gibt, erkennt der Abfrageplaner, dass er den Index verwenden kann, um die Anzahl der gescannten Zeilen zu reduzieren, bevor die Entfernung berechnet wird. Der räumliche Index ermöglicht es, nur Datensätze mit Geometrien abzurufen, deren Begrenzungsrahmen die erweiterte Ausdehnung überlappen und die daher innerhalb der erforderlichen Entfernung liegen könnten. Der tatsächliche Abstand wird dann berechnet, um zu bestätigen, ob der Datensatz in die Ergebnismenge aufgenommen werden soll.

Weitere Informationen und Beispiele finden Sie im PostGIS Workshop.

5.3. Beispiele für Spatial SQL

Die Beispiele in diesem Abschnitt verwenden eine Tabelle mit linearen Straßen und eine Tabelle mit polygonalen Gemeindegrenzen. Die Definition der Tabelle bc_roads lautet:

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

Die Definition der Tabelle bc_municipality lautet:

Ausgabe von Rastern
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.

Wie lang ist die Gesamtlänge aller Straßen, ausgedrückt in Kilometern?

Sie können diese Frage mit einem sehr einfachen SQL-Programm beantworten:

Code
SELECT sum(ST_Length(geom))/1000 AS km_roads FROM bc_roads;
Ausgabe von Rastern
km_roads
------------------
70842.1243039643

5.3.2.

Wie groß ist die Stadt Prince George (in Hektar)?

Diese Abfrage kombiniert eine Attributbedingung (für den Gemeindenamen) mit einer räumlichen Berechnung (der Polygonfläche):

Code
SELECT
  ST_Area(geom)/10000 AS hectares
FROM bc_municipality
WHERE name = 'PRINCE GEORGE';
Ausgabe von Rastern
hectares
------------------
32657.9103824927

5.3.3.

Welche ist die flächenmäßig größte Gemeinde der Provinz?

Diese Abfrage verwendet ein räumliches Maß als Ordnungswert. Es gibt mehrere Möglichkeiten, dieses Problem anzugehen, aber die effizienteste ist die folgende:

Code
SELECT
  name,
  ST_Area(geom)/10000 AS hectares
FROM bc_municipality
ORDER BY hectares DESC
LIMIT 1;
Ausgabe von Rastern
name           | hectares
---------------+-----------------
TUMBLER RIDGE  | 155020.02556131

Beachten Sie, dass wir zur Beantwortung dieser Abfrage die Fläche jedes Polygons berechnen müssen. Wenn wir dies häufig tun würden, wäre es sinnvoll, der Tabelle eine Flächenspalte hinzuzufügen, die aus Leistungsgründen indiziert werden könnte. Indem wir die Ergebnisse absteigend sortieren und den PostgreSQL-Befehl "LIMIT" verwenden, können wir einfach nur den größten Wert auswählen, ohne eine Aggregatfunktion wie MAX() zu verwenden.

5.3.4.

Wie lang ist die Länge der Straßen, die vollständig in jeder Gemeinde liegen?

Dies ist ein Beispiel für eine "räumliche Verknüpfung", bei der Daten aus zwei Tabellen (mit einer Verknüpfung) unter Verwendung einer räumlichen Interaktion ("enthalten") als Verknüpfungsbedingung zusammengeführt werden (anstelle des üblichen relationalen Ansatzes der Verknüpfung über einen gemeinsamen Schlüssel):

Code
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;
Ausgabe von Rastern
name                        | roads_km
----------------------------+----------
SURREY                      | 1539.476
VANCOUVER                   | 1450.331
LANGLEY DISTRICT            |  833.793
BURNABY                     |  773.769
PRINCE GEORGE               |  694.376
...

Diese Abfrage dauert eine Weile, da jede Straße in der Tabelle im Endergebnis zusammengefasst wird (etwa 250K Straßen für die Beispieltabelle). Bei kleineren Datensätzen (mehrere tausend Datensätze auf mehrere hundert) kann die Antwort sehr schnell sein.

5.3.5.

Erstellen Sie eine neue Tabelle mit allen Straßen innerhalb der Stadt Prince George.

Dies ist ein Beispiel für eine "Überlagerung", die zwei Tabellen aufnimmt und eine neue Tabelle aus räumlich beschnittenen oder geschnittenen Ergebnisgrößen ausgibt. Anders als bei der oben gezeigten "räumlichen Verknüpfung" werden bei dieser Abfrage neue Geometrien erstellt. Ein Overlay ist eine Art "Turbo" für räumliche Verknüpfungen und eignet sich für genauere Analysen:

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

Wie lang ist die "Douglas St" in Victoria in Kilometern?

Code
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';
Ausgabe von Rastern
kilometers
------------------
4.89151904172838

5.3.7.

Welches ist das größte Gemeindepolygon, das ein Loch hat?

Code
SELECT gid, name, ST_Area(geom) AS area
FROM bc_municipality
WHERE ST_NRings(geom) 
> 1
ORDER BY area DESC LIMIT 1;
Ausgabe von Rastern
gid  | name         | area
-----+--------------+------------------
12   | SPALLUMCHEEN | 257374619.430216