AddGeometryColumns 使わんでよくなったので感涙にむせんでるところ
ツイッターで http://blog.opengeo.org/2012/03/06/postgis-2-0-new-features-typmod/ が流れてきてました。
この文自体は読んでないのですが、PostGIS 2.0 から GEOMETRY型の Typmod になってるぜ、というものです。
AddGeometryColumns 使わなくていいようになった
Typmod というのは、つまるところ…
CREATE TABLE foo ( gid SERIAL PRIMARY KEY, geom GEOMETRY(Point, 4326) );
…でジオメトリカラムの生成が可能になるというもの。
以前は…
CREATE TABLE foo ( gid SEIRAL PRIMARY KEY ); SELECT AddGeometryColumn('foo','geom',4326,'POINT',2);
…と、後で追加を実行しなければならなかったのです。
AddGeometryColumn の引数の順序を覚えられなかった、アタマの悪い私にとっては、この時点で感激。
1.5 でも GEOGRAPHYの方は Typmod になってたりします。いやこのときもありがたい、GEOMETRYでもこうなってほしい、と思ったものです。
ALTER COLUMN でも活躍
さらに、ブログは続く。ALTER COLUMN が使える、と。
たとえば「空間参照系を EPSG:4326 にしてたけど EPSG:4612 にしたい」という場合などに役に立つ、と。
とりあえず1行入れてみます。
INSERT INTO foo (geom) SELECT 'SRID=4326;POINT(135 35)';
EWKTで見ると、空間参照系IDも出てくるので、空間参照系の確認に便利ですね。
db=# SELECT ST_AsEWKT(geom) FROM foo; st_asewkt
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- -
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
4326を4612にするのが「便利」かどうかは人によってとらえ方が違うと思います。
EPSG:4326 のフィーチャーを EPSG:4612 にして、geom にセットしようとすると、カラムの制限に引っかかります。
db=# UPDATE foo SET geom=ST_Transform(geom,4612); ERROR: Geometry SRID (4612) does not match column SRID (4326)
本来なら、カラムの制限を変更して、GEOMETRY_COLUMNSを変更して、geomの変換を行った結果をセットする、という3つの手続きを踏まないといけませんが、これがクエリひとつで可能です。
ALTER TABLE foo ALTER COLUMN geom TYPE GEOMETRY(point,4612) USING ST_Transform(geom,4612);
USINGの後に更新に使う関数を用意してやると、あわせて変換してくれます。ただし、PostgreSQL 8 以降でないと、USING句に対応していません。
db=# SELECT ST_AsEWKT(geom) FROM foo; st_asewkt
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- -
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
余談: AddGeometryColumns と GEOMETRY_COLUMNS
PostGIS 1.5 では、インストール時点で GEOMETRY_COLUMNS というテーブルが生成されました。これはカラムのメタデータを集めていて、PostGIS のクライアントは、このテーブルからカラムのジオメトリタイプなどを取得しています。QGIS もそう。
PostGIS内部では、PostgreSQLの機能を活用して、カラムに対して、SRIDの制限とかジオメトリタイプの制限とかを付けます。メタデータの方を見に行くよりPostgreSQLの機能使った方が安全性が上がるからだろうと思います。
1.5以前での AddGeometryColumns は、クライアントのための作業と内部のための作業を実行していました。つまり、カラムを生成して、カラムに制限を付け、GEOMETRY_COLUMNS にメタデータを入れる、という作業を実行していました。
PostGIS 2.0 では GEOMETRY_COLUMNS は、ビューに変更されました。このビューは、カラムの制限情報を読んで、これまでの GEOMETRY_COLUMNS テーブルと同じになるように加工したうえで返します。これで常にクライアント向けのメタ情報提供 (GEOMETRY_COLUMNS) と内部でのふるまい (カラムの制限) とで矛盾が生じないようになります。
ずっと前に、ジオメトリカラムを持つテーブルを削除する前に DropGeometryColumn でジオメトリカラムを削除しないと GEOMETRY_COLUMNS にゴミが残ってしまった、という経験を持っている人にとっては、「これからは DROP TABLE でばっさり切って問題ないんだぜ、いぇーい」となると思います。
MULTI- を取るのは完全には無理
shp2pgsqlで「ポリゴン」のシェープをSQLに変換した際、オプションで "-S" を指定しないと、MULTIPOLYGON になってしまう場合があります。これを POLYGON に変換できないか、というもの。
結論から言うと、完全には無理です。
CREATE TABLE mfoo ( gid SERIAL PRIMARY KEY, geom GEOMETRY(MULTIPOINT, 4326) ); INSERT INTO mfoo (geom) SELECT 'SRID=4326;MULTIPOINT(135 35, 134 35)';
で、MULTIPOINTを持つテーブルができました。
MULTI- から MULTI を取るのは ST_Dump です。これは、path と geom からなる複合型の集合で返ります。複数行になるというわけです。
db=# SELECT ST_AsEWKT( (ST_Dump(geom)).geom ) FROM mfoo; st_asewkt
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- -
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
こんなかんじ。
これを ALTER COLUMN TYPE で使えるかどうか、というと…
db=# ALTER TABLE mfoo ALTER geom TYPE geometry(point, 4326) USING ST_Dump(geom); ERROR: transform expression must not return a set
… PostgreSQL側で「変換式は集合を返したらいかん」と怒られました。
MULTI- を外すというのは、単純に外すのではなくて、要素ごとに分解して、それぞれで1行にしないといけないので、他のカラムの値を推測せよというのも無理な話ですから、あきらめざるを得ないです。
本当に MULTI- を外したいなら、別テーブルを作って INSERT でデータを作らないといけないでしょう。
ただし、全ての要素を取らなくてもいい、という場合には、ST_GeometryN を使うことで可能です。MULTI- の1つ目の要素を返させて、EWKT で出力する例を示します。
db=# SELECT ST_AsEWKT(ST_GeometryN(geom,1)) FROM mfoo; st_asewkt
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- -
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
これを ALTER COLUMN TYPE で使うことはできます。
ALTER TABLE mfoo ALTER geom TYPE geometry(point, 4326) USING ST_GeometryN(geom,1);
なお、2つ目以降はなくなります。
db=# SELECT gid, ST_AsEWKT(geom) FROM mfoo; gid | st_asewkt
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- -
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
PostGIS raster を試してみているが大量なファイルだとちょっとマズいのかも知れない
id:hogeman:20120307 とかを見て、奈良盆地は昔は池やってんでとかは昔に聞いたことあったけど奈良湖て言うんかー、という方向に向かうも踏みとどまる。
PostGIS raster についてやっておられてたので、マネしようとした次第。
ごそごそしてると広島県あたりの 10m DEM を発見。どう GDAL がサポートするラスタに変換したかは忘れた。とにかく 1125x750 の1バンドのラスタファイルがあった、と。
raster2pgsql -I -M -s 4612 -t 1125x750 (ラスタファイル名) ... (テーブル名) | psql
みたいなかんじ。
オプションについて -I は「インデクスを作る」 -M は「最後にVACUUM ANALYZEを実行」 -s 4612 は空間参照系の指定 -t 1125x750 はタイルファイルサイズ(ピクセル) です。
QGIS で見ようとする
QGISは 1.7.4 版を使用しています。
hogeman さんのブログで初めて知った "Load Postgis Raster to QGIS" (0.5.3版) プラグインを使用して、閲覧できましたが、一部つまづいたところがあったので、メモ。
- psycopg2 (たぶん PostgreSQL接続を行うPythonスクリプト) が必要。OSGeo4w の場合インストーラからインストール可。
- ダイアログが出て "(ユーザ名)@(ホスト名)" の「ユーザ名」の入力を求められたが、既にユーザ名は分かってるだろう、と空白にしてたら、全く接続できなかった。ダイアログで入力したユーザ名をそのまま使うらしい。
で、QGISにロードできたところで、プロパティを開いて、スタイルタブで、最小値と最大値(0と1000にしたけど高地ならまた違う)をユーザ定義して、コントラスト強調で「最大値・最小値まで引き伸ばす」にしてみると、おー、出た。
たくさんのファイルをロードしてみて大変なことになる
今度は、200ファイルぐらいを一気に raster2pgsql にかけて、データベースに叩き込んでみました。
raster2pgsql は、以前は Pythonスクリプトだったのですが、今はバイナリにしてて、スピードが段違いでした。
で、叩き込んでみたのを確かめようとしたところ、QGISちゃんが応答してくれないの…。
いちど止めて、あらためて、QGISの表示範囲を先ほどの1ファイルがフルに見える程度に絞ってロードを始めてみたところ、やはり応答してくれないの…。
PostgreSQL のログを見てると、何か一生懸命クエリを実行しています。
ひとつ抽出したのがこれです。
select rid, rast from public.dem where rast ~ st_setsrid(st_makebox2d(st_point(133.000000, 34.833333), st_point(133.125000, 34.916667)),4612)
"~" は ST_Covers(A,B) の演算子。クエリごとに box2d の範囲が違っています。
ロードできたが…
応答してくれない状態でブログ書いてたら、その脇でロードが終了しました。
で、なんとなく少しだけズームアウトしてみたら、また応答してくれなくなったorz