第4章 PostGISを使う: データ管理とクエリ

OpenGIS仕様は空間オブジェクトの表現について二つの標準を定義しています。Well-Knownテキスト (WKT)形式とWell-Knownバイナリ (WKB)形式です。WKTもWKBも、オブジェクトの型とオブジェクトを形成する座標に関する情報を持っています。

フィーチャーの空間オブジェクトのテキスト表現 (WKT)の例は、次の通りです。

OpenGIS仕様では、空間オブジェクトの内部保存書式は空間参照系識別子 (Spatial Referencing System IDentifier, SRID)を含むことも求められます。SRIDはデータベースへの挿入のために空間オブジェクトが生成される時に求められます。

これらの書式の入出力は次のインタフェースを用いて実現できます。

バイト配列 WKB = ST_AsBinary(geometry);
テキスト WKT = ST_AsText(geometry);
ジオメトリ = ST_GeomFromWKB (bytea WKB、SRID);
ジオメトリ = ST_GeometryFromText (テキスト WKT、SRID);

たとえば、OGC空間オブジェクトを生成して挿入する妥当なINSERTステートメントは次の通りです。

INSERT INTO geotable ( the_geom, the_name )
  VALUES ( ST_GeomFromText('POINT(-126.4 45.32)', 312), 'A Place');

OGC書式は2次元ジオメトリしかサポートされておらず、また、入出力の表現においてSRIDは*決して*埋め込まれません。

PostGIS拡張書式は現在のところOGC書式のスーパーセットとなっています (全ての妥当なWKB/WKTは妥当なEWKB/EWKTです)。しかし、特にもしOGCがPostGIS拡張と矛盾する新しい書式を出すことがあるなら、これは将来変更されるかも知れません。ゆえにこの機能に頼るべきではありません。

PostGIS EWKB/EWKT では 3DM, 3DZ, 4D の座標サポートが追加され、SRID情報が埋め込まれます。

拡張された空間オブジェクトのテキスト表現 (EWKT)の例は、次の通りです。

これらの書式の変換は次のインタフェースを用いて実現できます。

バイト配列 EWKB = ST_AsEWKB(geometry);
テキスト EWKT = ST_AsEWKT(geometry);
ジオメトリ = ST_GeomFromEWKB(bytea EWKB);
ジオメトリ = ST_GeomFromEWKT(text EWKT);

たとえば、PostGISの空間オブジェクトを作成し挿入する妥当なINSERTステートメントは次の通りです。

INSERT INTO geotable ( the_geom, the_name )
  VALUES ( ST_GeomFromEWKT('SRID=312;POINTM(-126.4 45.32 15)'), 'A Place' )

PostgreSQLの「標準的な形式」は単純なクエリ (全く関数呼び出しが無い)で得る表現であり、単純なINSERT, UPDATE, COPYで受け付けられることが保障されるものです。PostGISの"geometory"型の場合は次の通りです。

- 出力
  - バイナリ: EWKB
        ascii: HEXEWKB (EWKBのHEX表現)
- 入力
  - バイナリ: EWKB
        ascii: HEXEWKB|EWKT 

たとえば、次のステートメントは、標準的なASCII文字列による入出力の処理でEWKTを読み、HEXEWKBを返すものです。

=# SELECT 'SRID=4;POINT(0 0)'::geometry;

geometry
----------------------------------------------------
01010000200400000000000000000000000000000000000000
(1 row)

SQLマルチメディア・アプリケーション空間仕様は、円弧補完曲線を定義したSQL仕様の拡張です。

SQL-MMの定義では、3DM、3DZと4Dの座標を含みますが、SRID情報の埋め込みはできません。

Well-Known Text拡張はまだ完全にはサポートされていません。単純な曲線ジオメトリの例を次に示します。

	SQL-MM実装での全ての浮動小数点数の比較では、所定の丸め誤差があります。現在は1E-8です。

ジオグラフィ型はシンプルフィーチャーの最も簡単なもののみサポートします。標準的なジオメトリ型データで、SRIDが4326の場合は、ジオグラフィに自動でキャストされます。またデータ挿入においてEWKTとEWKBの取り決めを使うこともできます。

ジオグラフィ型のフィールドはgeography_columnsシステムビューに登録されます。

"geography_columns"ビューをチェックして、テーブルが一覧にあるか見て下さい。

CREATE TABLEの文法でジオグラフィカラムを持つテーブルを新規に生成できます。

CREATE TABLE global_points (
    id SERIAL PRIMARY KEY,
    name VARCHAR(64),
    location GEOGRAPHY(POINT,4326)
  );

locationカラムはジオグラフィ型で、二つの任意修飾子に対応していることにご注意ください。一つは、そのカラムで使用できる形状と次元を限定する型修飾子です。もう一つは、座標参照IDを特定の数に限定するSRID修飾子です。

型修飾子として受け付ける値は、POINT, LINESTRING, POLYGON, MULTIPOINT, MULTILINESTRING, MULTIPOLYGONです。型修飾子は、Z, M, ZMの後置子を付けることで次元制限にも対応します。たとえば、'LINESTRINGM'という型修飾子で、三つ目の次元をMとする3次元のラインストリングのみ受け付けることになります。同じように'POINTZM'で、4次元データを期待します。

SRIDを指定していない場合には、SRIDは4326 WGS84経度緯度が使われ、全ての計算はWGS 84を使って行われます。

テーブルを作ったら、次のようにしてGEOGRAPHY_COLUMNSを見ることができます。

-- メタデータビューの中身を見る
SELECT * FROM geography_columns;

ジオメトリカラムを使うのと同じようにテーブルへのデータの挿入ができます。

