PostGIS の LineString と Polygon で路線と公園を描きながら学んでみた
はじめに
GIS 初学者の筆者が PostGIS の LineString（線） と Polygon（面） を試してみました。地下鉄路線を線で、公園の敷地を面で表現し、長さや面積を計算してみます。
前回の記事では Point（点）で 2 点間の距離計算を扱いました。今回はその続きとして線と面を試してみました。知っておいた方が良いこととしてgeometry 型と geography 型の違い についても整理しました。
試したこと早見
今回手を動かして理解できたポイントは以下の通りです。
- LineString で地下鉄路線（南北線・東西線）を線として作成し、長さや頂点数を取得した
- Polygon で公園（大通公園・中島公園）を面として作成し、面積や周長を取得した
- geojson.io でジオメトリを地図上に可視化する方法を試した
- geometry 型と geography 型の違いを整理し、
::geographyキャストでメートル単位に変換できることを確認した - 前回の
ST_DistanceSphereとST_Distance+ geography 型の関係を比較した ST_Collect/ST_Dumpで Multi 系ジオメトリの集約・分解を試した
環境の準備
前回の記事で構築した Docker 環境をそのまま使用します。
# 前回のディレクトリに移動
cd learn-gis
# コンテナが起動していなければ起動
docker compose up -d
# psqlで接続
docker compose exec db psql -U gisuser -d gis_db
前回作成した landmarks テーブルが残っている状態で進めます。
LineString で線を扱う
LineString とは
まず LineString（ラインストリング） について調べました。2 つ以上の点を順に結んだ線を表すジオメトリ型で、道路、鉄道路線、河川など「線的な地物」の表現に使います。
WKT（Well-Known Text）での表記は以下の通りです。
LINESTRING(経度1 緯度1, 経度2 緯度2, 経度3 緯度3, ...)
ラインテーブルの作成
テーブルを作ります。GEOMETRY(LineString, 4326) で LineString 型かつ SRID 4326（WGS84）のみ受け付けるカラムを定義します。
CREATE TABLE subway_lines (
id SERIAL PRIMARY KEY,
name TEXT NOT NULL,
line_name TEXT NOT NULL,
geom GEOMETRY(LineString, 4326)
);
札幌地下鉄南北線を LineString で作成する
前回使った札幌の座標を流用して、地下鉄南北線の中心部区間を LineString で表現してみました。さっぽろ駅から中島公園駅までの 4 駅を線で結びます。
INSERT INTO subway_lines (name, line_name, geom) VALUES
('南北線（さっぽろ〜中島公園）', '南北線',
ST_SetSRID(ST_GeomFromText(
'LINESTRING(141.3507 43.0687, 141.3544 43.0607, 141.3533 43.0555, 141.3549 43.0487)'
), 4326));
-- さっぽろ駅 大通駅 すすきの駅 中島公園駅
もう 1 本追加して東西線を作成する
INSERT INTO subway_lines (name, line_name, geom) VALUES
('東西線（西11丁目〜バスセンター前）', '東西線',
ST_SetSRID(ST_GeomFromText(
'LINESTRING(141.3383 43.0588, 141.3509 43.0596, 141.3563 43.0603, 141.3634 43.0610)'
), 4326));
-- 西11丁目駅 西4丁目付近 大通駅 バスセンター前駅
WKT で確認してみる
SELECT name, ST_AsText(geom) AS wkt FROM subway_lines;
name | wkt
------------------------------------+---------------------------------------------------------------------------------
南北線（さっぽろ〜中島公園） | LINESTRING(141.3507 43.0687,141.3544 43.0607,141.3533 43.0555,141.3549 43.0487)
東西線（西11丁目〜バスセンター前） | LINESTRING(141.3383 43.0588,141.3509 43.0596,141.3563 43.0603,141.3634 43.061)
(2 rows)
LineString の基本操作を試す
| 関数 | 意味 |
|---|---|
ST_NPoints |
LineString の頂点数 |
ST_StartPoint |
始点（最初の頂点） |
ST_EndPoint |
終点（最後の頂点） |
ST_Length |
線の長さ |
SELECT
name,
ST_NPoints(geom) AS num_points,
ST_AsText(ST_StartPoint(geom)) AS start_point,
ST_AsText(ST_EndPoint(geom)) AS end_point,
ST_Length(geom) AS length_degree
FROM subway_lines;
ここで length_degree が 0.021... のような小さな値になっていて意味がわかりませんでした。これは「度」単位の数値であり、「geometry 型と geography 型」セクションで整理します。
name | num_points | start_point | end_point | length_degree
------------------------------------+------------+-------------------------+-------------------------+----------------------
南北線（さっぽろ〜中島公園） | 4 | POINT(141.3507 43.0687) | POINT(141.3549 43.0487) | 0.02111496568477214
東西線（西11丁目〜バスセンター前） | 4 | POINT(141.3383 43.0588) | POINT(141.3634 43.061) | 0.025204976230397646
(2 rows)
Polygon で面を扱う
Polygon とは
次に Polygon（ポリゴン） を試しました。閉じた領域を表すジオメトリ型で、行政区域、建物の敷地、公園の範囲など「面的な地物」の表現に使います。
WKT での表記は以下の通りです。
POLYGON((経度1 緯度1, 経度2 緯度2, 経度3 緯度3, ..., 経度1 緯度1))
ポリゴンテーブルの作成
CREATE TABLE parks (
id SERIAL PRIMARY KEY,
name TEXT NOT NULL,
geom GEOMETRY(Polygon, 4326)
);
大通公園を Polygon で作成してみる
大通公園は東西に約 1.5km 延びる細長い公園です。概略的な矩形で表現してみます。
INSERT INTO parks (name, geom) VALUES
('大通公園',
ST_SetSRID(ST_GeomFromText(
'POLYGON((141.3410 43.0600, 141.3565 43.0608, 141.3565 43.0598, 141.3410 43.0590, 141.3410 43.0600))'
), 4326));
-- ↑ NW ↑ NE ↑ SE ↑ SW ↑ 始点=終点（閉じる）
もう 1 つ追加して中島公園を作成する
中島公園は好きなのでキャプチャも貼っておきます。中島公園の座標は Google Maps で公園の四隅をクリックして取得しました。
実際に地図上で 4 つの角をクリックしてみると、「始点 = 終点で閉じる」という Polygon のルールが腑に落ちました。地図では 4 箇所をクリックするだけですが、WKT では最初の座標を末尾にもう一度書いて閉じます。そのため座標は 5 つになります。
INSERT INTO parks (name, geom) VALUES
('中島公園',
ST_SetSRID(ST_GeomFromText(
'POLYGON((141.3567 43.0476, 141.3565 43.0402, 141.3518 43.0401, 141.3512 43.0480, 141.3567 43.0476))'
), 4326));
穴あきポリゴンでドーナツ型を表現する
公園の中に池や建物がある場合、穴あき Polygon で表現できます。複数のリングを定義し、1 番目を外周、2 番目以降を穴とします。中島公園の池をくり抜こうかと思ったのですが、最初はシンプルに試したかったので座学で済ませました。
| リング | 回転方向 | 役割 |
|---|---|---|
| 1 番目のリング | 反時計回り | 外周（ポリゴンの境界） |
| 2 番目以降のリング | 時計回り | 穴（くり抜く領域） |
-- ドーナツ型Polygonの例（投入はしない。構文の解説のみ）
SELECT ST_AsText(
ST_GeomFromText(
'POLYGON(
(141.3567 43.0476, 141.3565 43.0402, 141.3518 43.0401, 141.3512 43.0480, 141.3567 43.0476),
(141.3530 43.0460, 141.3528 43.0430, 141.3520 43.0430, 141.3522 43.0460, 141.3530 43.0460)
)', 4326)
) AS donut_polygon;
-- ↑ 外周リング（反時計回り） ↑ 穴リング（時計回り）
Polygon の基本操作を試す
| 関数 | 意味 |
|---|---|
ST_NPoints |
頂点数（始点=終点を含むため、4 頂点の矩形は 5 になる） |
ST_Centroid |
ポリゴンの重心（中心点） |
ST_Area |
面積 |
ST_Perimeter |
周長（外周の長さ） |
SELECT
name,
ST_NPoints(geom) AS num_points,
ST_AsText(ST_Centroid(geom)) AS centroid,
ST_Area(geom) AS area_degree,
ST_Perimeter(geom) AS perimeter_degree
FROM parks;
ここでも面積が 0.0000... という小さな値、周長も度単位になっています。
name | num_points | centroid | area_degree | perimeter_degree
----------+------------+---------------------------------------------+------------------------+----------------------
大通公園 | 5 | POINT(141.34874999999997 43.0599) | 1.5500000000021898e-05 | 0.03304126286091674
中島公園 | 5 | POINT(141.3540209354025 43.044074220854164) | 3.8984999999998055e-05 | 0.025541044237441385
(2 rows)
geojson.io でジオメトリを地図上に表示する
Polygon や LineString は座標の羅列だけでは形状がまったく想像できませんでした。調べたところ、geojson.io を使えばブラウザ上で地図表示できることがわかりました。
手順
- PostGIS で GeoJSON を出力する
-- 中島公園のPolygonをGeoJSONとして出力
SELECT json_build_object(
'type', 'FeatureCollection',
'features', json_agg(
json_build_object(
'type', 'Feature',
'properties', json_build_object('name', name),
'geometry', ST_AsGeoJSON(geom)::json
)
)
) AS geojson
FROM parks
WHERE name = '中島公園';
- 出力された JSON をコピーする
- https://geojson.io にアクセスし、右側の JSON エディターに貼り付ける
- 左側の地図に Polygon が表示される
灰色の網掛けで中島公園が囲われていることを確認できました。
geojson.io はインストール不要で、GeoJSON を貼り付けるだけで地図表示できます。座標の入力ミス（lon/lat の逆転、Polygon が閉じていない等）にもすぐ気づけます。PostGIS のデバッグツールとして覚えておくと便利そうです。
geometry 型と geography 型について
前回 ST_DistanceSphere を使ったとき、結果は「メートル」で返ってきました。しかし今回 ST_Length や ST_Area を使うと「度」で返ってきます。調べてみると geometry 型と geography 型の違いが原因でした。
| geometry 型 | geography 型 | |
|---|---|---|
| 計算方式 | 平面（デカルト座標） として計算 | 球面（地球の曲面） として計算 |
| EPSG:4326 での単位 | 度（経度・緯度の数値差） | メートル |
| 関数の対応 | ほぼすべての ST_* 関数が使える | 一部の関数のみ対応 |
| パフォーマンス | 高速 | geometry 型より遅い |
| 用途 | 局所的な解析、投影座標系のデータ | 広域の距離・面積計算 |
ST_Length の違いを実際に確認してみる
同じ LineString に対して、geometry 型と geography 型で ST_Length を比較してみました。
SELECT
name,
-- geometry型: 平面計算 → 度が返る
ST_Length(geom) AS length_degree,
-- geography型: 球面計算 → メートルが返る
ROUND(ST_Length(geom::geography)::numeric, 0) AS length_meter
FROM subway_lines;
0.021度: 数値としては正しいが、人間には意味がわからない2290メートル: 約 2.3km でさっぽろ駅から中島公園駅までの距離はこのくらいです
name | length_degree | length_meter
------------------------------------+----------------------+--------------
南北線（さっぽろ〜中島公園） | 0.02111496568477214 | 2290
東西線（西11丁目〜バスセンター前） | 0.025204976230397646 | 2061
ST_Area の違いも確認してみる
Polygon の面積でも同じなのか試してみました。
SELECT
name,
-- geometry型: 度² が返る（人間にはわからない値）
ST_Area(geom) AS area_degree2,
-- geography型: 平方メートルが返る
ROUND(ST_Area(geom::geography)::numeric, 0) AS area_m2,
-- 平方メートル → ヘクタール（1ha = 10,000m²）
ROUND((ST_Area(geom::geography) / 10000)::numeric, 2) AS area_ha
FROM parks;
同じですね。
name | area_degree2 | area_m2 | area_ha
----------+------------------------+---------+---------
大通公園 | 1.5500000000021898e-05 | 140273 | 14.03
中島公園 | 3.8984999999998055e-05 | 352899 | 35.29
ST_DistanceSphere との関係
ここで前回使った ST_DistanceSphere との関係が気になったので調べました。ST_DistanceSphere は geometry 型を受け取りつつ内部で球面計算する特殊な関数です。2 つの方法を比較してみます。
| 方法 | 書き方 | 内部の計算 |
|---|---|---|
ST_DistanceSphere(geom, geom) |
geometry 型のまま渡す | 球体近似で距離計算 |
ST_Distance(geog, geog) |
geography 型にキャスト | 回転楕円体で距離計算（より精密） |
-- 前回の方法: ST_DistanceSphere（geometry型を渡すが球面計算する）
SELECT ROUND(ST_DistanceSphere(
ST_SetSRID(ST_MakePoint(141.3507, 43.0687), 4326),
ST_SetSRID(ST_MakePoint(141.3549, 43.0487), 4326)
)::numeric, 0) AS method1_m;
-- 今回の方法: geography型にキャストしてST_Distance
SELECT ROUND(ST_Distance(
ST_SetSRID(ST_MakePoint(141.3507, 43.0687), 4326)::geography,
ST_SetSRID(ST_MakePoint(141.3549, 43.0487), 4326)::geography
)::numeric, 0) AS method2_m;
どちらもメートルを返しますが、ST_Distance + geography 型のほうが高精度です。ST_DistanceSphere は手軽に球面距離を出したいときの便利関数です。
なぜ geometry 型が基本なのか整理してみた
「度」で返るのが不便なら geography 型だけ使えばよいのでは？と思いました。調べてみると、geometry 型が基本になっている理由がありました。実際はほとんどの空間操作では「度」が問題にならないようです。
| 操作の種類 | 例 | 単位が影響するか |
|---|---|---|
| 含まれるか？ | ST_Contains(区域, 地点) |
しない |
| 交差するか？ | ST_Intersects(路線, 区域) |
しない |
| 最寄りを探す | ORDER BY geom <-> target LIMIT 1 |
しない |
| バウンディングボックス検索 | geom && ST_MakeEnvelope(...) |
しない |
| 距離を測る | ST_Length, ST_Distance |
する |
| 面積を測る | ST_Area |
する |
「含まれるか」「交差するか」といったトポロジ判定は、日本国内の局所的なデータであれば geometry 型でも実用上問題ないようです。GIS の実務では距離・面積の計測よりトポロジ判定のほうが多いようです。
Multi 系ジオメトリで複数パーツを表現する
現実の地物は 1 つのジオメトリで表現できないケースがあることを知りました。
| 例 | 問題 | 解決 |
|---|---|---|
| 沖縄県（本島 + 離島群） | 1 つの Polygon では表現不可 | MultiPolygon |
| 分岐する河川 | 1 本の LineString では表現不可 | MultiLineString |
| 同じ名前の複数拠点 | 1 つの Point では表現不可 | MultiPoint |
WKT での表記
MultiPoint:
MULTIPOINT((141.35 43.06), (141.36 43.07))
MultiLineString:
MULTILINESTRING((141.35 43.06, 141.36 43.07), (141.37 43.08, 141.38 43.09))
MultiPolygon:
MULTIPOLYGON(
((141.34 43.05, 141.35 43.05, 141.35 43.06, 141.34 43.06, 141.34 43.05)),
((141.36 43.05, 141.37 43.05, 141.37 43.06, 141.36 43.06, 141.36 43.05))
)
ST_Collect で複数ジオメトリを 1 つにまとめる
ST_Collect は集約関数（SUM, COUNT 等と同様）で、複数のジオメトリを 1 つにまとめます。同じ型同士なら Multi 系に、異なる型が混在する場合は GeometryCollection になります。GROUP BY と組み合わせて「区ごとにポリゴンを集約」といった用途にも使えます。
-- 2つの公園を1つのMultiPolygonにまとめる
SELECT
ST_GeometryType(ST_Collect(geom)) AS collected_type,
ST_NumGeometries(ST_Collect(geom)) AS num_parts
FROM parks;
collected_type | num_parts
----------------------+-----------
ST_MultiPolygon | 2
ST_Dump で Multi 系を個々に分解する
ST_Dump は ST_Collect の逆で、Multi 系ジオメトリを個々のジオメトリに分解します。
| 関数 | 方向 | 用途 |
|---|---|---|
ST_Collect |
個別 → Multi | 複数ジオメトリを 1 つにまとめる |
ST_Dump |
Multi → 個別 | Multi 系を個々のジオメトリに分解する |
ST_NumGeometries |
— | Multi 系に含まれるジオメトリ数を返す |
ST_GeometryN(geom, n) |
— | n 番目のジオメトリを取り出す（1-indexed） |
SELECT
(dump).path[1] AS part_number,
ST_AsText((dump).geom) AS individual_geom
FROM (
SELECT ST_Dump(ST_Collect(geom)) AS dump FROM parks
) AS sub;
part_number | individual_geom
-------------+----------------------------------------------------------------------------------------
1 | POLYGON((141.341 43.06,141.3565 43.0608,141.3565 43.0598,141.341 43.059,141.341 43.06))
2 | POLYGON((141.3505 43.05,141.358 43.05,141.358 43.0455,141.3505 43.0455,141.3505 43.05))
まとめ
地下鉄南北線・東西線を題材に、LINESTRING(lon lat, ...) で線を表現する方法を学びました。ST_NPoints、ST_StartPoint、ST_EndPoint、ST_Length で頂点数や始点・終点・長さを取得できます。
大通公園・中島公園を題材に、POLYGON((lon lat, ...)) で面を表現する方法を学びました。始点と終点を同一座標で閉じるルールがあります。ST_Centroid、ST_Area、ST_Perimeter で重心・面積・周長を取得できます。
座標の羅列だけでは形状がわからないため、ST_AsGeoJSON で GeoJSON を出力し geojson.io に貼り付けて地図上で確認する方法も試しました。
EPSG:4326 の geometry 型で ST_Length / ST_Area を使うと「度」が返ります。::geography キャストでメートル単位に変換できます。前回使った ST_DistanceSphere は geometry 型を渡しつつ内部で球面計算する便利関数で、ST_Distance + geography 型のほうが高精度です。
Multi 系ジオメトリでは、ST_Collect で複数ジオメトリを 1 つに集約し、ST_Dump で個々に分解できます。2 つの公園を MultiPolygon にまとめる例で試しました。
おわりに
今回は手入力の座標で LineString と Polygon の基本を学びました。次回は国土数値情報などの実データを PostGIS にインポートし、ST_Contains や ST_Intersects を使った空間検索をやってみようと思います。
