PostGIS 2.0.0マニュアル日本語訳

ほとんどのLinuxのインストールと同様に、最初のステップでは、ソースコードのビルドに使われる Makefile を生成します。これは、シェルスクリプトで行います。

./configure

パラメータを付けない場合には、このコマンドは自動で、PostGISのソースコードのビルドを行うのに必要なコンポーネントやライブラリをシステム上で探します。./configureとするのが一般的な使い方ですが、標準的でない位置に必要なライブラリやプログラムを置いてある場合のために、いくつかのパラメータを受け付けます。

次のリストで、共通して使われるパラメータを示します。 完全なリストについては、--helpまたは--help=shortパラメータを使って下さい。

	PostGISをSVNレポジトリから得る場合には、はじめに次のスクリプトを実行します。 ./autogen.sh このスクリプトによってconfigureスクリプトが生成されます。これはPostGISのインストールに関するカスタマイズに使われます。 PostGISをアーカイブファイルとして入手する場合には、configureが既に生成されているので./autogen.shは不要です。

Makefileが生成されたら、PostGISのビルドは、次のコマンドを実行するだけです。

make

出力の最後の行が "PostGIS was built successfully. Ready to install." と出れば終わりです。

PostGIS 1.4.0版からは、全ての関数に文書から生成されるコメントが付きます。これらのコメントを後ほどインストールするには、次のコマンドを実行しますが、docbookが必要です。アーカイブファイルからインストールする場合は、postgis_comments.sql, raster_comments.sql, topology_comments.sqlは、docフォルダにあるので、コメントを作成する必要はありません。

make comments

PostGIS 2.0で導入されました。早見表や学習中の方のハンドアウトに適しているHTMLチートシートを生成します。xsltprocが必要で、topology_cheatsheet.html, tiger_geocoder_cheatsheet.html, raster_cheatsheet.html, postgis_cheatsheet.htmlの4ファイルが生成されます。

HTMLとPDFのビルド済みのものはPostGIS / PostgreSQL Study Guidesにあります。

make cheatsheets

PostGIS EXTENSIONはPostgreSQL 9.1以上で使えます。extensionsフォルダに移動してmake installを実行するとビルドできます。この処理は将来的に、9.1を対象としてコンパイルする時に自動化され、ビルド処理の中に完全に統合されます。

ソースレポジトリからビルドしている場合は、関数の記述を最初にビルドする必要があります。

make comments

アーカイブファイルからのビルドの場合は、ビルド済みのものがあるので、コメントのビルドは必須ではありません。

PostgreSQL 9.1を対象にビルドしている場合は、extensionsは自動的にmake install処理の一部としてビルドするべきです。必要ならextensionsフォルダからビルドできますし、他のサーバで必要ならファイルの複製ができます。

cd extensions
cd postgis
make clean
make 
make install
cd ..
cd postgis_topology
make clean
make 
make install
	  

EXTENSIONファイルは、OSに関係なく、常に同じ版のPostGISと同じです。PostGISバイナリを既にインストールしている限りは、EXTENSIONファイルを上をあるOSから別のものに複写して大丈夫です。

開発用と異なる別のサーバでEXTENSIONを手動でインストールしたい場合は、次のファイルをEXTENSIONフォルダからPostgreSQLインストール先のPostgreSQL / share / extensionフォルダに複写します。通常のPostGISで、サーバに無い場合に必要になるバイナリと同じです。

実行すると、PgAdmin->extensionでpostgis, postgis_topologyが有効なEXTENSIONとして見えます。

psqlを使う場合は、次のクエリを実行してEXTENSIONがインストールされていることを確認できます。

SELECT name, default_version,installed_version 
FROM pg_available_extensions WHERE name LIKE 'postgis%' ;
      name       | default_version | installed_version
-----------------+-----------------+-------------------
postgis          | 2.0.0     | 2.0.0
postgis_topology | 2.0.0      | 

クエリを行ったデータベースにEXTENSIONがインストールされている場合は、installed_versionカラムに記載が見えます。 レコードブロックが無い場合は、PostGIS EXTENSIONがインストールされていないことになります。PgAdmin III 1.4以上では、データベースブラウザツリーのextensionsセクションで提供されていて、右クリックでアップグレードまたアンインストールできます。

有効なEXTENSIONがある場合、pgAdmin EXTENSIONインタフェースまたは次のSQLの実行によって選択したデータベースにPostGIS EXTENSIONをインストールできます。

CREATE EXTENSION postgis;
CREATE EXTENSION postgis_topology;

この素晴らしいEXTENSION機能を使わずに2.0.0をインストールした場合でもEXTENSIONベースに変更することができます。まずpostgis_upgrade_20_minor.sql, raster_upgrade_20_minor.sql, topology_upgrade_20_minor.sqlのアップグレードスクリプトを実行して最新版にアップグレードします。

ラスタ機能無しでPostGISをインストールした場合には、最初にラスタ機能をインストールする必要があります(rtpostgis.sqlを使います)。

それから、次のコマンドを実行して、個々のEXTENSIONの関数をEXTENSION形式に更新します。

CREATE EXTENSION postgis FROM unpackaged;
CREATE EXTENSION postgis_topology FROM unpackaged;

PostGISのテストを行うには、次のコマンドを実行します。

make check

このコマンドで、実際のPostgreSQLデータベースに対して生成したライブラリを使用した、様々なチェックとレグレッションテストを行います。

	PostgreSQL, GEOS または Proj4 を標準の位置にインストールしていない場合には、環境変数LD_LIBRARY_PATHに、ライブラリの位置を追加する必要があるかも知れません。

	現在のところmake checkは、チェックを行う際に 環境変数PATHとPGPORTによっています。コンフィギュレーションパラメータ--with-pgconfigを使って特定したPostgreSQLではありません。PATHを編集して、コンフィギュレーションの際に検出したPostgreSQLと一致するようにして下さい。もしくは、間もなく襲ってくる頭痛の準備をしておいて下さい。

成功した場合は、テストの出力は次のようなかんじになります。

     CUnit - A Unit testing framework for C - Version 2.1-0
	 http://cunit.sourceforge.net/


Suite: print_suite
  Test: test_lwprint_default_format ... passed
  Test: test_lwprint_format_orders ... passed
  Test: test_lwprint_optional_format ... passed
  Test: test_lwprint_oddball_formats ... passed
  Test: test_lwprint_bad_formats ... passed