-- testテーブルにデータを追加する
INSERT INTO global_points (name, location) VALUES ('Town', 'SRID=4326;POINT(-110 30)');
INSERT INTO global_points (name, location) VALUES ('Forest', 'SRID=4326;POINT(-109 29)');
INSERT INTO global_points (name, location) VALUES ('London', 'SRID=4326;POINT(0 49)');

ジオメトリと同じ操作でインデックスを作成します。PostGISは、カラム型がジオグラフィであるかを見て、ジオメトリで使われる平面用インデックスの代わりに球面ベースのインデックスを作成します。

-- testテーブルに球面インデックスを作成
  CREATE INDEX global_points_gix ON global_points USING GIST ( location );

クエリと計測関数はメートル単位となります。そのため距離パラメータはメートル (面積の場合は平方メートル)単位となります。

-- 距離クエリの表示。ロンドンは1000km範囲外です
  SELECT name FROM global_points WHERE ST_DWithin(location, 'SRID=4326;POINT(-110 29)'::geography, 1000000);

ジオグラフィの威力については、シアトルからロンドンまで (LINESTRING(-122.33 47.606, 0.0 51.5))の飛行機がレイキャビク (POINT(-21.96 64.15))に最も近くなるときの距離を求めてみると分かります。

-- ジオグラフィを使った距離計算 (122.2km)
  SELECT ST_Distance('LINESTRING(-122.33 47.606, 0.0 51.5)'::geography, 'POINT(-21.96 64.15)'::geography);

-- ジオメトリを使った計算 (13.3 "度")
  SELECT ST_Distance('LINESTRING(-122.33 47.606, 0.0 51.5)'::geometry, 'POINT(-21.96 64.15)'::geometry);

異なる経度緯度座標系を試します。spatial_ref_sysテーブルで挙げられている経度緯度の空間参照系なら全て可能です。

-- NAD83経度緯度
SELECT 'SRID=4269;POINT(-123 34)'::geography;
                    geography
----------------------------------------------------
 0101000020AD1000000000000000C05EC00000000000004140
(1 row)

-- NAD27経度緯度
SELECT 'SRID=4267;POINT(-123 34)'::geography;

                    geography
----------------------------------------------------
 0101000020AB1000000000000000C05EC00000000000004140
(1 row)

-- メートル単位のNAD83 UTM、メートル単位の投影法ですのでエラーが出ます
SELECT 'SRID=26910;POINT(-123 34)'::geography;

ERROR:  Only lon/lat coordinate systems are supported in geography.
LINE 1: SELECT 'SRID=26910;POINT(-123 34)'::geography;

ジオグラフィ型は、レイキャビクとシアトル-ロンドン間の飛行機の大圏コースとの間の、球面上での本当に最も近い距離を計算します。

大圏コースマップ ジオメトリ型は、平面の世界地図上にプロットされたレイキャビクとシアトル-ロンドン間の直線とのデカルト距離という意味のない値を出します。結果の名目上の単位は「度」ですが、点間の本当の角度差にあっていませんので、「度」と言うこと自体不正確です。

ジオグラフィ型によって、経度緯度座標でデータを格納できるようになりましたが、ジオグラフィで定義されている関数が、ジオメトリより少ないのと、実行にCPU時間がかかる、というところが犠牲になっています。

選択した型が、期待する領域から出ないことを、ジオメトリ型にして使用する条件とすべきです。使用するデータは地球全体か、大陸か、州か、自治体か?

ジオグラフィとジオメトリ間のサポート状況の比較については「PostGIS関数対応マトリクス」をご覧下さい。ジオグラフィ関数の簡潔なリストと説明については「PostGISジオグラフィ対応関数」をご覧下さい。

spatial_ref_sysテーブルは、PostGISに含まれるもので、3000以上の空間参照系と、座標変換/投影変換を行うのに必要な詳細情報との一覧を持つ、OGC準拠のデータベーステーブルです。

PostGISのspatial_ref_sysテーブルには、projライブラリで使われる3000以上の一般に使われる空間参照系定義がありますが、全てを持っているわけではなく、proj4の構築に慣れているならカスタム投影を定義することができます。ほとんどの空間参照系は地域限定のもので、想定されている範囲の外で使うと意味が無いことに注意して下さい。

PostGISのコアセットに入っていない空間参照系を探すための素晴らしい資料がhttp://spatialreference.org/にあります。

よく共通的に使われる空間参照系は次の通りです (訳注: 日本では状況が異なります)。4326 - WGS 84 経度緯度、4269 - NAD 83 経度緯度、3395 - WGS 84 メルカトル図法、2163 - 米国ナショナルアトラス正積図法、NAD83とWGS84のUTM空間参照系 - UTMゾーンは計測に最も理想的なものの一つですが6度 (訳注: 経度)の領域しかカバーしません。

まざまな米国の州の平面空間参照系 (メートルまたはフィート単位) - 通常は州ごとに一つか二つあります。メートル単位のもののほとんどがPostGISのコアセットに入っていますが、多数のフィート単位のものやESRIが作ったものについてはspatialreference.orgから取得して下さい。

対象領域がどのUTMゾーンになるかを決めるには、utmzone PostGIS plpgsql helper functionを参照して下さい。

SPATIAL_REF_SYSテーブル定義は次の通りです。

CREATE TABLE spatial_ref_sys (
  srid       INTEGER NOT NULL PRIMARY KEY,
  auth_name  VARCHAR(256),
  auth_srid  INTEGER,
  srtext     VARCHAR(2048),
  proj4text  VARCHAR(2048)
)

SPATIAL_REF_SYSのカラムは次の通りです。

GEOMETRY_COLUMNSは、データベースのシステムカタログから読み取るビューです。構造は次の通りです。

\d geometry_columns

View "public.geometry_columns"
      Column       |          Type          | Modifiers
