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 SPATIALTABLE ( 
                  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 SPATIALTABLE ( 
                  THE_GEOM, 
                  THE_NAME 
        ) 
        VALUES ( 
                  GeomFromEWKT('SRID=312;POINTM(-126.4 45.32 15)'), 
                  'A Place' 
        )

The "canonical forms" of a PostgreSQL type are the representations you get with a simple query (without any function call) and the one which is guaranteed to be accepted with a simple insert, update or copy. For the postgis 'geometry' type these are:

- 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)
        

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 VARCHAR(256) NOT NULL, 
  F_TABLE_SCHEMA VARCHAR(256) NOT NULL, 
  F_TABLE_NAME 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 int4, PARK_NAME varchar(128), PARK_DATE date, PARK_TYPE varchar(2) );
SELECT AddGeometryColumn('parks', 'park_geom', 128, 'MULTIPOLYGON', 2 );

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

CREATE TABLE roads ( ROAD_ID int4, ROAD_NAME varchar(128) ); 
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()がエラーを投げました!」や「JTS Intersection()がエラーを投げました!」というようなメッセージに遭遇したら、それはたぶん、PostGISまたは使用しているライブラリの中のエラーを発見しました。PostGIS開発者に連絡するべきです。PostGIS関数が妥当である入力ジオメトリから妥当でないジオメトリが返る場合も同じです。

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

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

BEGIN;
INSERT INTO ROADS_GEOM (ID,GEOM,NAME ) VALUES (1,GeomFromText('LINESTRING(191232 243118,191108 243242)',-1),'Jeff Rd'); 
INSERT INTO ROADS_GEOM (ID,GEOM,NAME ) VALUES (2,GeomFromText('LINESTRING(189141 244158,189265 244817)',-1),'Geordie Rd'); 
INSERT INTO ROADS_GEOM (ID,GEOM,NAME ) VALUES (3,GeomFromText('LINESTRING(192783 228138,192612 229814)',-1),'Paul St'); 
INSERT INTO ROADS_GEOM (ID,GEOM,NAME ) VALUES (4,GeomFromText('LINESTRING(189412 252431,189631 259122)',-1),'Graeme Ave'); 
INSERT INTO ROADS_GEOM (ID,GEOM,NAME ) VALUES (5,GeomFromText('LINESTRING(190131 224148,190871 228134)',-1),'Phil Tce'); 
INSERT INTO ROADS_GEOM (ID,GEOM,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 id, AsText(geom) AS geom, name FROM ROADS_GEOM; 
id | geom                                    | 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 
  ID, NAME 
FROM ROADS_GEOM 
WHERE 
  GEOM ~= GeomFromText('LINESTRING(191232 243118,191108 243242)',-1);

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

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

SELECT 
  ID, NAME 
FROM ROADS_GEOM 
WHERE 
  GEOM && GeomFromText('POLYGON((191232 243117,191232 243119,191234 243117,191232 243117))',-1);

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

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

SELECT 
  AsText(GEOM) AS GEOM 
FROM ROADS_GEOM 
WHERE 
  GEOM && GeomFromText('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 [indexname] ON [tablename] 
  USING GIST ( [geometryfield] GIST_GEOMETRY_OPS ); 

空間インデックスの構築は、計算量を集中させて行われます。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空間インデックスの利点が出てくることだけは覚えておくことが重要です。distance()のような関数では演算の最適化を行うためにインデックスを使うことができません。 たとえば、次のクエリでは、大きなテーブルでは本当に遅くなります。

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

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

SELECT the_geom FROM geom_table
WHERE the_geom && 'BOX3D(90900 190900, 100100 200100)'::box3d
  AND 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)

postgis=# SELECT sum(length(the_geom))/1000 AS km_roads FROM bc_roads;
     km_roads
------------------
 70842.1243039643
(1 row)

postgis=# SELECT area(the_geom)/10000 AS hectares FROM bc_municipality 
          WHERE name = 'PRINCE GEORGE';
     hectares
------------------
 32657.9103824927
(1 row) 

postgis=# SELECT name, area(the_geom)/10000 AS hectares 
          FROM bc_municipality 
          ORDER BY hectares DESC 
          LIMIT 1;
     name      |    hectares
---------------+-----------------
 TUMBLER RIDGE | 155020.02556131
(1 row)

postgis=# SELECT m.name, sum(length(r.the_geom))/1000 as roads_km 
          FROM bc_roads AS r,bc_municipality AS m 
          WHERE r.the_geom && m.the_geom 
          AND 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
 ...

postgis=# CREATE TABLE pg_roads as
          SELECT intersection(r.the_geom, m.the_geom) AS intersection_geom, 
                 length(r.the_geom) AS rd_orig_length, 
                 r.* 
          FROM bc_roads AS r, bc_municipality AS m 
          WHERE r.the_geom && m.the_geom 
          AND intersects(r.the_geom, m.the_geom) 
          AND m.name = 'PRINCE GEORGE';

postgis=# SELECT sum(length(r.the_geom))/1000 AS kilometers 
          FROM bc_roads r, bc_municipality m 
          WHERE r.the_geom && m.the_geom 
          AND r.name = 'Douglas St' 
          AND m.name = 'VICTORIA';
    kilometers
------------------
 4.89151904172838
(1 row)

postgis=# SELECT gid, name, area(the_geom) AS area 
          FROM bc_municipality 
          WHERE 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"
  SIZE 2
  COLOR 255 0 0
 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を使用する方法も示しています。

...

...