Suite: Misc Suite
  Test: test_misc_force_2d ... passed
  Test: test_misc_simplify ... passed
  Test: test_misc_count_vertices ... passed
  Test: test_misc_area ... passed
  Test: test_misc_wkb ... passed
Suite: PointArray Suite
  Test: test_ptarray_append_point ... passed
  Test: test_ptarray_append_ptarray ... passed
Suite: PostGIS Computational Geometry Suite
  Test: test_lw_segment_side ... passed
  Test: test_lw_segment_intersects ... passed
  Test: test_lwline_crossing_short_lines ... passed
  Test: test_lwline_crossing_long_lines ... passed
  Test: test_lwline_crossing_bugs ... passed
  Test: test_lwpoint_set_ordinate ... passed
  Test: test_lwpoint_get_ordinate ... passed
  Test: test_point_interpolate ... passed
  Test: test_lwline_clip ... passed
  Test: test_lwline_clip_big ... passed
  Test: test_lwmline_clip ... passed
  Test: test_geohash_point ... passed
  Test: test_geohash_precision ... passed
  Test: test_geohash ... passed
  Test: test_isclosed ... passed
Suite: PostGIS Measures Suite
  Test: test_mindistance2d_tolerance ... passed
  Test: test_rect_tree_contains_point ... passed
  Test: test_rect_tree_intersects_tree ... passed
  Test: test_lwgeom_segmentize2d ... passed
Suite: WKT Out Suite
  Test: test_wkt_out_point ... passed
  Test: test_wkt_out_linestring ... passed
  Test: test_wkt_out_polygon ... passed
  Test: test_wkt_out_multipoint ... passed
  Test: test_wkt_out_multilinestring ... passed
:
:
--Run Summary: Type      Total     Ran  Passed  Failed
               suites       17      17     n/a       0
               tests       143     143     143       0
               asserts    1228    1228    1228       0


Creating spatial db postgis_reg
 Postgis 2.0.0SVN - 2011-01-11 15:33:37
   GEOS: 3.3.0-CAPI-1.7.0
   PROJ: Rel. 4.6.1, 21 August 2008

Running tests

 loader/Point.............. ok
 loader/PointM.............. ok
 loader/PointZ.............. ok
 loader/MultiPoint.............. ok
 loader/MultiPointM.............. ok
 loader/MultiPointZ.............. ok
 loader/Arc.............. ok
 loader/ArcM.............. ok
 loader/ArcZ.......... ok
 loader/Polygon.............. ok
 loader/PolygonM.............. ok
 loader/PolygonZ.............. ok
 regress. ok
 regress_index. ok
 regress_index_nulls. ok
 lwgeom_regress. ok
 regress_lrs. ok
 removepoint. ok
 setpoint. ok
 simplify. ok
 snaptogrid. ok
 affine. ok
 measures. ok
 long_xact. ok
 ctors. ok
 sql-mm-serialize. ok
 sql-mm-circularstring. ok
 sql-mm-compoundcurve. ok
 sql-mm-curvepoly. ok
 sql-mm-general. ok
 sql-mm-multicurve. ok
 sql-mm-multisurface. ok
 polyhedralsurface. ok
 out_geometry. ok
 out_geography. ok
 in_gml. ok
 in_kml. ok
 iscollection. ok
 regress_ogc. ok
 regress_ogc_cover. ok
 regress_ogc_prep. ok
 regress_bdpoly. ok
 regress_proj. ok
 dump. ok
 dumppoints. ok
 wmsservers_new. ok
 tickets. ok
 remove_repeated_points. ok
 split. ok
 relatematch. ok
 regress_buffer_params. ok
 hausdorff. ok
 clean. ok
 sharedpaths. ok
 snap. ok

Run tests: 55
Failed: 0

PostGISをインストールするには、次のコマンドを実行します。

make install

これにより、PostGISのインストールファイルが、--prefixパラメータで指定した、適切なサブディレクトリに複写されます。次に特筆すべきサブディレクトリを示します。

先にmake commentsを実行していた場合は、次のコマンドを実行すると、SQLファイルのインストール先にpostgis_comments.sql, raster_comments.sqlがインストールされます。

make comments-install

	postgis_comments.sql, raster_comments.sql, topology_comments.sqlは、xsltprocの外部依存ができたので、通常のビルドとインストールから切り離されました。

まず前述のPostGISのインストールを行います。

extrasフォルダが無い場合はhttp://www.postgis.org/download/postgis-2.0.0.tar.gzをダウンロードして下さい。

tar xvfz postgis-2.0.0.tar.gz

cd postgis-2.0.0/extras/tiger_geocoder/tiger_2010

tiger_loader.sqlのサーバ実行ファイルへのパス等を編集して下さい。

最初のTigerジオコーダのインストールで場合は、Windowsではcreate_geocode.batを、Linux/Unix/MacOS Xではcreate_geocode.shを、それぞれPostgreSQLの設定にあわせて編集します。コマンドラインから編集したスクリプトを実行します。このファイルを編集しない場合には、アイテムの位置が一般的なものになるだけです。Loader_Generate_Script実行後に生成されたスクリプトは編集できます。

データベースにtigerスキーマがあるか確認して下さい。データベースのsearch_pathに含まれているかを確認して、含まれていない場合は、

ALTER DATABASE geocoder SET search_path=public, tiger;

コマンドで追加して下さい。

住所正規化機能は、ややこしい住所を除いて、データなしでも大体動きます。次に示すテストを行って確認して下さい。

SELECT pprint_addy(normalize_address('202 East Fremont Street, Las Vegas, Nevada 89101')) As pretty_address;
pretty_address
---------------------------------------
202 E Fremont St, Las Vegas, NV 89101
			

2.0に含まれるTigerジオコーダがインストールされている場合には、どうしても必要な訂正がある際の臨時のアーカイブファイルからでも機能のアップグレードができます。

extrasフォルダが無い場合にはhttp://www.postgis.org/download/postgis-2.0.0.tar.gzをダウンロードします。

tar xvfz postgis-2.0.0.tar.gz

cd postgis-2.0.0/extras/tiger_geocoder/tiger_2010

