GDALさん…

ラスタをごにょごにょできるフリーのライブラリ+ユーティリティがGDAL。読みは「ぐだる」の模様。
これまでGDALで遊んでこなかったのですが、gdal添付ツールで遊んだところ、gdaldem で 陰影図ができたりと面白かったと ( http://www.gdal.org/gdal_utilities.html 参照)。
で、PostGIS 2.0 から標準添付になる PostGIS raster に蓄えといたら、あまり深く考えなくても任意箇所を抜き出してくれるだろう、と考えたと。

叩き込むまで

まず、PostGIS rasterを有効にするには、postgis.sql と spatial_ref_sys.sql を実行して、PostGISを叩き込む。そのうえで rtpostgis.sql を実行する。
その後、ラスタファイルをSQLに変換します。

raster2pgsql.py -r (入力ファイル) -t (テーブル名) -k 256x256 -s (ESPGコード) -I > (出力ファイル)

…だったと思う(うろおぼえ)。-k 256x256でタイル化して、1タイル1レコードで叩き込んでくれるようになっています。-I でインデクス作成。
これでとりあえずSQLができるので、見てみるものの、正直言って、何書いてるか分からん。末尾で ST_ConvexHull() の関数インデックスを作成していることだけ分かった。
もうちょっと踏み込んだ情報については、PostGISマニュアル(現時点ではSVN版)の"PostGIS Raster Frequently Asked Questions"を参照してください。コマンドとか書いてます。

GDAL 1.7 さん…

gdalinfo等で PostGIS raster を参照する場合は、ファイル名を次のようにします。

"PG:host=(ホスト) port=(ポート, 通常は 5432) dbname=(データベース名) user=(Postgresユーザ名) password=(ほげほげ) table=(テーブル名)"

としてやると、勝手に PostgreSQL に接続してくれます。これは便利だ。
しかし。
gdalinto を実行すると、データ無しなかんじになった。
で、pglog見ると、SELECTのWHERE条件式が ST_Contains([固定値BBOX], rast) とかなってる(rastはラスタデータを入れているカラム)。
…逆じゃねえかおい。
完全に[固定値BBOX]の内側になるrastは抜き出されるけど、かすってるのは抜き出されない。で、タイル化してるから余計にかすってるレコードばかりで、全滅したんだろうと。

GDAL 1.8 さん…

くしくもGDALのバージョンがあがったので、入れ替えてやってみた。
今度は、1.7のようなことなく、しっかりと抜き出してくれた。
そこで気をよくして、基盤地図情報5mメッシュの某県全体をダウンロードして、ラスタにして、PostGISに叩き込みました。ついでに、これ、MapServerに食わせてWMSにしてやれ、とか思いました。
で、やってみると、確かにそれっぽいデータを返してくれるのですが、異様に遅い。ていうか600*400程度(詳細覚えてない)のラスタ吐き出すのに、数十秒かかる。吐き出すラスタのサイズを考えると、ラスタ計算でここまで酷いわけがない。
まず疑ったのは、クエリプランナさんがインデクスを考慮に入れなかったのかも知れない、と。
ただ、1レコード 256x256 で、インデクスを考慮に入れないのは考えにくい。
とりあえず、ST_Containsで絞り込むクエリを実行させてみた。

EXPLAIN SELECT rid FROM ... WHERE ST_Contains(ST_SetSRID(BOX2D(...)));

そしたら、瞬殺だったんですよ。
これはプランナさんのせいじゃない、と思いました。ご、ごめんなさいプランナさん。
で、まあとりあえずpglogを見てみた、どんなクエリが出されてるんだ…?
…全レコードなめとるがな。
ST_MetaData() とか ST_BandMetaData() で、メタデータを読み取って、何か(ごめん読みきれてない)をしているのですが、このクエリには ST_Contains() とか、インデックスを活かした WHERE 条件が付いていません。そりゃ遅くもなる。
GDALのソースを当たってみましたが、やはり自動ではインデクスを活かした絞込みをしていませんでした。
これは…どうすりゃいいんだ…。ていうか、どうもできん。
とりあえず、改善されるのを待つ方向に。
まずはなによりインデクスを活かしてくれるようにしていただければさいわいです。