SpatiaLiteでジオメトリのカラムを作り、データを入れ、kmlで出力する

SpatiaLiteでジオメトリ型のカラムを作ったり、そこへデータを入れたり取得したりする方法をすぐ忘れるので、備忘録的に書いておきます。
環境として、Windows 10（64 bit）を前提として説明します。

SpatiaLiteとは何か

SQLiteで地理空間情報を扱えるようにしたものがSpatiaLiteです。また、その拡張モジュールをmod_spatialiteと呼びます。
PostgreSQLにおけるPostGISと同じような位置づけです。

spatialite_guiの利用

Windows 10ではGUIで利用できるspatialite_guiがあるので、それを利用します。

ダウンロード

The Gaia-SINS federated projects home-pageの下の方に、Windows向けバイナリのダウンロードへのリンクがあるので、"current stable version" をクリックしましょう。

ファイル一覧が出てくるので、 "spatialite_gui-4.3.0a-win-amd64.7z" をクリックしてダウンロードします。

最新の安定バージョンの公開日が2015年9月7日...最近は活発な更新が行われていないようですね。

ダウンロードした7zファイルは、7-Zipで解凍できます。
解凍すると "spatialite_gui.exe" という1ファイルだけが出てくるので、これをダブルクリックすれば起動します。

初回は新しいDBファイルを作成します。ツールバーの左から3つ目のアイコン "Creating a New (empty) SQLite DB" をクリックします。

表示されたファイルダイアログで、適当な名前をつけてファイルを保存しましょう。
2度目からは、一番左のアイコン "Connecting an exsiting SQLite DB" で、作成済みのファイルを開けばOKです。

ジオメトリ型カラムの作成

SpatiaLiteで気をつけなければらないのは、ジオメトリ型のカラムの扱いです。
例えば、3次元のPOLYGONを格納するカラムを持つテーブルを作成したいとしましょう。PostGISでは、以下のようにCREATE TABLE文一つで解決です。

-- PostGISではOK。SpatiaLiteではNG
CREATE TABLE polygons (
    name TEXT,
    geom GEOMETRY(POLYGONZ, 4326)
);

（ここで4326はSRIDで、Google MapやGoogle Earthで使われているWGS84の座標系を示す番号です）
SpatiaLiteでは、上記のような書き方はできません。
ジオメトリ型のカラムはCREATE TABLEをしたあとに、AddGeometryColumnという関数で追加する必要があります。

-- SpatiaLiteではこの書き方をする。
CREATE TABLE polygons (
    name TEXT
);
SELECT AddGeometryColumn(
    'polygons',         -- カラムを追加するテーブルの名前
    'geom',             -- 追加するカラムの名前
    4326,               -- SRID
    'POLYGONZ',         -- ジオメトリの型を指定
    'XYZ',              -- ジオメトリの次元を指定
    1                   -- 1にするとNOT NULL制約をつける。0だとつけない
);

上記SQLを画面右上のSQL入力欄にコピペして、右端の "Execute SQL Statement" ボタンをクリックすると、SQLが実行されます。
画面右下の結果セット表示欄に "1" が表示されれば正常完了です。
画面左側のツリーの "User Data" を右クリックして "refresh" すると、polygonsテーブルが追加されています。ノードを開いて、geomカラムに地球アイコンがついていることを確認しましょう。

WKTでのジオメトリ表現

WKTはWell-known textの略で、地図上で幾何学オブジェクト（点とか直線とか多角形とか）を表現するためのマークアップ言語です（Wikipediaの記事）。
例えば、x = 0、y = 1の二次元の点は以下のように表せます。

POINT (0 1)

3次元だと数字の並びが一つ増えます。幾何学の型にも "Z" が付きます。
(x, y, z) = (1, 2, 3)の点は以下のように表します。

POINT Z (1 2 3)

多角形はPOLYGONです。2次元空間において、左下隅が原点にくる辺の長さ1の正方形は以下のように表せます。

POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))

(x, y) = (0, 0)を始点として、反時計回りに頂点座標を順に記述します。
POLYGONの注意点として、始点と終点（今回は原点(0 0)）を一致させる必要があります。
最後の(0 0)を記述しないと、 "unclosed" とみなされます。

インサート

ジオメトリ型カラムに値を挿入するときは、WKT表現のテキストをバイナリのジオメトリデータに変換する関数GeomFromTextを利用します。

GeomFromText(<WKT文字列>, <SRID>)

例えば、以下のSQLはGooGle Earthのkmlの例にあるように、ペンタゴンを囲む五角形のポリゴン（高さ30mにある）を挿入します。

-- pentagon!
INSERT INTO polygons(name, geom)
VALUES
(
    'pentagon',
    GeomFromText(
        'POLYGON Z ((' ||
        '-77.0578845766096000 38.8725325989282000 30.0000000000000000, ' ||
        '-77.0546597375670000 38.8729101628170000 30.0000000000000000, ' ||
        '-77.0531553685479000 38.8705326779438000 30.0000000000000000, ' ||
        '-77.0555262249351000 38.8687578012560000 30.0000000000000000, ' ||
        '-77.0584405629039000 38.8699620650694000 30.0000000000000000, ' ||
        '-77.0578845766096000 38.8725325989282000 30.0000000000000000))',       -- WKT
        4326                                                                    -- SRID
    )
);

取得

ジオメトリ型のカラムをそのまま指定すると、バイナリデータのまま出力されるので人間が読めません。

-- geomカラムは人間に読めない
SELECT name, geom FROM polygons;

このとき、spatialite_guiではgeomカラムの値は "BLOB sz=196 GEOMETRY" みたいな表記になります。

そこで、人間がとりあえずこの値を見たいときは、このバイナリデータをWKTに変換する関数AsTextを使います。

SELECT AsText(geom) FROM polygons;
-- 出力
-- POLYGON Z((-77.057885 38.872533 30, -77.05466 38.87291 30, -77.053155 38.870533 30, -77.055526 38.868758 30, -77.058441 38.869962 30, -77.057885 38.872533 30))

KML出力

テーブルに格納したジオメトリをkmlで出力するには、AsKml関数を使います。

SELECT
    AsKml(
        name,       -- name
        name,       -- description
        geom,       -- geometry
        16      -- precision
    ) AS kml
FROM test_geom;

こうして得られたkml表現は、適当なファイルへコピーペーストして.kmlファイルとして保存すると、Google Earth Proのようなソフトに読み込んで表示させることができます。