空間データベースのレゾンデートルは、通常ならデスクトップ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) ) |
下に重なっているポリゴンを具体例としてしまします。ST_Relate は行列 212101212 を算出し、T*T***T** パターンと合致します。このパターンは 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)
実行可能クエリーと生成されるオーバーレイの上の部分が関係を見る元図です。クエリーの下の部分で、3x3 で 9要素の同種の例を生成しています。最後の外部/外部の要素は無限に広がるので有限範囲のビューポートで切り抜いています。
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をご覧下さい。
名前のある空間関係は、求める空間フィルター条件を十分に満たせない場合があります。この場合には、ST_Relate を使って完全な DE-9IM 交差行列による計算によって表現することができます。
特定の空間関係をテストするには、交差行列パターンを使います。これは、追加シンボル{T,*}で拡張された行列表現です。
T => インタセクションの次元は空ではないという意味です。すなわち{0,1,2}のいずれかです。
* => 何でも良い
交差行列のパターンを使って、特定の空間関係が簡潔に評価できます。ST_Relate と ST_RelateMatch は両方ともこのパターンと合致するかの判定ができます。
次に示す二つの道路区間は、ラインの一部を共有しています。ST_Crosses にラインを与えた場合には、内部の点でインターセクトする時にだけ TRUE を返すので、この関係を判定できません。1*1***1** パターンはラインのインターセクションを判定します。
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)
次の例では、湖の内部から水涯線上に伸びるふ頭を検索します。端点の一方が水涯線上にあることを条件とします。他のふ頭は該当しない例です。102101FF2 は、求める関係の完全なものを表現しています。
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 で答えを得ることができます。 コード
SELECT sum(ST_Length(geom))/1000 AS km_roads FROM bc_roads; 出力
km_roads ------------------ 70842.1243039643 |
|
|
5.3.2. |
プリンスジョージ市の大きさは ha 表記でいくらになるでしょう? |
|
このクエリでは、属性条件 (municipality name, 自治体名) に (ポリゴン面積の) 空間計算を併用しています。 コード
SELECT ST_Area(geom)/10000 AS hectares FROM bc_municipality WHERE name = 'PRINCE GEORGE'; 出力
hectares ------------------ 32657.9103824927 |
|
|
5.3.3. |
県内で最も大きな面積となる自治体はどこでしょう? |
|
このクエリでは、順序付けの値に空間計測関数を使っています。この問題に対しては複数の方法がありますが、最も効果的な方法は次の通りです。 コード
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. |
各自治体内に含まれる道路の総延長はいくらでしょう? |
|
これは、二つのテーブルからデータを持ち込んで (結合して) いるので「空間結合」の例です。しかし、結合の条件として共通キーの上で接続するという普通のリレーションのやり方でなく空間インタラクション条件 (「含む」) を使っています。 コード
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. |
プリンスジョージ市内の全ての道路からなるテーブルを作ります。 |
|
これは「オーバレイ」の例です。つまり、二つのテーブルを取得して、空間的に切り取られた結果からなる新しいテーブルを出力します。上で示した「空間結合」と違い、このクエリは実際に新しいジオメトリーを生成します。生成されたオーバレイはターボのかかった空間結合みたいなもので、より確かな解析作業に便利です。 コード
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 表記でいくらになるでしょう? |
|
コード
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. |
穴を持つ自治体ポリゴンのうち最も大きいのはどれでしょう? |
|
コード
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 |