QGISハンズオン中級（ラスタ編）

FOSS4G Advent Calendar 2012

この記事はFOSS4G Advent Calendar 2012の穴埋め用に用意しておいたものです。
30日目にして公開。
8月から下書き状態で眠っていた記事を、書き直したものなんですが。。。

ラスタ編だよ

半年ほどあいてしまいましたが、
FOSS4G Hokkaidoハンズオン（裏メニュー）
FOSS4G Hokkaidoハンズオン（裏メニュー）続き
の続きです。
ベクタ編に引き続きラスタ編もデータベースでどの辺を出来るかやってみます。
ラスタ編の資料はこちら。
https://dl.dropbox.com/u/1838363/QGIS_RAS.zip
ただしspatialiteはラスタに関しては解析機能がありませんので、PostgreSQL+PostGIS 2.0で行ってみます。

使った環境

この組み合わせで行っています。

Windows7(64bit)
PostgreSQL 9.2.1(64bit)
PostGIS 2.0.1

PostgreSQLのインストール方法は省略。Windowsであれば、インストーラをダウンロードしてくれば特に問題無いかと。
PostGIS 2.0.1のインストールは、PostgreSQLのインストールに続いてスタックビルダにて行っています。

準備

GUI使えるところはGUIで行きましょう。
PgAdmin?を起動して、ユーザーとデータベースを作っておきます。
ログインロール"waigania13"としておきました。

作ったユーザーをオーナーにしてデータベースを作成。

その際に、templateは"template_posgis_20"にしておくのをお忘れなく。

作業手順

さてここからが本番です。
QGISハンズオン中級(ラスタ編)の中身は適当なディレクトリに解凍しておいてください。PostgreSQLの各種コマンドはC:\Program Files\PostgreSQL\9.2\binにあるとして進めます。
ハンズオンで行なっている

ラスタの結合
ラスタの投影変換
ラスタ情報の確認
地形解析
ラスタのクリップ

辺りをなぞってみたいと思います。

10DEM Tiffファイルのインサート

raster2pgsqlで画像をデータベースに入れるためのSQL文を作成。demというテーブルに。QGIS_RAS\dataset\10DEMにある全てのtiffファイルを入れてしまいます。ラスタの結合はこれでOK。

raster2pgsql -c -I -s 4612 644130.tif dem > 644130.sql
raster2pgsql -a -s 4612 644131.tif dem > 644131.sql
raster2pgsql -a -s 4612 644140.tif dem > 644140.sql
raster2pgsql -a -s 4612 644141.tif dem > 644141.sql
raster2pgsql -a -s 4612 644150.tif dem > 644150.sql
raster2pgsql -a -s 4612 644151.tif dem > 644151.sql

なのですが、ここで問題が。出来たSQL文を実行すると、解像度が違うとか怒られます。
gdalinfoでタグの情報を見てみると、Pixel Sizeが微妙に違う(;_;)。644130.tifが、

Pixel Size = (0.000111111111111,-0.000111111110667)

644151.tifが、

Pixel Size = (0.000111111111111,-0.000111111112000)

精度的には無視しちゃっていい違いなのですが、ダメみたいですね。タグを書き換えたいとこですが、簡単な方法が思いつかず。

なので、ここは少々ずるをして、ハンズオン資料に入っている結合済みのラスタからスタートします。
QGIS_RAS\out_backup\6441.tifを使います。

raster2pgsql -c -I -s 4612 6441.tif dem > 6441.sql
psql -f 6441.sql -U waigania13 rastertest

psqlで少し中を見ておきましょう。

psql -U waigania13 rastertest

しておいて、\d demでカラムの表示、idとラスタデータのみ入っていることがわかります。

st_metadata()でメタデータも見ておきましょう。原点座標とピクセル数、解像度、SRIDといった情報が入っています。

QGISの準備

確認はQGISから行いましょう。
QGISからPostGISのラスターテーブルにアクセスするためには、"Load Postgis raster to QGIS"プラグインを使用します。インストールしておきましょう。

"レイヤ"メニューに、"Add a Postgis Raster Layer"が追加されます。

ここで一旦"レイヤ"→"PostGISレイヤの追加"に戻っていただいて、接続のための設定をしておく必要があります。

再度"Add a Postgis Raster Layer"を選択すると、ラスターテーブルを選べるようになります。

ここで、SRIDのところが"0"になってしまっているので、"4612"と入力して読み込んでみましょう。
バンド1を原色で表示した絵ですが、データベースに取り込まれていることを確認出来ます。

ラスタの投影変換

投影変換はst_transform()で行えます。
リサンプリング時のアルゴリズム等もオプションで指定出来ます。ここでは、JGD2000/平面直角座標系12系(SRID:2454)への変換として、SRIDだけ指定すると、

select st_transform(rast, 2454) from dem;

となります。
メタデータを覗いてみると、

投影変換されてそうですね。

ラスタのクリップ

ハンズオン資料(QGIS_RAS\dataset\shp\J_buf_z12)のポリゴンを使ってクリップを行ってみます。
まずはshapefileをデータベースへ取り込み。

shp2pgsql -c -I -s 2454 J_buf_z12.shp j_buf_z12 > j_buf_z12.sql
psql -U waigania12 -f j_buf_z12.sql rastertest

クリップはst_clip()にて行えます。
クリップした画像をdem_clipテーブルに保存するようにします。お互いのテーブルにデータは1行しかないので、特にwhere句も設定しません。

create table dem_clip(rid integer, rast raster);
insert into dem_clip(rid, rast)
select d.rid, st_clip(st_transform(d.rast, 2454), j.geom) 
from dem d, j_buf_z12 j;

QGISで結果を表示してみると、このような感じになります。

地形解析

st_hillshade()やst_slope()を用いて地形解析も行えます。
st_hillshade()を試してみましょう。

create table dem_hillshade(rid integer, rast raster);
insert into dem_hillshade(rid, rast) 
select rid, st_hillshade(rast, 1, '8BUI', radians(315), radians(45))
from dem;

QGISで結果を表示してみると、このような感じになります。

もうひとつおまけで、QGISのラスタ計算機のようなことをする場合、st_mapalgebraexpr()を使用出来ます。
標高1000m以上の地点を1、それ以外を0としたラスタを作るとすると、

create table dem_expr(rid integer rast raster);
insert into dem_expr(rid, rast) 
select rid, st_mapalgebraexpr(rast, 1, '1BB', 'CASE WHEN [rast] >= 1000 THEN 1 ELSE 0 END') 
from dem;

のように行えます。