地図とかの空間情報をSQLiteに格納するSpatiaLiteを使用してみる

##概要
SpatiaLiteは地図などの空間情報を格納できるSQLiteの拡張です。
ShpファイルやGeoJsonで記述されたGeometry情報をデータベースに格納し利用できます。

##配布ファイルの説明
###SpatiaLite
コマンドラインからデータベースの操作を行えます。

バイナリの配布
Windowsの場合は以下から任意のプラットフォームのバイナリが取得できます。
http://www.gaia-gis.it/gaia-sins/index.html

###mod_spatialite.dll/mod_spatialite.so
SQLiteの拡張モジュールです。
SQLite上で、load_extensionを利用してDLLまたはsoファイルを読み込むことでSpatiaLiteの機能が使用できるようになります。

バイナリの配布
Windowsの場合は以下から任意のプラットフォームのバイナリが取得できます。
http://www.gaia-gis.it/gaia-sins/index.html

ソースコードの配布
以下からlibspatialiteのソースコードをダウンロードしてビルドすることで作成できます。
https://www.gaia-gis.it/fossil/libspatialite/index

###spatialite_gui
GUIでspatialiteの機能が使用できます。

POLYGONなどのGeometry型の列については、画像としてその内容を確認できます。

バイナリの配布
Windowsの場合は以下から任意のプラットフォームのバイナリが取得できます。
http://www.gaia-gis.it/gaia-sins/index.html

##libspatialiteのビルド
もし、バイナリでmod_spatialiteを取得できない場合は、libspatialiteをビルドする必要があります。

libspatialite-4.2.0のビルドには下記のライブラリが必要です。
proj-4.8.0
proj.4は地図作成投影ライブラリです。
http://trac.osgeo.org/proj/

geos-3.4.2
Geometry EnginはJavaTopologySuiteのC++への移植で、Geometryを処理するための空間述語と機能のAPIです。
http://trac.osgeo.org/geos/

以上をインストールしたあとに、libsatialiteのソースがあるフォルダにて下記のコマンドを実行してください。

./configure --disable-freexl
make
make install

--disable-freexlでExcelとの連携機能を無効にしてます。
これを有効にした場合はfreexlが必要になります。
https://www.gaia-gis.it/fossil/freexl/index

###FreeBSDでエラーが出る場合
FreeBSDで下記のようなエラーがでる場合があります

/usr/bin/ld: cannot find -ldl

FreeBSDではdlopen は標準ライブラリなので-ldlは不要です。
src/Makefileで-ldlを検索して削除したのちにmakeをし直しましょう。

参考：
http://sourceforge.net/p/idjc/discussion/458834/thread/246f0841/

##mod_spatialiteをSQLiteから使用する。
先に述べたように、load_extensionでspatialiteが利用できますが、いくつか注意点があります。

・SQLiteのバージョンは3.7.17以降である必要があります。

・依存するDLLはすべて、SQLiteから認識できる場所に配置される必要があります。関連のDLLのあるフォルダを環境変数PATHで指定してください。

・SQLiteとmod_spatialiteのプロセスは同じビットでビルドしなければいけません。
　SQLiteが32bitで動作している場合はmod_spatialiteは32bitでビルドする必要がり、SQLiteが64bitで動作している場合はmod_spatialiteは64bitでビルドする必要があります。
　もし、Windowsで64bitバイナリのSQLiteが必要な場合は下記を記事通りに行えば作成できます。
http://qiita.com/akaneko3/items/0e99c3c1366dfbad006f

・load_extensionを利用してDLLまたはsoファイルを読み込むことで、以降spatialiteの機能が利用できるようになります。

SELECT load_extension('/usr/local/lib/mod_spatialite.so');

・SQLite3からデータベースを作成した場合、SpatiaLiteが管理するメタデータが作成されていません。このメタデータを作成するために下記のSQLを実行します。

SELECT InitSpatialMetaData();

##チュートリアル
ここでは、実際に国土数値情報の行政区域をspatialiteに格納して利用してみます。

###国土数値情報で全国の行政区域を取得する。
まず下記のページから全国の行政区域をダウンロードします。
http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03.html

###空のSQLiteファイルを作成する。
spatialite_guiを起動して空のSQLiteのデータベースファイルを作成します。

[Creating a New (empty) SQLite DB]アイコンを押して、ファイル名を入力します。

正常に作成されると、作成されたデータベースへのパスが表示されます。

ここまでは普通のデータベースと変わりありません。

###shpファイルをインポート
国土数値情報のshpファイルをデータベースに取り込みます。