-------------------+------------------------+-----------
 f_table_catalog   | character varying(256) |
 f_table_schema    | character varying(256) |
 f_table_name      | character varying(256) |
 f_geometry_column | character varying(256) |
 coord_dimension   | integer                |
 srid              | integer                |
 type              | character varying(30)  |

カラムの意味は次の通りです。

空間データを持つテーブルの生成は、1段階でできます。2次元ラインストリングでWGS84経度緯度のジオメトリカラムを持つroadsテーブルの生成の例を次に示します。

CREATE TABLE ROADS (ID serial, ROAD_NAME text, geom geometry(LINESTRING,4326) );

次の、3次元ラインストリングを追加する例で示す通り、標準的なALTER TABLEコマンドを使ってカラムを追加できます。

ALTER TABLE roads ADD COLUMN geom2 geometry(LINESTRINGZ,4326);

これが必要になる事例に、SQLビューとバルクインサートの二つがあります。バルクインサートの場合には、カラムに制約を与えるか、ALTER TABLEを実行することで、geometry_columnsテーブル内の登録を訂正することができます。ビューの場合には、CAST演算を使用します。カラムが型修飾子に基づく場合には、生成処理によって正しく登録されるので、何も行う必要がありません。ジオメトリに適用する空間関数を持たないビューも、基礎となるテーブルのジオメトリカラムと同じように登録されます。

-- 次のようなビューがあるとします
CREATE VIEW public.vwmytablemercator AS
        SELECT gid, ST_Transform(geom, 3395) As geom, f_name
        FROM public.mytable;

-- 正しく登録するには、
-- ジオメトリをキャストします。
--
DROP VIEW public.vwmytablemercator;
CREATE VIEW  public.vwmytablemercator AS
        SELECT gid, ST_Transform(geom, 3395)::geometry(Geometry, 3395) As geom, f_name
        FROM public.mytable;

-- ジオメトリタイプが確実に2次元ポリゴンだと知っているなら
-- 次のようにできます。
DROP VIEW public.vwmytablemercator;
CREATE VIEW  public.vwmytablemercator AS
        SELECT gid, ST_Transform(geom,3395)::geometry(Polygon, 3395) As geom, f_name
        FROM public.mytable;

-- 次のように、バルクインサートで派生テーブルを生成したとしましょう
SELECT poi.gid, poi.geom, citybounds.city_name
INTO myschema.my_special_pois
FROM poi INNER JOIN citybounds ON ST_Intersects(citybounds.geom, poi.geom);

-- 新しいテーブルに2次元インデックスを作ります
CREATE INDEX idx_myschema_myspecialpois_geom_gist
  ON myschema.my_special_pois USING gist(geom);

-- ポイントが3次元ポイントであったり、XYMポイントであったりした場合には、
-- 次のように、2次元インデックスでなくN次元インデックスを作ることになるかも
-- 知れません。
CREATE INDEX my_special_pois_geom_gist_nd
        ON my_special_pois USING gist(geom gist_geometry_ops_nd);

-- 新しいテーブルのジオメトリカラムをgeometry_columnsに手動登録するには、
-- 次のようにします。
-- カラムを型修飾子ベースにするために、基礎となるテーブル構造も変更することに
-- 注意して下さい。
SELECT populate_geometry_columns('myschema.my_special_pois'::regclass);

-- PostGIS 2.0を使っていて、何らかの理由で古い制約をもとにした定義を行う
-- (派生テーブルが同じタイプやSRIDを持たないといった場合)ことが必要な場合には、
-- 新しい任意変数use_typemodをfalseにします。
SELECT populate_geometry_columns('myschema.my_special_pois'::regclass, false); 

古い制約を基にした手法は現在も対応していますが、制約を基にしたジオメトリカラムで直接的にビューで使われている場合は、型修飾子のようには正しくgeometry_columnsに登録されません。次の例では、型修飾子を使ったカラム定義と、制約に基づくカラムの定義とを行っています。

CREATE TABLE pois_ny(gid SERIAL PRIMARY KEY, poi_name text, cat text, geom geometry(POINT,4326));
SELECT AddGeometryColumn('pois_ny', 'geom_2160', 2160, 'POINT', 2, false);

psqlで次を実行します。

\d pois_ny;

型修飾子と制約に基づくのとでは異なった定義になっているのが見えます。

Table "public.pois_ny"
  Column   |         Type          |                       Modifiers

-----------+-----------------------+------------------------------------------------------
 gid       | integer               | not null default nextval('pois_ny_gid_seq'::regclass)
 poi_name  | text                  |
 cat       | character varying(20) |
 geom      | geometry(Point,4326)  |
 geom_2160 | geometry              |
Indexes:
    "pois_ny_pkey" PRIMARY KEY, btree (gid)
Check constraints:
    "enforce_dims_geom_2160" CHECK (st_ndims(geom_2160) = 2)
    "enforce_geotype_geom_2160" CHECK (geometrytype(geom_2160) = 'POINT'::text
        OR geom_2160 IS NULL)
    "enforce_srid_geom_2160" CHECK (st_srid(geom_2160) = 2160)

geometry_columnsでは、両方とも正しく登録されています。

SELECT f_table_name, f_geometry_column, srid, type
        FROM geometry_columns
        WHERE f_table_name = 'pois_ny';

f_table_name | f_geometry_column | srid | type
-------------+-------------------+------+-------
pois_ny      | geom              | 4326 | POINT
pois_ny      | geom_2160         | 2160 | POINT

しかし、次のようにビューを作ろうとします。

CREATE VIEW vw_pois_ny_parks AS
SELECT *
  FROM pois_ny
  WHERE cat='park';

SELECT f_table_name, f_geometry_column, srid, type
        FROM geometry_columns
        WHERE f_table_name = 'vw_pois_ny_parks';