Windowsの場合はupgrade_geocoder.batスクリプト、Linux/Unix/MacOS Xの場合はupgrade_geocoder.shスクリプトの位置を特定します。 PostGISデータベースの資格情報を持つように編集し、コマンドラインからスクリプトを実行します。

より詳細なデータロードの説明はextras/tiger_geocoder/tiger_2010/READMEにあります。ここでは一般的な手順に触れるだけです。

ロード処理によって、センサスウェブサイトからそれぞれ求める州のデータをダウンロードし、ファイルを展開し、各州をそれぞれの州テーブルにロードします。 州テーブルはtigerスキーマで定義されているテーブルを継承しているので、全てのデータにアクセスするのは、これらのテーブルへのクエリで十分です。 また、ある州の再読み込みを行う必要がある場合や、それ以上の州が不要である場合に、 Drop_State_Tables_Generate_Scriptを使って州テーブルの集合を削除します。

必要なデータをロードできるようにするため、次のツールが必要です。

データロードについては、求める州のデータをロードするスクリプトをプラットフォームに合わせて生成するLoader_Generate_Scriptを参照して下さい。小さい単位でインストールできることに注意して下さい。求める州全てを一度にロードする必要はありません。必要に応じてロードすることができます。

求める州のダウンロードが終わった後に、次を実行します。

SELECT install_missing_indexes();

as described in Install_Missing_Indexes

動作すべきように動作するかを試すために、ロードした州の住所のジオコードをGeocodeで実行してみます。

EXTENSIONを使ってインストールした場合は、EXTENSIONモデルでアップグレードしなければなりません。 古いSQLスクリプトを使ってインストールした場合は、SQLスクリプトでアップグレードすべきです。適切な方を参照して下さい。

2.9.1.1. 9.1より前またはEXTENSIONを使わないソフトアップグレード

PostGISをEXTENSIONを使わずにインストールした人向けです。EXTENSIONを使っていてこの方法を使うと、次のようなメッセージが現れます。

can't drop ... because postgis extension depends on it

コンパイル後にpostgis_upgrade*.sqlを探して下さい。PostGISの版にあったものをインストールします。たとえば、PostGIS 1.3から1.5に上げるには postgis_upgrade_13_to_15.sqlを使います。PostGIS 1.*から2.0に移動したり、2.*からr7409以前に落とす場合は、ハードアップグレードして下さい。

psql -f postgis_upgrade_20_minor.sql -d your_spatial_database

ラスタ機能とトポロジ機能についても同じ手続きです。それぞれrtpostgis_upgrade*.sql, topology_upgrade*.sqlになります。次のようにします。

psql -f rtpostgis_upgrade_20_minor.sql -d your_spatial_database

psql -f topology_upgrade_20_minor.sql -d your_spatial_database

	アップグレードのための特定の版のpostgis_upgrade*.sqlが見つからない場合は、非常に古い版を使っています。ハードアップグレードが必要です。

PostGIS_Full_Version関数によって、"procs need upgrade"メッセージを使ってこの種のアップグレードを実行する必要性についての情報が得られます。

2.9.1.2. 9.1以上でEXTENSIONを使ったソフトアップグレード

EXTENSIONを使ってPostGISをインストールした場合には、EXTENSIONを使ってアップグレードする必要があります。EXTENSIONを使ったマイナーアップグレードはかなり楽です。

ALTER EXTENSION postgis UPDATE TO "2.0.0";
ALTER EXTENSION postgis_topology UPDATE TO "2.0.0";

次のようなエラー通知が表示されることがあります。

No migration path defined for ... to 2.0.0

この場合は、データベースをバックアップして、「EXTENSIONを使った空間データベースの生成」に記述されているように新しいデータベースを生成し、バックアップを新しいデータベースにリストアしなければなりません。 postgis extensionが既にインストールされているとのメッセージを得るかも知れませんが、無視して問題ありません。

	PostGISをバージョン指定なしにインストールした場合には、 しばしばリストアの前のPostGIS EXTENSIONの再インストールをとばすことができます。 バックアップはCREATE EXTENSION postgisだけで、リストアの間に最新版になります。

ハードアップグレードとは、PostGISで利用可能なデータの完全なダンプ/リロードを意味します。PostGISオブジェクトの内部格納状態が変更される場合や、ソフトアップグレードができない場合に、ハードアップグレードが必要です。付録のRelease Notesに、版ごとについて、ダンプ/リロード(ハードアップグレード)の要否を記載しています。

ダンプ/リロード作業はpostgis_restore.plスクリプトが補助します。このスクリプトは、PostGIS(古いものを含む)に属する定義を全て飛ばすように注意します。また、重複シンボルエラーや非推奨オブジェクトの持越し無く、スキーマとデータをPostGISをインストールしたデータベースにリストアできます

Windows用に関する追加情報は Windows Hard upgradeにあります。

手続きは次の通りです。

エラーは次の場合に起こりえます。

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

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

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

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

