2014/09/02
PostGIS
 >  SVG のパスを取り込む(1)
openSUSE に PostGIS ラスタを入れられなかった件を 2014/08/2408/25 に書いた。当面 geeko の PNG インポートはできないが、元は SVG - Scalable Vector Graphics なのでパスを解析してベクタとして取り込めるかもしれない。その試みの一回目。まだ直線のみの未完成段階なので、実行環境は 2014/08/13 の実機。完成したら当然 openSUSE の PostGIS に入れる。

http://en.wikipedia.org/wiki/File:OpenSUSE_official-logo-color.svg


↓ 今日できた所を R でプロットした。そのコードは最後に。先に SVG ファイルのダウンロードからの一連の手順を書く。


最初に SVG の中身を丸ごと PostgreSQL に入れる。SVG は XML の一種でテキストファイルなので COPY 文を使える。Wikipedia のロゴページ(冒頭の URL)から SVG をダウンロードして PostgreSQL ユーザが読み取れる場所に置き、適当な一時テーブル名を決めて ↓ のように一行ずつインポート。SVG 内にタブがあるので、まず存在しない制御文字(バックスペース)をデリミタに指定することで全行一列で認識させた。
DO $D$
BEGIN
EXECUTE format('
CREATE TABLE "201409"."02_svg_tmp"
(rid serial PRIMARY KEY, str text) ;
COPY "201409"."02_svg_tmp" (str) FROM %L (DELIMITER %L)'
, 'R:/OpenSUSE_official-logo-color.svg', E'\b') ;
END $D$ ;


上の rid 列は行番号。この順序で全行を連結して一つながりの文字列にし、正規表現で path 要素を取り出して別のテーブルに保存した。↓ path の d 属性が描画コマンドで、線種と座標が列挙されている。fill は塗り色で当面使わないが、後々のために入れた。
CREATE TABLE "201409"."02_svg_path" AS
WITH a AS (
SELECT * FROM "201409"."02_svg_tmp" ORDER BY rid
), b AS (
SELECT string_agg(str, ' ') str
FROM a
), c AS (
SELECT regexp_matches(str
, '<path fill="([^"]+)" d="([^"]+)"', 'g') reg
FROM b
)
SELECT row_number() over() :: int pid
, reg[1] fill
, regexp_replace(reg[2], '\s+', '', 'g') d
FROM c ;


上の pid 列が各 path 要素の ID。この SVG には14個の path があり、それぞれの d 属性を解釈して PostGIS のポリゴンにする。と言ってもベジェ曲線のように直接ジオメトリに変換できない線種もあるので模索が必要。とりあえず各 path を一線種ずつに分割し、座標を配列に入れて返す PL/pgSQL のストアド関数を作った。d 属性の詳細は下記二サイトが大変参考になった。

■ 京都大学 須崎純一先生 : SVGの基本演習
http://www.envinfo.uee.kyoto-u.ac.jp/user/susaki/envinfo/svg_basic.html

■ Mozilla Developer Network : SVG チュートリアル Paths
https://developer.mozilla.org/ja/docs/Web/SVG/Tutorial/Paths

↓ 作成したストアド関数の定義文。苦労したのは、今回の SVG が Adobe Illustrator によるもので d 属性の記述の仕方がやや特殊なのと、相対座標の処理が面倒だったこと。
CREATE OR REPLACE FUNCTION "201409"."02_svg_pathcmd"(text)
RETURNS TABLE (sid int, cid int, str text, cmd text, pxy numeric[])
LANGUAGE plpgsql IMMUTABLE AS $BODY$
DECLARE r record ;
cmd text ;
ar1 numeric[] ;
ar2 numeric[] ;
pxy numeric[] ;
BEGIN
FOR r IN
WITH b AS (
SELECT (regexp_matches($1, '(\w[\d\-\.\,]+)', 'g'))[1] str
), c AS (
SELECT row_number() over() :: int cid, b.str, left(b.str, 1) lab
FROM b
)
SELECT sum(CASE WHEN lab ~* 'M' THEN 1 END)
over(ORDER BY c.cid) :: int sid, *
FROM c ORDER BY 1, c.cid
LOOP
cmd = r.lab ;
ar1 = string_to_array(regexp_replace(right(r.str, -1)
, '(\d)-', '\1,-', 'g'), ',') ;
IF cmd ~* 'M' THEN
pxy = CASE WHEN cmd = 'M' THEN ar1
ELSE ARRAY[pxy[1] + ar1[1], pxy[2] + ar1[2]]
END ;
ar2 = NULL ;
ELSEIF cmd ~* '[HVLT]' THEN
ar2 = pxy || CASE
WHEN cmd = 'H' THEN ARRAY[ar1[1], pxy[2]]
WHEN cmd = 'h' THEN ARRAY[pxy[1] + ar1[1], pxy[2]]
WHEN cmd = 'V' THEN ARRAY[pxy[1], ar1[1]]
WHEN cmd = 'v' THEN ARRAY[pxy[1], pxy[2] + ar1[1]]
WHEN cmd ~ '[LT]' THEN ar1
ELSE ARRAY[pxy[1] + ar1[1], pxy[2] + ar1[2]]
END ;
pxy = ar2[3:4] ;
ELSEIF cmd ~* '[QS]' THEN
ar2 = pxy || CASE
WHEN cmd ~ '[QS]' THEN ar1
ELSE ARRAY[pxy[1] + ar1[1]
, pxy[2] + ar1[2]
, pxy[1] + ar1[3]
, pxy[2] + ar1[4]]
END ;
pxy = ar2[5:6] ;
ELSEIF cmd ~* 'C' THEN
ar2 = pxy || CASE
WHEN cmd = 'C' THEN ar1
ELSE ARRAY[pxy[1] + ar1[1]
, pxy[2] + ar1[2]
, pxy[1] + ar1[3]
, pxy[2] + ar1[4]
, pxy[1] + ar1[5]
, pxy[2] + ar1[6]]
END ;
pxy = ar2[7:8] ;
END IF ;
IF ar2 IS NOT NULL THEN
RETURN QUERY SELECT r.sid, r.cid, r.str, cmd, ar2 ;
END IF ;
END LOOP ;
END ;
$BODY$ ;