[Load Shape file]アイコンを押して、国土数値情報の「N03-14_140401.shp」を選択します。

この時、エンコードを聞かれるので「CP932」を選択してください。

正常に完了すると「N03-14_140401」というテーブルが作成されます。

###shpファイルから作成されたテーブルの確認
ここではshpファイルから作成されたテーブルの確認を行います。
作成されたテーブルの列を確認するにはテーブルを選択して、右クリックを押して「Show columns」を選択します。

N03_001～N03_004は県、区、市、N03_007にはコードが格納されているTEXT型の列です。
Geometryは行政区域の形がBlobとして格納されていてPOLYGONというspatiaLiteが指定した型になっています。

このテーブルの作成に実際しようされたSQL文はテーブルを右クリックして[Show CREATE statement]で確認できます。

CREATE TABLE "N03-14_140401" (
"PK_UID" INTEGER PRIMARY KEY AUTOINCREMENT,
"N03_001" TEXT,
"N03_002" TEXT,
"N03_003" TEXT,
"N03_004" TEXT,
"N03_007" TEXT, 
"Geometry" POLYGON)

では実際にテーブルの内容を見てみます。
そのために、SQL文を実行します。

SELECT * FROM "N03-14_140401" LIMIT 10;

実行すると、以下のような情報が取得できます。

PK_UID～N03_007はテキストでデータが表現されています。
しかし、GeometryはBLOBデータなので人が見て内容が確認できません。

この内容を確認する方法は２通りあります。
１つはBLOB exploreを使用する方法です。
Geometryの列の任意のセルを選択して右クリックを押し「BLOB explore」を選択します。

BLOB exploreではバイナリ、画像、SVG、GeoJSONなどの様々な形式でデータを表現することができます。

BLOB exploreを使用しない方法としては、SQLでGeometryの列を任意の形式で表示する方法があります。

以下の例ではテキスト形式でGeometry列を表示しています。

-- データが長いので最初の100文字まで表示
SELECT substr(ASTEXT(Geometry), 1, 100) FROM "N03-14_140401" LIMIT 10;

ASTEXTでPOLYGON情報を文字として表示できました。これは任意の形式で出力することが可能です。たとえば、ASGEOJSONとした場合は、GeoJSONとして出力されます。

SELECT substr(ASGEOJSON(Geometry), 1, 100) FROM "N03-14_140401" LIMIT 10;

ASTEXTやASGEOJSONのSpatiaLiteの提供する関数は下記の「Spatial SQL functions reference guide」に記載されています。

###自分で座標情報を持つテーブルを作成してみる
先の例ではshpファイルからテーブルを作りましたが、今回は、座標情報を持つテーブルを自分で作ってみます。

以下のテーブルでは東京タワーなどの地物情報の座標を持つテーブルを作成します。

まず、地物情報を含まない列のみを作成します。

CREATE TABLE "places" (
"PK_UID" INTEGER PRIMARY KEY AUTOINCREMENT,
"name" TEXT)

SQLを実行した後に、refreshを行うとテーブルが追加されます。

次に、地物情報を表す列をAddGeometryColumn()を用いて作成します。

Select AddGeometryColumn ('places', 'Geometry', 0, 'POINT', 'XY')

SQLを実行した後に、refreshを行うと列が追加されて、テーブルに地球のアイコンがつきます。

テーブルを作成できたら、次はデータを入力します。
以下の例では都庁の座標をテーブルに格納しています。

INSERT INTO places (name, Geometry) VALUES(
  '都庁',
  GeomFromText('POINT(139.692101 35.689634 )')
)

GeomFromTextは文字列をGeometryに変更する関数です。今回はテキスト形式から変換しましたが、GeomFromGeoJSON、GeomFromKmlなどを利用してGeoJSONやKmlから作成することができます。

なお、以下のように一気に作ると、一見それっぽく動作しますが、メタデータが関連付けられていないので、RTreeインデックスを作るときなどに正常に動作しません。

CREATE TABLE "places" (
"PK_UID" INTEGER PRIMARY KEY AUTOINCREMENT,
"name" TEXT,
"Geometry" POINT)

###座標情報に対する検索
ここでは座標情報について検索します。
placesテーブルで指定した座標が、どの行政区画に含まれているか検索してみます。

SELECT
  name, "N03-14_140401".*
FROM
 "N03-14_140401"
INNER JOIN places ON
 Contains("N03-14_140401".Geometry, places.Geometry)

Contains関数は第一引数の範囲に第二引数の領域が含まれているか調べます。
この結果は次のようになります。

