第4章 PostGISを使う

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

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

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

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

bytea WKB = asBinary(geometry); 
text WKT = asText(geometry); 
geometry = GeomFromWKB(bytea WKB, SRID); 
geometry = GeometryFromText(text WKT, SRID);

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

INSERT INTO geotable ( the_geom, the_name )
  VALUES ( 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)の例は、次の通りです。

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

bytea EWKB = asEWKB(geometry); 
text EWKT = asEWKT(geometry); 
geometry = GeomFromEWKB(bytea EWKB); 
geometry = GeomFromEWKT(text EWKT);

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

INSERT INTO geotable ( the_geom, the_name ) 
  VALUES ( 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拡張はまだ完全にはサポートされていません。単純な曲線ジオメトリの例を次に示します。

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テーブルは、次のように定義されています。

CREATE TABLE geometry_columns ( 
  f_table_catalog    VARRCHAR(256) NOT NULL, 
  f_table_schema     VARCHAR(256) NOT NULL,
  f_table_nam        VARCHAR(256) NOT NULL, 
  f_geometry_column  VARCHAR(256) NOT NULL, 
  coord_dimension    INTEGER NOT NULL, 
  srid               INTEGER NOT NULL, 
  type               VARCHAR(30) NOT NULL 
)

カラムは次のとおりです。

空間データを持つテーブルを生成するには、次の通り二段階で行います。

次はテーブルを作成して空間カラムを作る例です (128というSRIDがあると仮定します)。

CREATE TABLE parks ( 
  park_id    INTEGER, 
  park_name  VARCHAR,
  park_date  DATE,
  park_type  VARCHAR
);
SELECT AddGeometryColumn('parks', 'park_geom', 128, 'MULTIPOLYGON', 2 );

もうひとつ、ジェネリックな「ジオメトリ」型とSRID不明を示す-1を使った例を挙げます。

CREATE TABLE roads ( 
  road_id INTEGER,
  road_name VARCHAR
);
SELECT AddGeometryColumn( 'roads', 'roads_geom', -1, 'GEOMETRY', 3 );

GEOSライブラリで実装される関数のほとんどが、ジオメトリがOpenGIS単純フィーチャー仕様の定義から見て妥当であることが仮定されています。ジオメトリの妥当性のチェックには、IsValid()関数を使います。例を次に示します。

gisdb=# select isvalid('LINESTRING(0 0, 1 1)'), 
        isvalid('LINESTRING(0 0,0 0)'); 

 isvalid | isvalid
---------+--------- 
       t |       f

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

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

妥当なジオメトリでPostGIS関数を呼んだ時に、"GEOS Intersection() threw an error!"または"JTS Intersection() threw an error!"のような、変なエラー・メッセージに遭遇するなら、たぶんPostGISか使用ライブラリでエラーが発生しています。PostGIS開発者と連絡をとるべきです。PostGIS関数が妥当な入力から無効なジオメトリーを返すならば、同じことがいえます。

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

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

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

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

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

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

-a, -c, -dおよび-pは互いに排他的であることに注意して下さい。

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

# shp2pgsql shaperoads myschema.roadstable > roads.sql 
# psql -d roadsdb -f roads.sql

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

# shp2pgsql shaperoads myschema.roadstable | psql -d roadsdb

データベースからデータを引き出す最もストレートな手段は、次のように、SQLのSELECTクエリを使ってカラムを可読なテキストファイルとして出力することです。

db=# SELECT road_id, 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コマンドラインからジオメトリとボックスの特定を行うときは、"GeomFromText()"関数で、明示的に文字列表現をジオメトリに変換しなければならないことに注意して下さい。たとえば、次のようになります。

SELECT road_id, road_name 
  FROM roads 
  WHERE roads_geom ~= GeomFromText('LINESTRING(191232 243118,191108 243242)',-1);

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

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

SELECT road_id, road_name 
FROM roads 
WHERE roads_geom && GeomFromText('POLYGON((...))',-1);

上のクエリでは、比較するためにポリゴンのバウンディングボックスを用いています。

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

SELECT AsText(roads_geom) AS geom 
FROM roads 
WHERE 
  roads_geom && SetSRID('BOX3D(191232 243117,191232 243119)'::box3d,-1);

ここで、BOX3Dの投影を指定するためにSRIDを使っていることに注意して下さい。SRIDを-1に設定しているのは、SRIDを指定しないということを示しています。

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

pgsql2shp [<オプション>] <データベース> [<スキーマ>.]<テーブル>

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

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

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

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

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

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

空間インデックスの構築は、計算量を集中させて行われます。100万行のテーブルで、300MHzのSolaris機ではGiSTインデックスの構築に概ね1時間かかりました。インデックスを構築したあとは、クエリプランの最適化に使うため、次のようにPostgreSQLにテーブル統計情報の収集をさせることが重要です。

VACUUM ANALYZE [テーブル名] [(カラム名)];
-- 次のクエリはPostgreSQL 7.4以前でのみ必要です
SELECT UPDATE_GEOMETRY_STATS( [テーブル名], [(カラム名)] );

GiSTインデックスはPostgreSQLのR木インデックスと比べて二つの利点を持っています。まず、GiSTインデックスは「NULLセーフ」、すなわちNULL値を含むインデックスカラムで利用できることです。次に、GiSTインデックスはGISオブジェクトがPostgreSQLで8Kのページサイズを超えるサイズを扱う際に重要な「不可逆」の概念を持っていることです。不可逆にすることによって、PostgreSQLは、インデックスにおけるオブジェクトの「重要な」部分、GISオブジェクトの場合にはバウンディングボックスになりますが、これのみを納めることができます。R木インデックスで8Kを超えるGISオブジェクトのインデックスを構築しようとすると、失敗します。

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

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

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

SELECT the_geom 
FROM geom_table 
WHERE ST_Distance(the_geom, GeomFromText('POINT(100000 200000)', -1)) < 100

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

SELECT the_geom 
FROM geom_table 
WHERE the_geom && 'BOX3D(90900 190900, 100100 200100)'::box3d 
  AND
ST_Distance(the_geom, GeomFromText('POINT(100000 200000)', -1)) < 100

このクエリは、同じジオメトリを選択しますが、より効果的な方法で行われます。the_geomにGiSTインデックスがあると仮定すると、クエリプランナは、distance()関数の結果を計算する前に行を減らすためにインデックスを使うことができると認識します。&&演算子で使われているBOX3Dジオメトリは、指定位置を中心とした一辺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 ST_Intersects(r.the_geom, m.the_geom)
  AND m.name = 'PRINCE GEORGE';

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

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)

