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