型修飾子によるgeomのビューカラムは正しく登録されますが、制約に基づくものは正しく登録されません。

f_table_name   | f_geometry_column | srid |   type
------------------+-------------------+------+----------
 vw_pois_ny_parks | geom              | 4326 | POINT
 vw_pois_ny_parks | geom_2160         |    0 | GEOMETRY

これは、将来的にPostGISの版で変更されるかもしれませんが、今のところは、制約に基づくビューカラムを正しく登録させるには、次のようにします。

DROP VIEW vw_pois_ny_parks;
CREATE VIEW vw_pois_ny_parks AS
SELECT gid, poi_name, cat,
  geom,
  geom_2160::geometry(POINT,2160) As geom_2160
  FROM pois_ny
  WHERE cat = 'park';
SELECT f_table_name, f_geometry_column, srid, type
        FROM geometry_columns
        WHERE f_table_name = 'vw_pois_ny_parks';

f_table_name   | f_geometry_column | srid | type
------------------+-------------------+------+-------
 vw_pois_ny_parks | geom              | 4326 | POINT
 vw_pois_ny_parks | geom_2160         | 2160 | POINT

PostGISはOpen Geospatial Consortium (OGC)のOpenGIS仕様に準拠しています。多くのPostGISメソッドは、操作対象のジオメトリが単純かつ妥当であることが求められます (正確に言うとそう仮定します)。たとえば、ポリゴンの外に穴があるようなものの面積を計算したり、単純でない境界線を持つポリゴンを作ったりするのは、意味がありません。

OGS仕様に沿うと、単純なジオメトリとは、自己インタセクトや自己接触のような、異常な幾何点を持たないことです。主に0次元または1次元のジオメトリ ([MULTI]POINT, [MULTI]LINESTRING))に適用します。 他方、ジオメトリの妥当性は、主に2次元ジオメトリ ([MULTI]POLYGON)に適用し、妥当なポリゴンを特徴づける位置指定子の集合を定義します。個々のジオメトリクラスには、単純性と妥当性をさらに詳細に述べる特定の条件があります。

POINTは0次元ジオメトリオブジェクトとして常に単純です。

MULTIPOINTは、二つの座標値 (POINT)が同じでないなら単純です。

LINESTRINGは、2度同じPOINTを通らない (終点は除きます。この場合は線型環と呼ばれ、さらに言うと閉じていると思われます)なら単純です。

(a)	(b)
(c)	(d)

(a)と(c)は単純なLINESTRINGです。(b)と(d)は単純ではありません。

MULTILINESTRINGは、 全ての要素が単純で、かつ任意の2要素のインタセクトが要素の境界であるPOINTでだけ発生する場合に限って単純です。

(e)

(f)

(g)

(e)と(f)は単純なMULTILINESTRINGです。(g)は単純ではありません。

定義からPOLYGONは常に単純です。バウンダリ内の環 (外環と内環からなる)のうち二つがクロスしていないなら妥当です。POLYGONの境界は、POINTとインタセクトするかも知れませんが、接点にしかなりません (すなわち線上にない)。POLYGONはカットラインまたはスパイクを持たなくても良く、内環は外環の中に完全に含まれていなければなりません。

(h)	(i)	(j)
(k)	(l)	(m)

(h)と(i)は妥当なPOLYGONです。(jからm)は単一のPOLYGONとしては表現できませんが、(j)と(m)は妥当なMULTIPOLYGONとして表現できます。

MULTIPOLYGONは、全ての要素が妥当で、二つのポリゴン要素について内側がインタセクトしていない場合は妥当です。ポリゴン要素の任意の二つの境界は接触してもよいですが、有限な数のPOINTでなければなりません。

(n)

(o)

(p)

(n)と(o)は妥当でないMULTIPOLYGONです。(p)は妥当です。

GEOSライブラリを使って実装されている関数のほとんどは、ジオメトリがOpenGISシンプルフィーチャー仕様で定義されているように妥当であると仮定しています。ジオメトリが単純であるか、また妥当であるか、のチェックとしてST_IsSimple()とST_IsValid()が使えます。

-- 一般的に、線フィーチャーの妥当性のチェックは
-- 常にTRUEを返すので意味がありません
-- しかし、この例では、PostGISがOGCのIsValidの定義を拡張して
-- *一意な頂点*が2より少ないラインストリングについてFALSEを
-- 返すようにしています
gisdb=# SELECT
   ST_IsValid('LINESTRING(0 0, 1 1)'),
   ST_IsValid('LINESTRING(0 0, 0 0, 0 0)');

 st_isvalid | st_isvalid
------------+-----------
      t     |     f

デフォルトでは、PostGISはジオメトリ入力に関するこの妥当性チェックを適用しません。複雑なジオメトリの妥当性のチェックはCPU時間を多く必要とするためです。データソースが信用できない場合は、手動でこのチェックを強制するための制約を付けることができます。

ALTER TABLE mytable
  ADD CONSTRAINT geometry_valid_check
        CHECK (ST_IsValid(the_geom));

妥当な入力ジオメトリでPostGIS関数を呼んだのに「GEOS Intersection()がエラーを投げました!」というようなエラーメッセージに遭遇したなら、PostGISまたは使用しているライブラリの中のエラーを発見しました。PostGIS開発者に連絡するべきです。PostGIS関数が妥当である入力ジオメトリから妥当でないジオメトリが返る場合も同じです。

	厳格にOGCジオメトリに準拠すると、Z値やM値を持てません。ST_IsValid()は高次を考慮に入れません。AddGeometryColumn()を実行するとジオメトリの次元をチェックする制約が加わるので、そこで2を指定すれば十分です。

代表的な空間述語 (ST_Intersects, ST_Contains, ST_Crosses, ST_Touches, ...)は、求められる空間フィルタを適切に提供するには十分でないことが時々あります。