-- Example
WITH a AS (
SELECT *, "201409"."02_svg_pathcmd"(d) rec
FROM "201409"."02_svg_path"
)
SELECT pid, fill, (rec).* FROM a ;


↑ ストアド関数の使用例。列 pid が path 要素の ID、列 sid が各 path に含まれるポリゴン(2 以下は中マドの役割になる)。列 cid が一つずつの線で、cmd がその種類(大文字は絶対座標、小文字は相対座標)。pxy が絶対座標に変換したもの。例えば cmd = C なら三次ベジェ曲線で pxy は 8 個の要素を持つ(起点のXY、その制御点のXY、終点の制御点のXY、終点のXY)。

ベジェ曲線の処理は明日以降とし、とりあえず制御点を除く全ての点と、各線の起終点を真っすぐ結んだ線を PostGIS で出力してみた。↓ ST_Point の第二引数(経度)に -1 を掛けているのは、SVG と PostGIS で縦軸の方向が逆なため。
WITH a AS (
SELECT pid, "201409"."02_svg_pathcmd"(d) rec
FROM "201409"."02_svg_path"
), b AS (
SELECT pid, (rec).* FROM a
), c AS (
SELECT pid, pxy, array_length(pxy, 1) len FROM b
), d AS (
SELECT pid, ST_Point(pxy[1], pxy[2] * -1) p1
, ST_Point(pxy[len - 1], pxy[len] * -1) p2
FROM c
)
SELECT pid, p1, ST_MakeLine(p1, p2)
FROM d ;


上の結果を R で読み込みプロットしたのが、冒頭の直線だけの geeko。三つのパッケージ(RPostgreSQL、rgeos、maptools)を使い、PostGIS の点と線を ST_Collect で一行に束ね、ST_AsText で WKT 形式にして読み込み、rgeos の readWKT 関数で R のジオメトリオブジェクトに変換して maptools でプロット。↓
library(RPostgreSQL)
library(maptools)
library(rgeos)
con = dbConnect(PostgreSQL(), dbname='xxx', port=xxx, user='xxx', password='xxx')
df1 = dbGetQuery(con, '
WITH a AS (
SELECT "201409"."02_svg_pathcmd"(d) rec
FROM "201409"."02_svg_path"
), b AS (
SELECT (rec).* FROM a
), c AS (
SELECT pxy, array_length(pxy, 1) len FROM b
), d AS (
SELECT ST_Point(pxy[1], pxy[2] * -1) p1
, ST_Point(pxy[len - 1], pxy[len] * -1) p2
FROM c
)
SELECT ST_AsText(ST_Collect(p1))
, ST_AsText(ST_Collect(ST_MakeLine(p1, p2)))
FROM d ;
')
invisible(dbDisconnect(con))
poi = readWKT(df1[1,1])
lin = readWKT(df1[1,2])
windows(height=(340-46)/96, width=(444-8)/96)
par(plt=rep(c(1,99)/100, 2), xaxs='i', yaxs='i')
plot(poi, cex=0.3, pch=16, col=3)
plot(add=T, lin, col=3)


↓ プロット結果の再掲。次回はベジェ曲線を処理して、何とかオリジナルの形に近づけたい。それができたら、各 path をきちんと閉じてポリゴンにしたり、中マドのあるポリゴンを PostGIS のマルチポリゴンに置き換える予定。そうしてやっと、R で可視化する際に色を塗れる。
<< 3次ベジェ曲線の近似
順列生成の PL/pgSQL と再帰クエリ >>