空間データベースのレゾンデートルは、通常ならデスクトップGISの機能が必要なクエリをデータベース内で実行することです。PostGISを使うには、使用可能な空間関数は何かを知り、またクエリ内でどう使うかを知って、適切なインデックスで能率を向上させることが求められます。
空間関係は、二つのジオメトリについて、一方がもう一方にどのような相互関係になっているかを示すものです。ジオメトリのクエリにおける基本的な機能です。
OpenGIS Simple Features Implementation Specification for SQLによると「二つのジオメトリの比較の基本的なアプローチは、二つのジオメトリの内部、境界、外部のインタセクションの比較と、『インタセクション行列』の要素に基づく2ジオメトリの関係の分類です」。
点集合トポロジでは、2次元空間に埋め込まれたジオメトリの中にあるポイントは、次に示す三つの集合に分類されます。
ジオメトリの境界は、一次元低いジオメトリです。POINTでは、次元が0になり、境界は空集合です。LINESTRINGの境界は二つの端点です。POLYGONの境界は、外環と内環の線です。
ジオメトリの内部は、ジオメトリの境界以外のポイントです。POINTでは、内部はポイント自体です。LINESTRINGの内部は端点の間のポイントの集合です。POLYGONの内部は、ポリゴン内部の面です。
ジオメトリの外部はジオメトリが組み込まれた空間の残りです。言い換えると、ジオメトリの内部にも境界にもない点の全てです。これは2次元の閉じていない面になります。
次元拡張9交差モデル (Dimensionally Extended 9-Intersection Model, DE-9IM)は、二つのジオメトリの空間関係を九つの交差の次元を指定することで記述します。交差次元は3×3の交差行列で正式に表現することができます。
ジオメトリgに対する内部、境界、外部はI(g)、B(g)、E(g)と表記します。また、dim(s)はsの集合を{0,1,2,F}の値で示すます。
0 => 点
1 => 線
2 => 面
F => 空集合
この表記法を使うと、二つのジオメトリaとbの交差行列は次の通りです。
| 内部 (Interior) | 境界 (Boundary) | 外部 (Exterior) | |
|---|---|---|---|
| 内部 (Interior) | dim( I(a) ∩ I(b) ) | dim( I(a) ∩ B(b) ) | dim( I(a) ∩ E(b) ) |
| 境界 (Boundary) | dim( B(a) ∩ I(b) ) | dim( B(a) ∩ B(b) ) | dim( B(a) ∩ E(b) ) |
| 外部 (Exterior) | 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.
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;
matrix | matches_overlap_pattern | overlaps -----------+-------------------------+---------- 212101212 | t | t (1 row)
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.
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;
左から右に、上から下に読みます。交差行列の文字列表現は'212101212'です。
詳細情報については次をご覧下さい。
共通の空間関係を簡単に決定できるように、PGC SFSは名前付き空間関係述語の集合を定義しています。PostGISではST_Contains、ST_Crosses、ST_Disjoint、ST_Equals、ST_Intersects、ST_Overlaps、ST_Touches、ST_Withinが提供されています。非標準の空間関係述語ST_Covers、ST_CoveredBy、ST_ContainsProperlyも定義されています。
空間述語は通常SQLのWHERE節やJOIN節内で条件に使用されます。名前付き空間述語は、インデックスが有効なら自動的に空間インデックスを使うので、バウンディングボックス演算子&&を使う必要はありません。例えば次のようになります。
SELECT city.name, state.name, city.geom FROM city JOIN state ON ST_Intersects(city.geom, state.geom);
詳細や図についてはPostGIS Workshopをご覧下さい。
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.
特定の空間関係をテストするには、交差行列パターンを使います。これは、追加シンボル{T,*}で拡張された行列表現です。
T => インタセクションの次元は空ではないという意味です。すなわち{0,1,2}のいずれかです。
* => 何でも良い
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.
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;
overlap | matrix | ST_Crosses | matches 1*1***1** ----------------------------+-----------+------------+------------------- LINESTRING(90 110,130 110) | 1F1FF0102 | f | t (1 row)
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.
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;
wharf | matching wharf | matrix
-------+-----------------------------------+-----------
1 | LINESTRING(120 180,145 85,110 80) | 102101FF2
(1 row)
空間条件を使用するクエリを構築する時、最良の効果を得るには、空間インデックスが存在する場合に (「空間インデックス」参照)これを確実に使用することが重要です。そのためには、WHERE節やON節で、空間演算子またはインデックス対応関数を使用しなければなりません。
空間演算子には、バウンディングボックス演算子 (最もよく使われるのは&&です。「バウンディングボックス演算子」参照)、および近傍クエリで使用される距離演算子 (最もよく使われるのは<->です。「距離演算子」参照)が含まれます。
インデックス対応関数は、自動的にバウンディングボックス演算子を空間条件に追加します。インデックス対応関数は空間関係述語を含みます。空間関係述語には、ST_Contains, ST_ContainsProperly, ST_CoveredBy, ST_Covers, ST_Crosses, ST_Intersects, ST_Overlaps, ST_Touches, ST_Within, ST_Within, ST_3DIntersectsがあり、距離述語にはST_DWithin, ST_DFullyWithin, ST_3DDFullyWithin, ST_3DDWithin があります。
ST_Distanceといった関数は、演算の最適化のためにはインデックスを使用しません。例えば、次のクエリは、大きなテーブルでは非常に遅くなります。
SELECT geom FROM geom_table WHERE ST_Distance(geom, 'SRID=312;POINT(100000 200000)') < 100
このクエリはgeom_tableテーブル内の、(100000, 200000)のポイントから100単位内にある全てのジオメトリを選択します。テーブル内の個々のポイントと指定したポイントとの距離を計算しているため、非常に遅くなります。すなわち、1回のST_Distance()の計算で、テーブルの全ての行について計算することになります。
インデックス対応関数ST_DWithinを使用すると、処理行数を実質的に減らすことができます。次のようにします。
SELECT geom FROM geom_table WHERE ST_DWithin(geom, 'SRID=312;POINT(100000 200000)', 100)
このクエリは、同じジオメトリを選択しますが、より効率的な方法を取ります。 ST_DWithin()が内部で&&演算子をクエリジオメトリのバウンディングボックスを拡大したボックスで使うことによって可能となります。geom上に空間インデックスが存在するなら、クエリプランナは距離計算の前に対象行数を減らすためにインデックスを使えることを認識します。空間インデックスによって、バウンディングボックスが拡張された範囲とオーバラップするジオメトリだけを検索して、そのため、求めようとする距離内にあるかも知れないジオメトリを検索することができます。その後で、結果集合内のレコードを含めるかどうかを確認するための実際の距離計算が行われます。
詳細情報と例についてはPostGIS Workshopをご覧下さい。
本節の例では、線の道路のテーブルとポリゴンの市区町村境界テーブルとを使います。bc_roadsテーブルの定義は次の通りです。
Column | Type | Description ----------+-------------------+------------------- gid | integer | Unique ID name | character varying | Road Name geom | geometry | Location Geometry (Linestring)
bc_municipalityテーブルの定義は次の通りです。
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. |
道路の総延長はkm表記でいくらになるでしょう? |
|
この問題は、次のようなとても単純なSQLで答えを得ることができます。 Code
SELECT sum(ST_Length(geom))/1000 AS km_roads FROM bc_roads; 出力:
km_roads ------------------ 70842.1243039643 |
|
|
5.3.2. |
プリンスジョージ市の大きさはha表記でいくらになるでしょう? |
|
このクエリでは、属性条件 (municipality name, 自治体名)に (ポリゴン面積の)空間計算を併用しています。 Code
SELECT ST_Area(geom)/10000 AS hectares FROM bc_municipality WHERE name = 'PRINCE GEORGE'; 出力:
hectares ------------------ 32657.9103824927 |
|
|
5.3.3. |
県内で最も大きな面積となる自治体はどこでしょう? |
|
このクエリでは、順序付けの値に空間計測関数を使っています。この問題に対しては複数の方法がありますが、最も効果的な方法は次の通りです。 Code
SELECT name, ST_Area(geom)/10000 AS hectares FROM bc_municipality ORDER BY hectares DESC LIMIT 1; 出力:
name | hectares ---------------+----------------- TUMBLER RIDGE | 155020.02556131 このクエリの答えを出すためには、全てのポリゴンの面積を求める必要があることに注意して下さい。このクエリを多く実行する場合、性能向上のためにテーブルにareaカラムを追加して、別のインデックスを追加することができるようにするのは、意義のあることです。結果を距離について降順に並べ替え、PostgreSQLの"LIMIT"コマンドを用いることで、max()のような集約関数を使わずに、簡単に最も大きい値を集約関数を得ることができます。 |
|
|
5.3.4. |
各自治体内に含まれる道路の総延長はいくらでしょう? |
|
これは、二つのテーブルからデータを持ち込んで (結合して)いるので「空間結合」の例です。しかし、結合の条件として共通キーの上で接続するという普通のリレーションのやり方でなく空間インタラクション条件 (「含む」)を使っています。 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; 出力:
name | roads_km ----------------------------+---------- SURREY | 1539.476 VANCOUVER | 1450.331 LANGLEY DISTRICT | 833.793 BURNABY | 773.769 PRINCE GEORGE | 694.376 ... このクエリは、テーブル内の全ての道路の合計を最終結果 (この例での話ですが約250Kmの道です)にまとめられるので、少し時間がかかります。より小さいオーバレイ (数百の道路で数千のレコード)の場合、応答はもっと早くなりえます。 |
|
|
5.3.5. |
プリンスジョージ市内の全ての道路からなるテーブルを作ります。 |
|
これは「オーバレイ」の例です。つまり、二つのテーブルを取得して、空間的に切り取られた結果からなる新しいテーブルを出力します。上で示した「空間結合」と違い、このクエリは実際に新しいジオメトリを生成します。生成されたオーバレイはターボのかかった空間結合みたいなもので、より確かな解析作業に便利です。 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. |
ビクトリア州の「ダグラス通り」の長さはkm表記でいくらになるでしょう? |
|
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'; 出力:
kilometers ------------------ 4.89151904172838 |
|
|
5.3.7. |
穴を持つ自治体ポリゴンのうち最も大きいのはどれでしょう? |
|
Code
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 |