たとえば、道路網を表現する線型のデータセットがあるとします。ビジネスルールを無視しているかも知れませんが、点で交差するだけでなく線上で交差する道路区間を全て判別することがGIS解析者の仕事となるかも知れません。この場合、 ST_Crossesでは重要な空間フィルタとして十分ではありません。線型のフィーチャーでは、点でクロス (訳注: インタセクトより厳しい条件です)している場合のみtrueが返ります。 空間的にインタセクトしている (ST_Intersects)と判別された二つの区間の組み合わせについて、実際のインタセクト (ST_Intersection)を取り、 インタセクトのST_GeometryTypeが''LINESTRING'であるかを見る ([MULTI]POINTや[MULTI]LINESTRING等からなるGEOMETRYCOLLECTIONが帰ってくる場合にしっかり対応するため)、という2段階の方法で解くことはできます。 よりエレガントかつ速い解法が本当に望ましいでしょう。

二つ目の[理論的な]例として、GIS解析者が全ての湖の境界に線でインタセクトする波止場やドックの全ての位置を特定しようとするとします。ここで、波止場の一端だけ陸にあるとします。言いかえると、波止場が湖の中にあるが完全に中に入りきってはいなくて、湖と線でインタセクトして、波止場の一方の端が完全に湖に入っていて、かつもう一方の端点が境界線上にあります。解析者は空間述語の組み合わせで、非常に望まれるフィーチャーを引き出す必要があるかも知れません。 ST_Contains(lake, wharf) = TRUE ST_ContainsProperly(lake, wharf) = FALSE ST_GeometryType(ST_Intersection(wharf, lake)) = 'LINESTRING' ST_NumGeometries(ST_Multi(ST_Intersection(ST_Boundary(wharf), ST_Boundary(lake)))) = 1 …(複雑なので以下略)

ここで、Dimensionally Extended 9 Intersection Mode、略してDE-9IMを見てみましょう。

4.3.6.1. 理論

OpenGIS Simple Features Implementation Specification for SQLによると「二つのジオメトリの比較の基本的なアプローチは、二つのジオメトリの内部、境界、外部のインタセクションの比較と、『インタセクション行列』の要素に基づく2ジオメトリの関係の分類です」。

境界 (Boundary)

ジオメトリの境界は一つ低い次元のジオメトリの集合です。0次元のPOINTの境界は空集合です。LINESTRINGの境界は二つの端点です。 POLYGONの境界は外環と内環を形成する線です。

内部 (Interior)

ジオメトリの内部は境界を取り去った際に残るジオメトリです。POINTの内部はPOINT自身です。LINESTRINGの内部は二つの端点の間の実際の点の集合です。POLYGONの内部はポリゴンの内側の範囲であるサーフェスです。

外部 (Exterior)

ジオメトリの外部は、内部と境界を除いた全領域のサーフェスです。

ジオメトリaがあり、aの内部、境界、外部をそれぞれI(a)、B(a)、E(a)とします。数学的な行列表現は次のようになります。

	内部 (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) )

ここで、dim(a)はaの次元で、ST_Dimensionで規定されますが、{0,1,2,T,F,*}の値域を持ちます。

• 0 => 点

• 1 => 線

• 2 => 面

• T => {0,1,2}

• F => 空集合

• * => 何でも良い

二つのオーバラップするポリゴンについて可視化すると、次のようになります。

	内部 (Interior) 境界 (Boundary) 外部 (Exterior) 内部 (Interior) dim(...) = 2 dim(...) = 1 dim(...) = 2 境界 (Boundary) dim(...) = 1 dim(...) = 0 dim(...) = 1 外部 (Exterior) dim(...) = 2 dim(...) = 1 dim(...) = 2

左から右、上から下に読むと、次元行列は'212101212'と表現されます。

一つ目の例である、2線が線上でインタセクトする場合の関係行列は'1*1***1**'となります。

-- 線上でクロスする道路区間の判別
SELECT a.id
FROM roads a, roads b
WHERE a.id != b.id
AND a.geom && b.geom
AND ST_Relate(a.geom, b.geom, '1*1***1**');

二つ目の例である、一部が湖の水涯線上にある波止場についての関係行列は'102101FF2'となります。

-- 一部が湖の水涯線上にある波止場の判別
SELECT a.lake_id, b.wharf_id
FROM lakes a, wharfs b
WHERE a.geom && b.geom
AND ST_Relate(a.geom, b.geom, '102101FF2');

詳細情報ついては、次のページをご覧ください。

• OpenGIS Simple Features Implementation Specification for SQL (1.1版, 2.1.13.2節)

• Dimensionally Extended Nine-Intersection Model (DE-9IM)

• GeoTools: Point Set Theory and the DE-9IM Matrix

• Hui XiongさんによるEncyclopedia of GIS

データをテキスト表現に変換できるなら、フォーマットされたSQLを使うのがデータをPostGISに入れる最も簡単な方法です。Oracleや他のSQLデータベースを使うように、SQL端末モニタにSQLの"INSERT"ステートメントで一杯になった大きなテキストファイルをパイプで送ることで、大量のデータをロードできます。

データアップロードファイル (たとえばroads.sql)は次のようになるでしょう。

BEGIN;
INSERT INTO roads (road_id, roads_geom, road_name)
  VALUES (1,'LINESTRING(191232 243118,191108 243242)','Jeff Rd');
INSERT INTO roads (road_id, roads_geom, road_name)
  VALUES (2,'LINESTRING(189141 244158,189265 244817)','Geordie Rd');
INSERT INTO roads (road_id, roads_geom, road_name)
  VALUES (3,'LINESTRING(192783 228138,192612 229814)','Paul St');