都庁は新宿区にあるので期待どおりの結果といえるでしょう。
先の例ではテーブルに格納されているデータで検索しましたが、直接WHERE区に指定することも可能です。

SELECT
  *
FROM
 "N03-14_140401"
WHERE
  Contains("N03-14_140401".Geometry,GeomFromText('POINT(139.692101 35.689634 )'))

###MBRと大まかな検索
先ほどの検索は各座標について正確に検索していましたが、大まかな位置を早く取得したい場合もあります。

この場合は、Minimum Bounding Rectangle - MBRを使用します。
MBRはGEOMETRYのおおまかな大きさを表すものです。

以下のような複雑な図形を赤い枠線の図形とみなします。

これにより、不正確ではあるが、早い検索が行えます。
GEOMETRYからMBRを取得するには Envelope()を用います。

SELECT
   ASTEXT(Envelope(Geometry))
FROM
 "N03-14_140401"
LIMIT 10

この結果は以下のように単純な図形になります。

では、先ほどと同じ都庁を含む座標をMBRを利用して検索してみましょう。

SELECT
  *
FROM
 "N03-14_140401"
WHERE
  MBRContains("N03-14_140401".Geometry,GeomFromText('POINT(139.692101 35.689634 )'))

検索速度が向上したかわりに以下のように、不要なデータも抽出されてしまっています。

###RTreeインデックスを利用した検索
通常、DBを扱う場合、インデックスを利用することで検索速度が向上します。
残念なことに、文字や数値に付与できるインデックスをGeometryに付与することはできません。
しかしながら、SpatiaLiteはSQLiteのRTreeを利用して、RTreeインデックスをサポートしています。
RTreeインデックスのアルゴリズムについては下記を参考にしてください。

Wonderful R*Tree Spatial Index
http://www.gaia-gis.it/gaia-sins/spatialite-cookbook/html/rtree.html

RTreeの概要としては、GeometryをMBRをつかって四角形(rectangle)にし、その四角形の交わりに応じてツリーを作り、検索しやすくしています。
rectangleを使ったツリーを作成しているのでRTree indexと言われます。

####RTreeインデックスの作成と利用
では、RTreeインデックスを実際使用してみます。
RTreeインデックスを作成するにはCreateSpatialIndex関数を使用します。第一引数にテーブル名、第二引数に列名を指定します。

SELECT CreateSpatialIndex("N03-14_140401", "Geometry") ;

正常に作成できた場合は、1が返り、refresh後にSpatialIndexに新しいテーブルが表示されることが確認できます。

RTreeインデックスを作成することにより、以下の４つのテーブルが作成されます。

・idx_N03-14_140401_Geometry
・idx_N03-14_140401_Geometry_node
・idx_N03-14_140401_Geometry_parent
・idx_N03-14_140401_Geometry_rowid

idx_N03-14_140401_Geometry以外は内部で使用するテーブルになっており、ユーザーが直接操作することを期待していません。

idx_N03-14_140401_GeometryはVIRTUAL TABLEになっており、実際には内部テーブルを利用して結果を返します。

では、idx_N03-14_140401_Geometryの内容を確認してみます。

SELECT * FROM "idx_N03-14_140401_Geometry" LIMIT 10;

このように各レコードのMRBが格納されています。

ユーザはこのテーブルを利用することで、大量のデータの中からデータをフィルターをし、結果を求めることができます。

SELECT 
  * 
FROM "N03-14_140401"
WHERE ROWID IN(
  SELECT
    pkid
  FROM
   "idx_N03-14_140401_Geometry"
  WHERE
    MBRContains(
      BuildMBR(xmin,ymin,xmax,ymax),
      GeomFromText('POINT(139.692101 35.689634 )')
  )
) AND Contains(Geometry, 
  GeomFromText('POINT(139.692101 35.689634 )')
)

通常の検索より早く正確なデータが取得できたかと思います。これはデータ数が多くなれば、なるほどその差が大きく現れます。

このようにRTreeインデックスは通常のインデックスと違い、実テーブルにSQLを実行しただけでは適用されずに、作成されたVIRTUAL TABLEを利用してサブクエリ―やJOINでデータをフィルタリングする役割になります。

別テーブルでインデックスを管理していることで、実テーブルが更新された時に不一致がでるかもしれないという疑問をもつでしょう。
しかし、CreateSpatialIndexを実行した時にトリガーが作成されているので、ユーザは一切手を加えずにインデックスとテーブルの同期がとれます。

このトリガーを確認するには次のようなSQLを実行します。

