2013/10/17
PostGIS > PL/Python でラスタをファイル出力（3）
PostGIS > PL/Python でラスタをファイル出力（3）
1. 今日やったこと
一昨日・昨日の続き。PostGIS ラスタを SQL で PNG ファイルに出力する態勢が整ったので、R で PNG を読み込み、一枚に結合プロットした。R の処理を PL/R のストアド関数として定義した。結果、一連の処理が SQL 一文で済むようになった。実行環境は 2013/10/12 を、また PostgreSQL 9.3 への PL/R インストールは 2013/10/13 をそれぞれ参照。実行例。昨日、テスト用に下のようなラスタテーブルを準備した。ラスタの中身は、埼玉県内の市区町村ポリゴンに適当な色を付けたもの。
今日作った PL/R のストアド関数を下のように実行すれば、PDF で出力される。クエリでテーブル名を指定したり、必要行を選択する。関数本体は 2. で書く。
SELECT 自作関数名('出力ファイルパス', 'ラスタを呼び出すクエリ') ;下が出力結果。一行ごとに別々だったラスタが PNG ファイル出力を経由して一枚に結合され、経緯度の軸を付けてプロットされた。
https://kenpg.up.seesaa.net/image/20131016_test1.pdf
市区町村境界に白い線が描かれているように見えるが、これは Acrobat の表示による。各市区町村が別々のデータなので、自動的にデータの境界に白線が引かれる。下のように拡大して見ると、市区町村間に隙間はないが、表示の上で白線が出ることが分かる。
以下 PL/R のストアド関数の詳細を示す。
2. 自作した PL/R のストアド関数
定義文の全体と、中で呼び出すファイル出力関数の定義文は下記のとおり。後者は昨日の再掲。PL/R から PL/Python を呼び出している訳で、いわば R と Python 両方の機能を PostgreSQL のストアド関数上で連携させている。CREATE OR REPLACE FUNCTION 関数名(text, text)
RETURNS void LANGUAGE plr STABLE AS $BODY$
# PNG 残さないフラグ
del = TRUE
# PNG を PDF と同じ場所に出力するためのパス設定
pth = sub("[^/\\]+$", "", arg1)
sql = paste(sep="", '
WITH p AS (
SELECT ファイル出力関数名(
text \'', pth, '\' || row_number() over() || \'.png\'
, ST_AsPNG(rast)) fpath, (ST_MetaData(rast)).*
FROM (', arg2, ') foo
)
SELECT fpath
, upperleftx x1
, upperlefty + height * scaley y1
, upperleftx + width * scalex x2
, upperlefty y2
FROM p')
pg.thrownotice(sql)
graphics.off()
pdf(arg1)
df1 = pg.spi.exec(sql)
plot(range(c(df1[,"x1"], df1[,"x2"])), range(c(df1[,"y1"], df1[,"y2"]))
, type="n"
, xlab = "Longitude (degrees)"
, ylab = "Latitude (degrees)"
, asp = 1
)
for(i in 1:nrow(df1)){
png = df1[i, "fpath"]
pg.thrownotice(png)
rasterImage(png::readPNG(png, native=TRUE)
, df1[i, "x1"], df1[i, "y1"], df1[i, "x2"], df1[i, "y2"]
, interpolate = FALSE)
if(del) file.remove(png)
}
dev.off()
$BODY$ ;
CREATE FUNCTION ファイル出力関数(fpath text, bytes bytea)
RETURNS text LANGUAGE plpython3u STABLE AS
$BODY$
f = open(fpath, 'wb+')
f.write(bytes)
return fpath
$BODY$ ;
ラスタテーブルと関数に適当な名前を設定して実行した例を pgAdminⅢ のクエリツールで示す。メッセージ部分に PL/R から実行したクエリ内容が表示され、続いてそのクエリ内で出力された PNG ファイルパスが表示される。
3. ラスタの四隅座標を得るクエリ
このストアド関数のポイントが、上記メッセージ部分に示されたクエリ。ラスタを PNG ファイル出力すると同時に四隅座標を出力している。四隅座標は ST_MetaData 関数を使う。これは record 型を返し、その中に左端・上端の座標が入っている。また右端・下端の座標も計算して得られる。当該部分を再掲すると下のとおり。SELECT upperleftx x1 -- 左端
, upperlefty + height * scaley y1 -- 上端
, upperleftx + width * scalex x2 --右端
, upperlefty y2 -- 下端
FROM (
SELECT (ST_MetaData(ラスタ)).* FROM (ラスタの入っているテーブル)
) foo ;
4. R で PNG をプロットしたり貼り合わせる方法
今回の処理のひな型は次のようになる。plot(c(左端, 右端), c(下端, 上端), type="n" # 全体の範囲で空プロット
rasterImage(png::readPNG("PNGファイルその一", native = TRUE)
, 左端, 下端, 右端, 上端, interpolate = FALSE)
rasterImage(png::readPNG("PNGファイルその二", native = TRUE)
, 左端, 下端, 右端, 上端, interpolate = FALSE) ……
まず全体の範囲で空プロットする。次に PNG 一つずつ png ライブラリの readPNG 関数で読み込み、rasterImage で四隅座標を指定して描画する。貼り合わせといってもこれだけの素朴な方法。なお readPNG の native = TRUE 指定は、グレースケールの PNG を読み込むのに必須。また rasterImage の interpolate = FALSE は、画像を滑らかに見せる必要がないため指定した。
一つ一つの PNG は下のようになっていて、背景色を透過指定している。というか PostGIS ラスタの ST_AsPNG 関数が上手くできていて、NODATA 値のピクセルを自動的に透過色に指定してくれるのである（RGB 三色かグレースケールの場合。アルファチャンネルがある場合は別）。そして R の readPNG & rasterImage の組み合わせも上手くできていて、PNG の透過指定をきちんと反映してくれる。その結果、四隅座標さえ適切に指定すれば、市区町村などの「貼り合わせ」ができる訳。
5. 複数テーブルのラスタも一枚に出力可能
描画したいラスタが複数のテーブルに分かれていても、今回の関数に渡すクエリで UNION ALL すれば一枚にプロットできる。クエリの結果順に描画されるので、もし重なる部分があれば UNION ALL の後ろほど上に来る。SELECT 自作関数名('出力ファイルパス', '
SELECT rast FROM テーブルその一
UNION ALL
SELECT rast FROM テーブルその二 -- こちらが上になる
') ;今回の関数は、描画したいラスタの解像度や色がバラバラでも問題なく一枚にプロットできる。また PL/R で追加的なプロットをすることも可能。それら一連の処理が SQL によるストアド関数実行だけで済むようになった。
2. の繰り返しになるが、今回の方法は PL/R から PL/Python を呼び出すことにより、R と Python 両方の機能を PostgreSQL のストアド関数上で連携させている。自分にとって PostgreSQL は、データベースだけでなく複数のスクリプト言語の統合プラットフォームみたいになってきた。