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](http://www.gaia-gis.it/gaia-sins/)の下の方に、Windows向けバイナリのダウンロードへのリンクがあるので、"current stable version" をクリックしましょう。
![image.png](https://qiita-image-store.s3.amazonaws.com/0/122210/b3550e3a-778d-de99-3509-644a6fb0c42c.png)

ファイル一覧が出てくるので、 "spatialite_gui-4.3.0a-win-amd64.7z" をクリックしてダウンロードします。
![image.png](https://qiita-image-store.s3.amazonaws.com/0/122210/ba12db63-5680-2ddc-1f48-df92fcd31a6f.png)

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

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

![image.png](https://qiita-image-store.s3.amazonaws.com/0/122210/ce2d9171-402a-09ee-3773-d44629ef763d.png)

初回は新しいDBファイルを作成します。ツールバーの左から3つ目のアイコン "Creating a New (empty) SQLite DB" をクリックします。
![image.png](https://qiita-image-store.s3.amazonaws.com/0/122210/f73342a5-816d-4633-022f-8796aead1cdf.png)

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


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

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

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

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

```sql
-- 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カラムに地球アイコンがついていることを確認しましょう。
![image.png](https://qiita-image-store.s3.amazonaws.com/0/122210/98a9470f-dda6-7f5c-4060-6d76fc84b4b0.png)

# WKTでのジオメトリ表現

WKTはWell-known textの略で、地図上で幾何学オブジェクト（点とか直線とか多角形とか）を表現するためのマークアップ言語です（[Wikipediaの記事](https://ja.wikipedia.org/wiki/Well-known_text)）。
例えば、`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にある）を挿入します。

```sql
-- 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
	)
);
```

# 取得

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

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

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

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

```sql
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`関数を使います。

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

こうして得られたkml表現は、適当なファイルへコピーペーストして.kmlファイルとして保存すると、[Google Earth Pro](https://www.google.com/intl/ja/earth/desktop/)のようなソフトに読み込んで表示させることができます。