select name  from sqlite_master where type = 'trigger' and tbl_name='N03-14_140401';

今回は次のようなトリガーが自動で作成されていることが確認できます。

・ggi_N03-14_140401_Geometry
・ggu_N03-14_140401_Geometry
・gii_N03-14_140401_Geometry
・giu_N03-14_140401_Geometry
・gid_N03-14_140401_Geometry

####RTreeインデックスの削除
作成したRTreeインデックスを削除するには、DisableSpatialIndexを使用します。

SELECT DisableSpatialIndex("N03-14_140401", "Geometry") ;

refreshを実行するとSpatial Indexが削除されていることが確認できます。

##pythonからの利用
pythonからSpatiaLiteを使用する場合、SQLiteが適切に利用できる状況なら簡単に利用できます。

Pythonが使っているSqlite3のバージョンは次の手順で調べることができます。

import sqlite3
print (sqlite3.sqlite_version_info)

もし、バージョンを満たさない場合は、SQLITEを更新してください。
Windowsの場合は以下のフォルダのDLLを書き換えるといいでしょう。

C:\Python27\DLLs\sqlite3.dll

以下はPythonでのコードになります。

#!/usr/bin/python
# -*- coding: utf-8 -*-
import sqlite3
import os

# mod_spatialiteのあるフォルダをPATHに加える
os.environ["PATH"] = os.environ["PATH"] + ';C:\\tool\\spatialite\\mod_spatialite-4.2.0-win-x86'
cnn = sqlite3.connect('database/gyouseikuiki.sqlite')

# mod_spatialiteの読み込み
cnn.enable_load_extension(True)
cnn.execute("SELECT load_extension('./mod_spatialite-4.2.0-win-x86/mod_spatialite.dll');")
sql = """
SELECT
  N03_001,
  N03_002,
  N03_003,
  N03_004,
  N03_007, 
  AsGeoJson(Geometry)
FROM
 "N03-14_140401"
WHERE
  MBRContains("N03-14_140401".Geometry,GeomFromText('POINT(139.692101 35.689634 )'))
"""

ret = cnn.execute(sql)
for r in ret:
  print('----------------------------')
  print(r[0].encode('utf_8') )
  print(r[1].encode('utf_8') )
  print(r[2].encode('utf_8') )
  print(r[3].encode('utf_8') )
  print(r[4].encode('utf_8') )
  print(r[5].encode('utf_8') )

環境変数をmod_spatialiteのあるパスに通し、enable_load_extensionで拡張DLLの使用を許可してから、mod_spatialiteをロードします。

あとは、通常通りSQLが実行できます。

###enable_load_extensionでエラーになる場合
enable_load_extensionを使用した時に以下のようなエラーが発生する場合があります。

"'sqlite3.Connection' object has no attribute 'enable_load_extension'"

原因として、MaxOSやDebianなどのいくつかのOSにプリインストールされているPythonはコンパイル時に、拡張ライブラリの読み込み機能を無効にしているためです。

これに対応するにはPythonをインストールしなおすしかありません。
インストールする際、Pythonのソースコードのフォルダにある、setup.pyを手で修正する必要があります。

SQLITE_OMIT_LOAD_EXTENSIONを検索して、その行をコメントアウトして、./configure, make, make install を実行します。

もし、make install でエラーが発生する場合は、make install時にエラーがでるなら「make -i altinstall」をためしてみてください。

###Peeweeを利用する方法
ORMであるPeeweeライブラリでspatialiteを使用する方法については下記を参考にしてください。

##phpからの利用
Windows+PHPの場合、簡単に使うことができません。
詳細は下記を参照してください。

PHPで地図とかの空間情報をDBに格納しようとしてSpatialiteを使うと苦労する
http://qiita.com/mima_ita/items/1a90de89f0a194ba843e

##その他Tip
###測地系の変換をSpatialiteを使って行う
Transform関数を使うと測地系の変換が行えます。
以下の例ではJGD2000/平面直角座標系6を世界測地系に変換しています。

select AsText(Transform(GeomFromText('POINT(-4408.916645 -108767.765479)', 2448), 4326))

##参考
SpatiaLite 4.2.0 SQL functions reference list
http://www.gaia-gis.it/gaia-sins/spatialite-sql-4.2.0.html

The SpatiaLite Cookbook
http://www.gaia-gis.it/gaia-sins/spatialite-cookbook/index.html

A quick tutorial to SpatiaLite - a Spatial extension for SQLite(古いが有用)
http://www.gaia-gis.it/gaia-sins/spatialite-tutorial-2.3.1.html