bytea WKB = ST_AsBinary(geometry);
text WKT = ST_AsText(geometry);
geometry = ST_GeomFromWKB(bytea WKB, SRID);
geometry = ST_GeometryFromText(text 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情報が埋め込まれます。

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

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

bytea EWKB = ST_AsEWKB(geometry);
text EWKT = ST_AsEWKT(geometry);
geometry = ST_GeomFromEWKB(bytea EWKB);
geometry = 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"型の場合は次の通りです。

- Output
  - binary: EWKB
	ascii: HEXEWKB (EWKB in hex form)
- Input
  - binary: EWKB
	ascii: HEXEWKB|EWKT 

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

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

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

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

TSQL-MMの定義では、3dm、3dzと4dの座標を含みますが、SRID情報の埋め込みはできません。

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

	PostGIS 1.4より前では、曲線ポリゴンで複合曲線をサポートしていませんでしたが、PostGIS 1.4以降は曲線ポリゴンでの複合曲線の使用をサポートするようになりました。

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

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

新しいジオグラフィカラムはgeometry_columnsに登録されません。システムカタログを見るgeography_columnsという新しいビューに登録されるので、AddGeom... といった関数を使わずに、自動管理されます。

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

CREATE TABLEでジオグラフィカラムを持つテーブルを作ることができます。ジオメトリと違って、AddGeometryColumns()でメタデータにカラム情報を登録する処理を別に行う必要がなくなりました。

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

locationカラムはGEOGRAPHY型で、2つのオプション修飾子をサポートすることにご注意ください。ひとつは、そのカラムで使用できる形状と次元を限定する型修飾子です。もうひとつは、座標参照IDを特定の数に限定するSRID修飾子です。

型修飾子で許される値はPOINT, LINESTRING, POLYGON, MULTIPOINT, MULTILINESTRING, MULTIPOLYGONです。この修飾子ではZ, M, ZMの拡張子によって次元を限定します。たとえば、'LINESTRINGM'という値では、3次元で、第3次元が計測値(measure)となります。同じように見て 'POINTZM' は、4次元となります。

SRID修飾子には、現在は4326(WGS84)のみ許されるという制限があります。SRIDを指定したくない場合は0(未定義の回転楕円体)を使います。この場合の全ての計算は、WGS84とします。

将来的には、他のSRIDによって、WGS84以外の回転楕円体に関する計算ができるようにします。

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

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

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

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

GEOMETRYと同じ操作でインデクスを作成します。PostGISは、カラム型がGEOGRAPHYであるかを見て、GEOMETRYで使われる平面用インデクスの代わりに球面ベースのインデクスを作成します。

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

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

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

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

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

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

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

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

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

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

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

spatial_ref_sysテーブルは、PostGISに含まれるもので、3000以上の空間参照系を持つ、OGC準拠のデータベーステーブルです。これらの間での変換/投影変換に必要です。

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

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のカラムは次の通りです。

2.0.0より前のPostGISでは、直接編集可能なgeometry_columnsはテーブル、時々、実際のジオメトリカラムとの同期が取れていませんでした。PostGIS 2.0.0では、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 int4
	  	, ROAD_NAME varchar(25), geom geometry(LINESTRING,4326) );

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

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

後方互換のため、今でも、管理関数を使って空間テーブルを2段階で生成することもできます。

次はテーブルを作成して空間カラムを作る例です(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 );

うひとつ、ジェネリックな"geometry"型とSRID不明を示す0を使った例を挙げます。

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

AddGeometryColumn()アプローチでは、ジオメトリカラムを作成して、新しいカラムをgeometry_columnsテーブルに登録します。ソフトウェアでgeometry_columnsを使う場合には、クエリに必要なジオメトリカラムの全てがこのビューに登録されている必要があります。PostGIS 2.0からは、geometry_columnsは編集可能でなく、全てのジオメトリカラムは自動登録されます。

しかし、 カラムが生成時に特定のタイプで定義されなかった場合には、一般的なジオメトリカラムとして登録することがあります。

こういうことが発生して、AddGeometryColumnが使えない事例は、SQLビューとバルクインサートの2つがあります。これらの場合には、カラムに制約を与えることで、geometry_columnsテーブル内の登録を訂正することができます。 PostGIS 2.0以上では、カラムが型修飾に基づく場合には、生成処理によって正しく登録されるので、何も行う必要がありません。

-- 次のように作ったビューがあるとしましょう
CREATE VIEW  public.vwmytablemercator AS
	SELECT gid, ST_Transform(geom,3395) As geom, f_name
	FROM public.mytable;
	
-- PostGIS 2.0以上で正しく登録するには
-- ジオメトリへのキャストが必要です。
--
DROP VIEW public.vwmytablemercator;
CREATE VIEW  public.vwmytablemercator AS
	SELECT gid, ST_Transform(geom,3395)::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次元または3Mのポイントの場合、
-- 2次元インデクスのかわりにn次元インデクスを
-- 作成します
CREATE INDEX my_special_pois_geom_gist_nd 
	ON my_special_pois USING gist(geom gist_geometry_ops_nd);

-- 手動で新テーブルのジオメトリカラムをgeometry_columnsに登録するには、
-- このアプローチはPostGIS 2.0でもPostGIS 1.4以上でも動作します。
-- PostGIS 2.0ではテーブル構造をカラムをTypmodに基づくものに変更します。
-- PostGIS 2.0より前ではこれはビューを登録するのにも使えます。
SELECT populate_geometry_columns('myschema.my_special_pois'::regclass); 

-- PostGIS 2.0を使用していて、どんな理由にしても、
-- 古い制約を基にした定義の挙動が必要だとします
-- (全ての子テーブルが同じタイプとSRIDを持っているわけではない継承テーブルといった場合)
-- 新しい、省略可能引数のuse_typmodを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 varchar(20)
   , 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は、2つの座標値(POINT)が同じでないなら単純です。

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

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

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

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

(e)

(f)

(g)

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

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

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

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

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

(n)

(o)

(p)

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

GEOSライブラリを使って実装されている関数のほとんどは、ジオメトリがOpenGIS Simple Feature Specificationで定義されているように妥当であると仮定しています。ジオメトリが単純であるか、また妥当であるか、のチェックとして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()がエラーを投げました!"や"JTS Intersection()がエラーを投げました!"というようなメッセージに遭遇したら、それはたぶん、PostGISまたは使用しているライブラリの中のエラーを発見しました。PostGIS開発者に連絡するべきです。PostGIS関数が妥当である入力ジオメトリから妥当でないジオメトリが返る場合も同じです。

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

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