INSERT INTO roads (road_id, roads_geom, road_name)
  VALUES (4,'LINESTRING(189412 252431,189631 259122)','Graeme Ave');
INSERT INTO roads (road_id, roads_geom, road_name)
  VALUES (5,'LINESTRING(190131 224148,190871 228134)','Phil Tce');
INSERT INTO roads (road_id, roads_geom, road_name)
  VALUES (6,'LINESTRING(198231 263418,198213 268322)','Dave Cres');
COMMIT;

データファイルは、次に示す"psql"というSQL端末モニタを使って、簡単にPostgreSQLにパイプで送ることができます。

psql -d [データベース名] -f roads.sql

shp2pgsqlデータローダは、ESRIシェープファイルをPostGIS/PostgreSQLデータベースに、ジオメトリまたはジオグラフィとして挿入するための適切なSQLに変換します。ローダには、次に示すコマンドラインフラグによって区別される、いくつかの操作モードがあります。

また、コマンドラインローダのほとんどのオプションに対応するshp2pgsql-guiグラフィカルユーザインタフェースがあります。一度限りのロードしか行わずスクリプト化しない場合やPostGISに不慣れな方にとって便利になるかもしれません。これはpgAdminIIIのプラグインとして構築することもできます。

ローダを使って入力ファイルを生成してアップロードするセッション例は次の通りです。

# shp2pgsql -c -D -s 4269 -i -I shaperoads.shp myschema.roadstable > roads.sql
# psql -d roadsdb -f roads.sql

変換とアップロードはUNIXのパイプを使うと一回で実行できます。

# shp2pgsql shaperoads.shp myschema.roadstable | psql -d roadsdb

データベースからデータを引き出す最も直接的な手段は、次のように、SQLのSELECTクエリを使って返ってくるレコードとカラムの数を減らし、結果のカラムを可読テキストファイルにダンプするやり方です。

db=# SELECT road_id, ST_AsText(road_geom) AS geom, road_name FROM roads;

road_id | geom                                    | road_name
--------+-----------------------------------------+-----------
          1 | LINESTRING(191232 243118,191108 243242) | Jeff Rd
          2 | LINESTRING(189141 244158,189265 244817) | Geordie Rd
          3 | LINESTRING(192783 228138,192612 229814) | Paul St
          4 | LINESTRING(189412 252431,189631 259122) | Graeme Ave
          5 | LINESTRING(190131 224148,190871 228134) | Phil Tce
          6 | LINESTRING(198231 263418,198213 268322) | Dave Cres
          7 | LINESTRING(218421 284121,224123 241231) | Chris Way
(6 rows)

しかし、返ってくる結果の数を削るために、なんらかの制限をかけることが重要となるときがあるでしょう。属性ベースの制限の場合、非空間テーブルで使う通常の文法と同じSQLを使うだけです。空間ベースの制限の場合、次の演算子が使用可能であり、便利です。

次に、これらの演算子をクエリで使うことができます。SQLコマンドラインからジオメトリとボックスの指定を行うときは、明示的に文字列表現をジオメトリに変換しなければならないことに注意して下さい。たとえば、次のようになります。ただし312は架空の空間参照系番号で、ここでのデータに合致しています。

SELECT road_id, road_name
  FROM roads
  WHERE roads_geom='SRID=312;LINESTRING(191232 243118,191108 243242)'::geometry;

上のクエリは"ROADS_GEOM"テーブルから、その値と等価である単一のレコードを返します。

道路がポリゴンで定義した面を通過するかどうかをチェックするには次のようにします。

SELECT road_id, road_name
FROM roads
WHERE ST_Intersects(roads_geom, 'SRID=312;POLYGON((...))');

最も一般的な空間クエリは「フレームベース」のクエリでしょう。これは、表示するためのデータの価値のある「マップフレーム」を取得するために、データブラウザやウェブマッパのようなクライアントソフトウェアに使われます。

"&&"演算子を使うとき、比較フィーチャーをBOX3DかGEOMETRYかに指定することができます。ただし、GEOMETRYを指定すると、それのバウンディングボックスが比較に使われます。

次に示すクエリのように、フレームにBOX3Dオブジェクトを使います。

SELECT ST_AsText(roads_geom) AS geom
FROM roads
WHERE
  roads_geom && ST_MakeEnvelope(191232, 243117,191232, 243119,312);

エンベロープの投影を指定するためにSRID 312を使っていることに注意して下さい。

pgsql2shpテーブルダンパは、データベースに直接接続して、テーブル (あるいはクエリによって定義されたもの)をシェープファイルに変換するものです。基本的な文法は次の通りです。

pgsql2shp [<オプション>] <database> [<スキーマ>.]<table>

pgsql2shp [<オプション>] <データベース> <クエリ>

コマンドラインオプションは次の通りです。

GiSTは「汎用的な検索木 (Generalized Search Tree)」の意味で、インデックスの一般化された形式です。GISインデックスに加えて、GiSTは通常のB木インデックスに従わない全ての種類の不規則なデータ構造 (整数配列, スペクトラルデータ等)の検索速度を向上させるために使います。

ひとたびGISデータテーブルが数千行を超えたら、空間検索の速度向上のためインデックスを構築したくなるでしょう (これは属性検索でない場合です。属性でしたら通常のインデックスを属性フィールドに追加します)。

GiSTインデックスをジオメトリカラムに追加するための文は次の通りです。

CREATE INDEX [インデックス名] ON [テーブル名] USING GIST ( [ジオメトリカラム名] ); 

上の文では常に2次元インデックスを構築します。n次元インデックスをジオメトリ型で使うには、次の文でインデックスを生成できます。

CREATE INDEX [インデックス名] ON [テーブル名] USING GIST ( [ジオメトリカラム名] gist_geometry_ops_nd);