MapServerでPostGISを使うには、MapServerのコンフィギュレーション方法についての知識が必要ですが、この文書の範囲外です。この節では、PostGIS特有の問題とコンフィギュレーション詳細について記載します。

MapServerでPostGISを使うには、次のものが必要です。

MapServerは、他のPostgreSQLクライアントのように、libpqを使ってPostGIS/PostgreSQLデータにアクセスします。これは、システムがPostgreSQLクライアントライブラリであるlibpqを持っている限りは、PostGISサーバへのネットワークアクセスを持つあらゆる機械にMapServerをインストールできるということを示しています。

USING疑似SQL節を使ってMapServerがより複雑なクエリの結果を理解できるようにするための情報を追加します。より詳しく言うと、ビューまたは副問い合わせが元テーブル (DATA定義で"FROM"の右にあるもの)として使われる時、MapServerが自動的に一意な識別子がそれぞれの行にあるか、また、SRIDがテーブルにあるかを判別するのは困難です。USING節によって、MapServerがこれらの情報を得ることができます。例を次に挙げます。

DATA "the_geom FROM (
  SELECT 
    table1.the_geom AS the_geom, 
    table1.oid AS oid, 
    table2.data AS data 
  FROM table1 
  LEFT JOIN table2 
  ON table1.id = table2.id
) AS new_table USING UNIQUE oid USING SRID=-1"

簡単な例から始めて、ステップアップしていきましょう。次のMapServerレイヤ定義を考えて下さい。

LAYER 
  CONNECTIONTYPE postgis 
  NAME "roads"
  CONNECTION "user=theuser password=thepass dbname=thedb host=theserver" 
  DATA "the_geom FROM roads" 
  STATUS ON 
  TYPE LINE 
  CLASS 
    COLOR 0 0 0 
  END 
END

このレイヤは"roads"テーブルにある道路ジオメトリの全部を黒線で表示するものです。

では、少なくとも1:100000にズームするまでは高速道路だけを表示したい、としましょう。次の二つのレイヤで、その効果が実現できます。

LAYER 
  CONNECTION "user=theuser password=thepass dbname=thedb host=theserver" 
  DATA "the_geom FROM roads"
  MINSCALE 100000 
  STATUS ON 
  TYPE LINE 
  FILTER "road_type = 'highway'" 
  CLASS 
    COLOR 0 0 0 
  END 
END 
LAYER 
  CONNECTION "user=theuser password=thepass dbname=thedb host=theserver"
  DATA "the_geom FROM roads" 
  MAXSCALE 100000 
  STATUS ON 
  TYPE LINE
  CLASSITEM road_type 
  CLASS 
    EXPRESSION "highway" 
    STYLE
      WIDTH 2 
      COLOR 255 0 0  
    END
  END 
  CLASS  
    COLOR 0 0 0 
  END 
END

一つ目のレイヤはスケールが1:100000以上であるときに使われ、道路タイプが"highway"である道路のみ黒線で表示されます。FILTERオプションによって、道路タイプが"highway"の場合のみ表示することになります。

二つ目のレイヤはスケールが1:100000未満である時に使われ、"highway"は赤い二重細線で表示され、他の道路は黒線で表示されます。

さて、MapServerの機能を使うだけで、二つのおもしろいことを実行しました。しかし、DATAのSQLステートメントは、単純なままです。道路名が (どういう理由かは知りませんが)他のテーブルに収められていて、それのデータを取得するためにテーブルを連結して、道路のラベルを取る必要がある、とします。

LAYER 
  CONNECTION "user=theuser password=thepass dbname=thedb host=theserver" 
  DATA "the_geom FROM (SELECT roads.oid AS oid, roads.the_geom AS the_geom, 
        road_names.name as name FROM roads LEFT JOIN road_names ON 
        roads.road_name_id = road_names.road_name_id) 
        AS named_roads USING UNIQUE oid USING SRID=-1" 
  MAXSCALE 20000 
  STATUS ON 
  TYPE ANNOTATION 
  LABELITEM name
  CLASS 
    LABEL 
      ANGLE auto 
      SIZE 8 
      COLOR 0 192 0 
      TYPE truetype 
      FONT arial
    END
  END 
END

このANNOTAIONレイヤでは、縮尺が1:20000以下のときに、全ての道路に緑色のラベルを表示します。また、この例は、DATA定義で、SQLのJOINを使用する方法も示しています。

...

...