たとえば、道路網を表現する線型のデータセットがあるとします。ビジネスルールを無視しているかも知れませんが、点で交差するだけでなく線上で交差する道路区間を全て判別することがGIS解析者の仕事となるかも知れません。この場合、ST_Crossesでは重要な空間フィルタとして十分ではありません。線型のフィーチャーでは、点で交差している場合のみtrueが返ります。 空間的にインタセクトしている(ST_Intersects)と判別された2つの区間の組み合わせについて、実際のインタセクト(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 Modelまたは略してDE-9IMを見てみましょう。

4.3.6.1. Theory

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

境界(Boundary)

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

内部(Interior)

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

外部(Exterior)

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

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

i￡°	内部	境界	外部
内部	dim( I(a) ∩ I(b) )	dim( I(a) ∩ B(b) )	dim( I(a) ∩ E(b) )
境界	dim( B(a) ∩ I(b) )	dim( B(a) ∩ B(b) )	dim( B(a) ∩ E(b) )
外部	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 => 空集合

• * => 何でも良い

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

i￡°
	i￡° 内部 境界 外部 内部 dim(...) = 2 dim(...) = 1 dim(...) = 2 境界 dim(...) = 1 dim(...) = 0 dim(...) = 1 外部 dim(...) = 2 dim(...) = 1 dim(...) = 2

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

1つ目の例である、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**');

2つ目の例である、一部が湖の水涯線上にある波止場についての関係行列は'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 (version 1.1, section 2.1.13.2)

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

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

• Encyclopedia of GIS By Hui Xiong

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

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

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

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

psql -d [database] -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コマンドラインからジオメトリとボックスの特定を行うときは、"GeomFromText()"関数で、明示的に文字列表現をジオメトリに変換しなければならないことに注意して下さい。 たとえば、次のようになります(312は架空の空間参照系番号で、ここでのデータに合致しています)。

SELECT road_id, road_name
  FROM roads
  WHERE ST_OrderingEquals(roads_geom , ST_GeomFromText('LINESTRING(191232 243118,191108 243242)',312) ) ;

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

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

SELECT road_id, road_name
FROM roads
WHERE roads_geom && ST_GeomFromText('POLYGON((...))',312);

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

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

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

SRID 312を使っていますが、エンベロープの投影を指定してます。

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

pgsql2shp [<options>] <database> [<schema>.]<table>

pgsql2shp [<options>] <database> <query>

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

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

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

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

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

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

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

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

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

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

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

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

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

SELECT the_geom
FROM geom_table
WHERE ST_DWithin(the_geom,  ST_MakeEnvelope(90900, 190900, 100100, 200100,312), 100)

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

	ふるまいの変更
	PostGIS 1.3.0では、ST_DisjointとST_Relateの注目すべき例外がありますが、ほとんどのジオメトリ関係関数は暗黙的なバウンディングボックスオーバラップ演算子を含んでいます。

本節の例では、線形の道、ポリゴンの自治体境界、の2つのテーブルを使います。テーブルの定義をしまします。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_Contains(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)

raster2pgsqlは、GDALがサポートするラスタ書式をPostGISラスタテーブルにロードするのに適切なSQLにするバイナリファイルです。ラスタのオーバビューの生成だけでなく、ラスタファイルのフォルダのロードも可能です。

raster2pgsqlは、ほとんどの場合、PostGISの一部としてコンパイルされます(GDALライブラリをコンパイルしていない場合)が、バイナリファイルによってサポートされるラスタタイプは、GDALでコンパイルされたのと同じです。raster2pgsqlがサポートするラスタタイプの一覧を得るには、-Gスイッチを使います。この一覧は、インストールした PostGIS が提供する ST_GDALDrivers と同じになるはずです。

	このツールの古い版では、Pythonスクリプトでした。実行ファイルに置き換えられています。Pytonスクリプトが必要な場合は、GDAL PostGIS Raster Driver Usageに、Pythonの例があります。raster2pgsql の Python スクリプトは、今後の PostGIS raster では動作しないかも知れませんし、サポートされませんんおで、ご注意ください。

使用例:

raster2pgsql raster_options_go_here raster_file someschema.sometable > out.sql

ローダを用いて入力ファイルを100x100のタイルで生成して、データベースにアップロードする例は、次の通りです。

	public.demelevationでなくdemelevationというようにスキーマ名を外すことができます。ラスタテーブルはデータベースまたユーザの指定するデフォルトのスキーマに生成されます。

raster2pgsql -s 4236 -I -C -M *.tif -F -t 100x100 public.demelevation > elev.sql
psql -d gisdb -f elev.sql

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

raster2pgsql -s 4236 -I -C -M *.tif -F -t 100x100 public.demelevation | psql -d gisdb

マサチューセッツ州平面のメートル単位の空中写真タイルをaerialという名前のスキーマにロードします。 元の画像と2, 4レベルのオーバビューのテーブルとを生成します。 データ格納にCOPYを使用し(データベースに仲介ファイルなくまっすぐ入ります)。 -eでトランザクションを指定しません(待たずにテーブルのデータを見たい場合には良いです)。ラスタを128x128ピクセルのタイルに分解してラスタ制約を適用します。INSERTモードでなくCOPYモードを使用します。-Fで、カラム名をタイル切り出し元ファイルのファイル名にします。

raster2pgsql -I -C -e -Y -F -s 26986 -t 128x128  -l 2,4 bostonaerials2008/*.jpg aerials.boston | psql -U postgres -d gisdb -h localhost -p 5432

-- サポートされているラスタタイプの一覧:
raster2pgsql -G

-Gコマンドの出力は次のようになります。

Available GDAL raster formats:
  Virtual Raster
  GeoTIFF
  National Imagery Transmission Format
  Raster Product Format TOC format
  ECRG TOC format
  Erdas Imagine Images (.img)
  CEOS SAR Image
  CEOS Image
  JAXA PALSAR Product Reader (Level 1.1/1.5)
  Ground-based SAR Applications Testbed File Format (.gff)
  ELAS
  Arc/Info Binary Grid
  Arc/Info ASCII Grid
  GRASS ASCII Grid
  SDTS Raster
  DTED Elevation Raster
  Portable Network Graphics
  JPEG JFIF
  In Memory Raster
  Japanese DEM (.mem)
  Graphics Interchange Format (.gif)
  Graphics Interchange Format (.gif)
  Envisat Image Format
  Maptech BSB Nautical Charts
  X11 PixMap Format
  MS Windows Device Independent Bitmap
  SPOT DIMAP
  AirSAR Polarimetric Image
  RadarSat 2 XML Product
  PCIDSK Database File
  PCRaster Raster File
  ILWIS Raster Map
  SGI Image File Format 1.0
  SRTMHGT File Format
  Leveller heightfield
  Terragen heightfield
  USGS Astrogeology ISIS cube (Version 3)
  USGS Astrogeology ISIS cube (Version 2)
  NASA Planetary Data System
  EarthWatch .TIL
  ERMapper .ers Labelled
  NOAA Polar Orbiter Level 1b Data Set
  FIT Image
  GRIdded Binary (.grb)
  Raster Matrix Format
  EUMETSAT Archive native (.nat)
  Idrisi Raster A.1
  Intergraph Raster
  Golden Software ASCII Grid (.grd)
  Golden Software Binary Grid (.grd)
  Golden Software 7 Binary Grid (.grd)
  COSAR Annotated Binary Matrix (TerraSAR-X)
  TerraSAR-X Product
  DRDC COASP SAR Processor Raster
  R Object Data Store
  Portable Pixmap Format (netpbm)
  USGS DOQ (Old Style)
  USGS DOQ (New Style)
  ENVI .hdr Labelled
  ESRI .hdr Labelled
  Generic Binary (.hdr Labelled)
  PCI .aux Labelled
  Vexcel MFF Raster
  Vexcel MFF2 (HKV) Raster
  Fuji BAS Scanner Image
  GSC Geogrid
  EOSAT FAST Format
  VTP .bt (Binary Terrain) 1.3 Format
  Erdas .LAN/.GIS
  Convair PolGASP
  Image Data and Analysis
  NLAPS Data Format
  Erdas Imagine Raw
  DIPEx
  FARSITE v.4 Landscape File (.lcp)
  NOAA Vertical Datum .GTX
  NADCON .los/.las Datum Grid Shift
  NTv2 Datum Grid Shift
  ACE2
  Snow Data Assimilation System
  Swedish Grid RIK (.rik)
  USGS Optional ASCII DEM (and CDED)
  GeoSoft Grid Exchange Format
  Northwood Numeric Grid Format .grd/.tab
  Northwood Classified Grid Format .grc/.tab
  ARC Digitized Raster Graphics
  Standard Raster Product (ASRP/USRP)
  Magellan topo (.blx)
  SAGA GIS Binary Grid (.sdat)
  Kml Super Overlay
  ASCII Gridded XYZ
  HF2/HFZ heightfield raster
  OziExplorer Image File
  USGS LULC Composite Theme Grid
  Arc/Info Export E00 GRID
  ZMap Plus Grid
  NOAA NGS Geoid Height Grids

データベース内でラスタやラスタテーブルを生成したい場合が多くあります。これを行うための関数が多数あります。一般的な手順は次の通りです。

raster_columnsは、ラスタタイプのデータベースにおける全てのラスタテーブルカラムのカタログです。テーブルの制約を使ったビューなので、他のデータベースのバックアップからラスタテーブルをリストアしたとしても、情報は常に矛盾がありません。raster_columnsカタログには次のカラムがあります。

ローダを使わずにテーブルを生成したり、ロード時に-Cフラグを忘れたりした場合には、事後にAddRasterConstraintsで制約を強制でき、raster_columnsカタログは、ラスタタイルの共通の情報を登録します。

raster_overviews は、オーバビューで使われるラスタテーブルカラムに関する情報のカタログで、オーバビューを用いる際に知っておくと便利な情報も持ちます。オーバビューテーブルはraster_columnsとraster_overviewsの両方のカタログに入れられます。オーバビューもラスタのひとつであるのは確かですが、より高い解像度テーブルの解像度を落としたカリカチュアになるという特殊な目的を満たすためでもあるからです。ラスタをロードする際に-lスイッチを使うと、オーバビューが主ラスタテーブルと一緒に生成されます。

オーバビューテーブルには、他のラスタテーブルと同じ制約と、オーバビュー特有の制約となる追加情報があります。

	raster_overviewsの情報はraster_columnsとは重複しません。raster_columnsにあるオーバビューテーブルに関する情報が必要な場合は、raster_overviewsとraster_columnsを結合すると、必要な情報の集合を完全に取得することができます。

オーバビューの主たる理由は次の2つです。

raster_overviewsカタログには、次の情報のカラムがあります。

本節では、PHP の PostgreSQL ドライバと ST_AsGDALRaster 等の関数を使って、HTML img タグに埋め込むことができる PHP リクエストストリームにラスタの 1, 2, 3 バンドを出力する方法を示します

サンプルクエリでは、 指定したWGS84バウンディングボックスにインタセクトするタイルを取って、 ST_Unionでインタセクトしたタイルを結合して全てのバンドを返し、ST_Transformでユーザ指定投影法に変換し、ST_AsPNGを使ってPNGで結果を出力するためのラスタ関数群全体をまとめる方法を示します。

後に示すスクリプトは、

http://mywebserver/test_raster.php?srid=2249

で、マサチューセッツ州平面(フィート単位)のラスタ画像を取得します。

<?php
/** test_raster.php コンテント **/
$conn_str ='dbname=mydb host=localhost port=5432 user=myuser password=mypwd';
$dbconn = pg_connect($conn_str);
header('Content-Type: image/png');  
/** 投影法が指定されたらそれを使います。
    指定されないなら、マサチューセッツ州平面(メートル) を使います。 **/
if (!empty( $_REQUEST['srid'] ) && is_numeric( $_REQUEST['srid']) ){
		$input_srid = intval($_REQUEST['srid']);
}
else { $input_srid = 26986; }
/** PostgreSQL 9.0以上ではbytea_outputの設定が必要でしょう。
    8.4 では不要です。 **/
$sql = "set bytea_output='escape';
SELECT ST_AsPNG(ST_Transform(
			ST_AddBand(ST_Union(rast,1), ARRAY[ST_Union(rast,2),ST_Union(rast,3)])
				,$input_srid) ) As new_rast
 FROM aerials.boston 
	WHERE 
	 ST_Intersects(rast, ST_Transform(ST_MakeEnvelope(-71.1217, 42.227, -71.1210, 42.218,4326),26986) )"; 
$result = pg_query($sql);
$row = pg_fetch_row($result);
pg_free_result($result);
if ($row === false) return;
echo pg_unescape_bytea($row[0]);
?>

本節では、Npgsql PostgreSQL .NETドライバとST_AsGDALRaster等の関数を使って、HTML imgタグに埋め込むことができるように、ラスタの 1, 2, 3 バンドを出力する方法を示します。

この例ではNpgsql .NET PostgreSQLドライバが必要です。最新版はhttp://npgsql.projects.postgresql.org/ あります。最新版をダウンロードして、ASP.NET の binフォルダに入れるだけでうまくいきます。

サンプルクエリでは、指定したWGS84バウンディングボックスにインタセクトするタイルを取り、ST_Unionでインタセクトしたタイルを結合して全てのバンドを返し、ST_Transformでユーザ指定投影法に変換し、ST_AsPNGを使ってPNGで結果を出力するためのラスタ関数群全体をまとめる方法を示します。

この例はC#で実装している点を除いては「ST_AsPNG を他の関数とあわせて使った PHP 出力例」と同じです。

後に示すスクリプトは、

http://mywebserver/TestRaster.ashx?srid=2249

で、マサチューセッツ州平面(フィート単位)のラスタ画像を取得します。

 -- web.config 接続文字列部 --
<connectionStrings>
    <add name="DSN" 
        connectionString="server=localhost;database=mydb;Port=5432;User Id=myuser;password=mypwd"/>
</connectionStrings>

// Code for TestRaster.ashx
<%@ WebHandler Language="C#" Class="TestRaster" %>
using System;
using System.Data;
using System.Web;
using Npgsql;

public class TestRaster : IHttpHandler
{
	public void ProcessRequest(HttpContext context)
	{
		
		context.Response.ContentType = "image/png";
		context.Response.BinaryWrite(GetResults(context));
		
	}

	public bool IsReusable {
		get { return false; }
	}

	public byte[] GetResults(HttpContext context)
	{
		byte[] result = null;
		NpgsqlCommand command;
		string sql = null;
		int input_srid = 26986;
        try {
		    using (NpgsqlConnection conn = new NpgsqlConnection(System.Configuration.ConfigurationManager.ConnectionStrings["DSN"].ConnectionString)) {
			    conn.Open();

                if (context.Request["srid"] != null)
                {
                    input_srid = Convert.ToInt32(context.Request["srid"]);  
                }
                sql = @"SELECT ST_AsPNG(
                            ST_Transform(
			                ST_AddBand(
                                ST_Union(rast,1), ARRAY[ST_Union(rast,2),ST_Union(rast,3)])
				                    ,:input_srid) ) As new_rast 
                        FROM aerials.boston 
	                        WHERE 
	                            ST_Intersects(rast, 
                                    ST_Transform(ST_MakeEnvelope(-71.1217, 42.227, -71.1210, 42.218,4326),26986) )";
			    command = new NpgsqlCommand(sql, conn);
                command.Parameters.Add(new NpgsqlParameter("input_srid", input_srid));
           
			
			    result = (byte[]) command.ExecuteScalar();
                conn.Close();
			}

		}
        catch (Exception ex)
        {
            result = null;
            context.Response.Write(ex.Message.Trim());
        }
		return result;
	}
}