空間インデックスの構築は、計算量を集中させて行われます。また、この時には、テーブルへの書き込みアクセスがブロックされます。そのため、本番システムではより遅いCONCURRENTLYを選択するかも知れません。次のようにします。

CREATE INDEX CONCURRENTLY [インデックス名] ON [テーブル名] USING GIST ( [ジオメトリカラム名] ); 

インデックス構築後に、時々PostgreSQLにテーブルの統計情報を集めさせると助かります。クエリプランの最適化に使われます。

VACUUM ANALYZE [テーブル名] [(カラム名)];

BRINは"Block Range Index"のことで、PostgreSQL 9.5で導入された一般性のあるインデックスの形式です。BRINは不可逆インデックスで、主に読み取り効率と書き込み効率とを妥協させるために使われます。主要な目標は、非常に大きなテーブルの取り扱いで、そのテーブルのために、テーブル内の物理的な位置との自然な相関があるカラムを一部に持ちます。GISインデックスだけでなく、BRINは規則的だったりそうでない多種多様なデータ構造 (整数、配列等)における検索速度の向上のために使われます。

ひとたびGISデータテーブルが数千行を超えたら、空間検索の速度向上のためインデックスを構築したくなるでしょう (これは属性検索でない場合です。属性でしたら通常のインデックスを属性フィールドに追加します)。GiSTインデックスは、サイズがデータベースに割り当てられているRAM容量を超えずに、ストレージサイズが十分であり、かつ書き込みワークロードのペナルティに余裕がある限りは、本当に高性能です。その他の場合にはBRINインデックスが他の選択肢として考えられます。

BRINインデックスの原理は、テーブルブロックの集合の全ての行に含まれる全てのジオメトリを余すところなく取り込んだバウンディングボックスだけを格納することです。明らかに、このインデックス手法は、ブロックレンジの結果バウンディングボックスが相互に排他的になるように、データが物理的に並ぶ場合にだけ効率が良くなります。結果インデックスは本当に小さいのですが、多くの場合で、GiSTインデックスより効率が落ちます。

BRINインデックスの構築はGiSTインデックスの構築より少ない資源集約で済みます。BRINインデックスがGiSTインデックスが求めるよりも十分の一以下であるのが本当に一般的です。BRINインデックスは、多数のテーブルブロックに対して一つのバウンディングボックスを格納するだけなので、この手の他のインデックスよりもディスク使用量が数千分の一程度で済むのが一般的です。

レンジ内で要約するブロック数を選択できます。この数字を減らすと、インデックスは大きくなりますが、効率向上の助けになる可能性があります。

BRINインデックスを「ジオメトリ」カラムに追加するための文は次の通りです。

CREATE INDEX [インデックス名] ON [テーブル名] USING BRIN ( [ジオメトリカラム名] ); 

上の文では常に2次元インデックスを構築します。3次元インデックスは次の文で生成できます。

CREATE INDEX [インデックス名] ON [テーブル名] USING BRIN ( [ジオメトリカラム名] brin_geometry_inclusion_ops_3d);

また、4次元演算子クラスを使う4次元インデックスを使うこともできます。

CREATE INDEX [インデックス名] ON [テーブル名] USING BRIN ( [ジオメトリカラム名] brin_geometry_inclusion_ops_4d);

上の書式では、デフォルトの数値やレンジ内ブロック数を使っていますが、値は128です。ブロック数を指定するには、レンジ内の要約をすることになります。次の文のようにして作ります。

CREATE INDEX [インデックス名] ON [テーブル名] USING BRIN ( [ジオメトリカラム名] ) WITH (pages_per_range = [数字]); 

また、BRINインデックスは、多数の行で一つのインデックス値を格納することを心に留めておいて下さい。テーブルに違う次元のジオメトリを格納する場合には、インデックスの効率が悪くなります。この効率欠落を回避するには、格納したジオメトリの次元数の最小値となる演算子クラスを選択します。

「ジオグラフィ」型もまたBRINインデックスに対応しています。BRINインデックスを「ジオグラフィ」カラムに構築するための文は次の通りです。

CREATE INDEX [インデックス名] ON [テーブル名] USING BRIN ( [ジオグラフィカラム名] ); 

上の書式では常に回転楕円対面上の地理空間オブジェクトの2次元インデックスを構築します。

現在のところは「包括対応」だけをここで考えています。これは、&&、~と@演算子だけが2次元の場合に使われることを意味します (「ジオメトリ」「ジオグラフィ」の両方)。&&&演算子は3次元ジオメトリで使えます。しばらくはKNN検索に対応しません。

SP-GiSTは、「空間分割された一般探索木」を表します。四分木、k次元木、基数木 (トライ木)のような分割探索木に対応するインデックスの総称的な形式です。このデータ構造の一般的な機能は、検索空間を反復して分割することですが、分割は等しいサイズである必要はありません。SP-GiSTは、GISインデックスだけでなく、電話回線のルーティングや、IPルーティング、部分文字列検索等といった、様々な種類のデータを探索する速度の向上に使われます。

GiSTインデックスのためのケースですので、空間オブジェクトを覆うバウンディングボックスを保存するという意味で、SP-GiSTインデックスは不可逆です。SP-GiSTインデックスは、GiSTインデックスの選択肢の一つとして考えることができます。能率試験によって、SP-GiSTインデックスは,、多数のオーバラップするオブジェクトがある「スパゲッティデータ」と呼ばれる状態のときに、特に有利であることが分かりました。

一度GISデータテーブルが数千行を超えると、データの空間探索の速度向上にSP-GiSTインデックスを使うと良いかも知れません。「ジオメトリ」カラムにSP-GiSTインデックスを構築するための文は次の通りです。