ひとつのイメージを取って指定したファイルに出力する、Javaコンソールアプリケーションです。

最新のPostgreSQL JDBC ドライバはhttp://jdbc.postgresql.org/download.htmlから取得できます。

後に示すコードをコンパイルします。コマンドは次の通りです。

set env CLASSPATH .:..\postgresql-9.0-801.jdbc4.jar
javac SaveQueryImage.java
jar cfm SaveQueryImage.jar Manifest.txt *.class

次のようにコマンドラインから呼び出します。

java -jar SaveQueryImage.jar "SELECT ST_AsPNG(ST_AsRaster(ST_Buffer(ST_Point(1,5),10, 'quad_segs=2'),150, 150, '8BUI',100));" "test.png" 

 -- Manifest.txt --
Class-Path: postgresql-9.0-801.jdbc4.jar
Main-Class: SaveQueryImage

// Code for SaveQueryImage.java
import java.sql.Connection;
import java.sql.SQLException;
import java.sql.PreparedStatement;
import java.sql.ResultSet;
import java.io.*;

public class SaveQueryImage {
  public static void main(String[] argv) {
      System.out.println("Checking if Driver is registered with DriverManager.");
      
      try {
        //java.sql.DriverManager.registerDriver (new org.postgresql.Driver());
        Class.forName("org.postgresql.Driver");
      } 
      catch (ClassNotFoundException cnfe) {
        System.out.println("Couldn't find the driver!");
        cnfe.printStackTrace();
        System.exit(1);
      }
      
      Connection conn = null;
      
      try {
        conn = DriverManager.getConnection("jdbc:postgresql://localhost:5432/mydb","myuser", "mypwd");
        conn.setAutoCommit(false);

        PreparedStatement sGetImg = conn.prepareStatement(argv[0]);
		
        ResultSet rs = sGetImg.executeQuery();
		
		FileOutputStream fout;
		try
		{
			rs.next();
			/** ユーザが指定するファイルへの出力 **/
			fout = new FileOutputStream(new File(argv[1]) );
			fout.write(rs.getBytes(1));
			fout.close();
		}
		catch(Exception e)
		{
			System.out.println("Can't create file");
			e.printStackTrace();
		}
		
        rs.close();
		sGetImg.close();
        conn.close();
      } 
      catch (SQLException se) {
        System.out.println("Couldn't connect: print out a stack trace and exit.");
        se.printStackTrace();
        System.exit(1);
      }   
  }
}

これは、サーバディレクトリ内でレコードごとにファイルを生成するPythonストアド関数です。

// PLPython PostgreSQLストアドプロシージャです。
// PL/Pythonが必要です。
CREATE OR REPLACE FUNCTION write_file (param_bytes bytea, param_filepath text)
RETURNS text
AS $$
f = open(param_filepath, 'wb+')
f.write(param_bytes)
return param_filepath
$$ LANGUAGE plpythonu;

-- 5つの画像をPostgreSQLサーバに可変サイズで描きます。
-- PostgreSQLデーモンのアカウントにフォルダへの書き込み権限が必要ですので
-- ご注意ください。
-- 生成されたファイル名をエコーバックします。
 SELECT write_file(ST_AsPNG(
	ST_AsRaster(ST_Buffer(ST_Point(1,5),j*5, 'quad_segs=2'),150*j, 150*j, '8BUI',100)),
	 'C:/temp/slices'|| j || '.png')
	 FROM generate_series(1,5) As j;
	 
     write_file
---------------------
 C:/temp/slices1.png
 C:/temp/slices2.png
 C:/temp/slices3.png
 C:/temp/slices4.png
 C:/temp/slices5.png	 

悲しいことにPSQLには、組み込み機能を用いてバイナリを出力するのは簡単ではありません。これは、ちょっとしたハックです。Clever Trick Challenge -- Outputting bytea with psqlで提案されている概要に基づいています。レガシー気味になっているPostgreSQLラージオブジェクトのサポートに乗っているものです。最初の起動時に、PSQLはデータベースに接続します。