CREATE INDEX [インデックス名] ON [テーブル名] USING SPGIST ( [ジオメトリカラム] ); 

上の文では、2次元インデックスを構築します。ジオメトリ型の3次元インデックスは、次のように、3次元演算子クラスを使用して生成します。

CREATE INDEX [インデックス名] ON [テーブル名] USING SPGIST ([ジオメトリカラム] spgist_geometry_ops_3d);

空間インデックスの構築は、計算量を集中させて行われます。また、この時には、テーブルへの書き込みアクセスがブロックされます。そのため、本番システムでは、より遅いCONCURRENTLYを選択するかも知れません。次のようにします。

CREATE INDEX CONCURRENTLY [インデックス名] ON [テーブル名] USING SPGIST ( [ジオメトリカラム] ); 

インデックス構築後に、時々PostgreSQLにテーブルの統計情報を集めさせると助かります。クエリプランの最適化に使われます。

VACUUM ANALYZE [テーブル名] [(カラム名)];

SP-GiSTインデックスは次の演算子を含むクエリの実行速度を向上させられます。

現時点ではkNN探索に対応していません。

通常、インデックスは見えないところでデータアクセスの速度向上を行います。すなわち、ひとたびインデックスが構築されたら、クエリプランナは透過的に、クエリプランの速度を向上させるためにインデックス情報を使うべき時を判断します。残念なことに、PostgreSQLクエリプランナは、GiSTインデックスの使用について十分に最適化できず、時々、検索で空間インデックスを使用すべきなのに、テーブル全体を順に走査することがあります。

空間インデックスが使用されていない (または属性インデックスがその問題のために使用されていない)場合、次の二つのことができます。

クエリを作成するとき、&&のようなバウンディングボックスを基準とした演算子によってのみGiST空間インデックスの利点が出てくることだけは覚えておくことが重要です。ST_Distance()のような関数では演算の最適化を行うためにインデックスを使うことができません。たとえば、次のクエリでは、大きなテーブルでは本当に遅くなります。

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

このクエリは、geom_tableにおける (100000, 200000)の点から距離が100単位以内にある全てのジオメトリを選択します。このクエリは、テーブル内にあるそれぞれの点と指定した点との距離を計算する、すなわち、それぞれの行で一つのST_Distance()計算を行うため、遅くなります。1ステップでインデックスを利用するST_DWithin関数を使って、求められる距離計算の量を減らすことで回避できます。次のようにします。

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

このクエリは、同じジオメトリを選択しますが、より効果的な方法で行われます。the_geomにGiSTインデックスがあると仮定すると、クエリプランナはST_Distance()関数の結果を計算する前にインデックスを使って行を減らせると認識します。 &&演算子で使われるST_MakeEnvelopeジオメトリは、元のポイントに中央寄せした200単位の正方形です。これは「クエリボックス」です。&&演算子は、結果セットを「クエリボックス」にオーバラップするバウンディングボックスを持つジオメトリだけに素早く減らすためにインデックスを使います。「クエリボックス」がジオメトリテーブル全体の範囲より十分に小さいと仮定すると、行われなければならない距離計算の量は劇的に減少します。

本節の例では、線型の道、ポリゴンの自治体境界、の二つのテーブルを使います。テーブルの定義をしまします。bc_roadsについては次の通りです。

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

bc_municipalityテーブルの定義については次の通りです。

Column     | Type              | Description
-----------+-------------------+-------------------
gid        | integer           | Unique ID
code       | integer           | Unique ID
name       | character varying | City / Town Name
the_geom   | geometry          | Location Geometry (Polygon)

SELECT sum(ST_Length(the_geom))/1000 AS km_roads FROM bc_roads;

km_roads
------------------
70842.1243039643
(1 row)

SELECT
  ST_Area(the_geom)/10000 AS hectares
FROM bc_municipality
WHERE name = 'PRINCE GEORGE';

hectares
------------------
32657.9103824927
(1 row)

SELECT
  name,
  ST_Area(the_geom)/10000 AS hectares
FROM
  bc_municipality
ORDER BY hectares DESC
LIMIT 1;

name           | hectares
---------------+-----------------
TUMBLER RIDGE  | 155020.02556131
(1 row)

SELECT
  m.name,
  sum(ST_Length(r.the_geom))/1000 as roads_km
FROM
  bc_roads AS r,
  bc_municipality AS m
WHERE
  ST_Contains(m.the_geom, r.the_geom)
GROUP BY m.name
ORDER BY roads_km;

name                        | roads_km
----------------------------+------------------
SURREY                      | 1539.47553551242
VANCOUVER                   | 1450.33093486576
LANGLEY DISTRICT            | 833.793392535662
BURNABY                     | 773.769091404338
PRINCE GEORGE               | 694.37554369147
...

CREATE TABLE pg_roads as
SELECT
  ST_Intersection(r.the_geom, m.the_geom) AS intersection_geom,
  ST_Length(r.the_geom) AS rd_orig_length,
  r.*
FROM
  bc_roads AS r,
  bc_municipality AS m
WHERE
  m.name = 'PRINCE GEORGE'
        AND ST_Intersects(r.the_geom, m.the_geom);

SELECT
  sum(ST_Length(r.the_geom))/1000 AS kilometers
FROM
  bc_roads r,
  bc_municipality m
WHERE
        r.name = 'Douglas St'
        AND m.name = 'VICTORIA'
        AND ST_Intersects(m.the_geom, r.the_geom);

kilometers
------------------
4.89151904172838
(1 row)

SELECT gid, name, ST_Area(the_geom) AS area
FROM bc_municipality
WHERE ST_NRings(the_geom) > 1
ORDER BY area DESC LIMIT 1;

gid  | name         | area
-----+--------------+------------------
12   | SPALLUMCHEEN | 257374619.430216
(1 row)