この方法はPythonの場合と違い、ローカル機にファイルが生成されます

SELECT oid, lowrite(lo_open(oid, 131072), png) As num_bytes
 FROM 
 ( VALUES (lo_create(0), 
   ST_AsPNG( (SELECT rast FROM aerials.boston WHERE rid=1) ) 
  ) ) As v(oid,png);
-- 次のような出力が得られます --
   oid   | num_bytes
---------+-----------
 2630819 |     74860
 
-- 続いて、oidに注意して'C:/temp/aerial_smap.png'を
-- ローカル機のファイルパスに置き換えたうえで
-- 次を実行します。
 \lo_export 2630819 'C:/temp/aerial_samp.png'
 
-- データベース上のラージオブジェクトストレージから
-- ファイルを削除します。
SELECT lo_unlink(2630819);
   			

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

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

MapServerは、他のPostgreSQLクライアントのように、libpqインタフェースを使って、PostGIS/PostgreSQLデータにアクセスします。よってMapServerはPostGISサーバにアクセスするネットワークを持つ計算機にインストールでき、PostGISをデータソースとして使用することができます。システム間の接続は速いほど良いです。

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

DATA "geom FROM (
  SELECT 
    table1.geom AS geom, 
    table1.gid AS gid, 
    table2.data AS data 
  FROM table1 
  LEFT JOIN table2 
  ON table1.id = table2.id
) AS new_table USING UNIQUE gid USING SRID=4326"

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

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

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

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

LAYER 
  CONNECTIONTYPE postgis 
  CONNECTION "user=theuser password=thepass dbname=thedb host=theserver" 
  PROCESSING "CLOSE_CONNECTION=DEFER"
  DATA "geom from roads"
  MINSCALE 100000 
  STATUS ON 
  TYPE LINE 
  FILTER "road_type = 'highway'" 
  CLASS 
    COLOR 0 0 0 
  END 
END 
LAYER 
  CONNECTIONTYPE postgis 
  CONNECTION "user=theuser password=thepass dbname=thedb host=theserver"
  PROCESSING "CLOSE_CONNECTION=DEFER"
  DATA "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  
    STYLE
      COLOR 0 0 0 
    END
  END 
END

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

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

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

LAYER 
  CONNECTIONTYPE postgis
  CONNECTION "user=theuser password=thepass dbname=thedb host=theserver" 
  DATA "geom FROM (SELECT roads.gid AS gid, roads.geom AS 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 gid USING SRID=4326" 
  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を使用する方法も示しています。

...

...

現版のPostgreSQL(8.0を含む)では、TOASTテーブルに従うクエリオプティマイザの弱さに苦しみます。 TOASTテーブルは、(長いテキスト、イメージ、多数の頂点を持つ複合ジオメトリといった)通常のデータページに適合しない、(データサイズという意味では)巨大な値を納めるための「拡張部屋」の一種です。詳細情報はthe PostgreSQL Documentation for TOASTをご覧ください。

(高解像度で全てのヨーロッパの国の境界を含むテーブルのような)大きなジオメトリがあるうえ、行がそう多くないテーブルを持つようになると、この問題が出てきます。テーブル自体は小さいのですが、多くのTOASTスペースを使います。例として、テーブル自体は概ね80行で3データページしか使わなくてもTOASTテーブルで8225ページを使うとします。

ここで、ジオメトリ演算子の&&を使って、ほとんどマッチしないようなバウンダリボックスを検索するクエリを出してみます。クエリオプティマイザにはテーブルは3ページ80行しかないように見えます。オプティマイザは、小さなテーブルを順に走査する方がインデクスを使うよりも早いと見積もります。そして、GiSTインデクスは無視すると決めます。通常なら、この見積もりは正しいです。しかし、この場合は&&演算子が全てのジオメトリをディスクから呼び出しでバウンディングボックスと比較しなければならなくなり、ゆえに、全てのTOASTページもまた呼び出す必要があります。

このバグに苦しむかどうかを見るには、PostgreSQLの"EXPLAIN ANALYZE"コマンドを使います。詳しい情報と技術に関する詳細については、postgres performance mailing list のスレッド(http://archives.postgresql.org/pgsql-performance/2005-02/msg00030.php)をご覧下さい。

PostgreSQLコミュニティでは、TOASTを意識したクエリ見積もりを作ることで、この問題を解決しようとしています。今のところは、二つの応急処置があります。

ひとつは、クエリプランナにインデクスの使用を強制することです。クエリを発行する前に"SET enable_seqscan TO off;"をサーバに送信します。これは基本的にクエリプランナに対して可能な限り順に走査することを避けるよう強制します。そのためGiSTインデクスを通常使うようになります。しかし、このフラグは接続するたびに設定しなければならず、他のケースにおいてはクエリプランナに誤った見積もりをさせることになるので、 "SET enable_seqscan TO on;"をクエリの後に送信すべきです。

もうひとつは、順に走査することをクエリプランナが考える程度に早くすることです。これは、バウンダリボックスの「キャッシュ」を行う追加カラムを作成し、このカラムにマッチさせるようにすることで達成することができます。ここでの例では次のようになります。

SELECT AddGeometryColumn('myschema','mytable','bbox','4326','GEOMETRY','2'); 
UPDATE mytable SET bbox = ST_Envelope(ST_Force_2d(the_geom));

そして、次のように、&&演算子をgeom_columnに対して行っていたものをbboxに変更します。

SELECT geom_column 
FROM mytable 
WHERE bbox && ST_SetSRID('BOX3D(0 0,1 1)'::box3d,4326);

もちろん、mytableの行を変更または追加したら、bboxを「同期」するようにしなければなりません。最もすっきりした方法はトリガです。もしくは、アプリケーションを変更してbboxカラムの現状を保持するか、テーブル更新後にいつもUPDATEクエリを実行するかでも対応できます。

これらの設定はpostgresql.conf内にあります。

checkpoint_segments

constraint_exclusion

shared_buffers

work_mem(並べ替えや複雑なクエリに使われるメモリ)

maintenance_work_mem (VACUUM, CREATE INDEX等で使われるメモリ)