地図の塗り分けに関する話で「4色問題(定理)」というのがあります。簡単に説明すると「どんな地図でも4色あれば同じ色が隣接しないように塗り分けられる」というものです。今回はそれをSpatialiteを使って実装する方法を紹介.....のつもりでしたが途中で挫折しました。またその気になったときに思い出せるように途中経過をメモっておきます。
ちなみに4色問題自体は解決(?)しています。あと、これは証明方法の実装ではありません。
http://ja.wikipedia.org/wiki/%E5%9B%9B%E8%89%B2%E5%AE%9A%E7%90%86
プログラムの準備
- OSGeo4Wでgdalをインストール
- spatialiteをとってくる
http://www.gaia-gis.it/gaia-sins/windows-bin-x86/spatialite-4.1.1-win-x86.zip
データの準備
- ポリゴンデータをspatialiteの形式で用意します。
- shpからでも以下のコマンドで変換できます。
ogr2ogr -gt 65536 -append -nlt polygon -s_srs epsg:4612 -t_srs epsg:4612 -f SQLite veg.sqlite veg.shp -dsco SPATIALITE=YES
- カラムに整数型のcolorを追加しておきます。（color列には、この後のプログラムで1〜4の塗り分けのための数字が入ります）
spatialite.exeを使って以下のようにします。データベースはveg.sqlite、テーブルはp574036という名前の場合です。
spatialite> alter table p574036 add color integer; spatialite> update p574036 set color=0;
プログラム
こちらを参考にpython＆spatialiteに移植しました。キモは隣接するポリゴンの抽出を、その都度spatialiteのST_TOUCHESを利用している部分です。
http://itpro.nikkeibp.co.jp/article/COLUMN/20070312/264600/
# coding:utf-8
from pyspatialite import dbapi2 as db
def isusedcolor(area,color):
query = 'select a.color from p574036 a, p574036 b where b.OGC_FID='
query += '%d ' % area
query += 'and ST_TOUCHES(a.geometry,b.geometry)'
#print query
cur.execute(query)
m = cur.fetchall()
for checkareacolor in m:
if color==checkareacolor[0]:
return 1
return 0
def drawcolor(area,color):
query = 'update p574036 set color='
query += '%d ' % color
query += 'where OGC_FID='
query += '%d' % area
cur.execute(query)
con.commit()
def getcolor(area):
query = 'select color from p574036 '
query += 'where OGC_FID='
query += '%d' % area
cur.execute(query)
color = cur.fetchone()
return color[0]
def colorsearch():
area = 1
color = 1
while True:
if color >= 5:
area = area -1
color = getcolor(area)
drawcolor(area,0)
color = color + 1
continue
if isusedcolor(area,color)==1:
color = color+1
else:
#print area,
#print color
drawcolor(area,color)
area = area + 1
color = 1
if area > AREA_LIMIT:
break
con = db.connect("veg_p57.sqlite")
cur = con.cursor()
query = 'select count(*) from p574036'
cur.execute(query)
row = cur.fetchone()
AREA_LIMIT = row[0]
colorsearch()
cur.close()
con.close()
実行
上のプログラムをfourcolor.pyで保存して以下コマンドで実行します。
python fourcolor.py
うまくいけば、color列に1〜4の数値が入ります。
QGISでveg.sqliteを開きcolor列で塗り分ければできあがり！
...となるはずでしたが...いやできるんです...ポリゴンが少なければ。
課題 orz
上記方法は、地図が4色で塗れるまで、総当りでやってみるという方法です。
なので、ポリゴン数が増えると探索時間が爆発して、プログラムが終わりません。
ということで、別の方法を考えたいのですが、ググったらA○○GIS用のプログラムを見つけました。
http://arcscripts.esri.com/details.asp?dbid=14822
これ使ったら1日たっても終わらなかった塗り分けが、あっというまにできました。
詳細は載っていませんが、おそらくグラフ理論だとかを駆使した4色問題の証明方法に関わる
もっと根本的な手法なんでしょう。
ハイ、ここで挫折しました。
参考
http://itpro.nikkeibp.co.jp/article/COLUMN/20070312/264600/
http://www.gaia-gis.it/spatialite-2.4.0-4/splite-python.html
http://www.gesource.jp/programming/python/code/0013.html
http://forums.arcgis.com/threads/29947-quot-Four-Color-a-Map-quot-for-ArcGIS-10
http://arcscripts.esri.com/details.asp?dbid=14822
http://fourcolormap.com/
http://www.vb-helper.com/howto_map_coloring.html
これからこの本読んで解決できるように勉強します。。。
付録
4色にこだわらず、なんとなく少ない数で塗るプログラム
# coding:utf-8
from pyspatialite import dbapi2 as db
con = db.connect("veg_p57.sqlite")
cur = con.cursor()
query = 'select count(*) from p574036'
cur.execute(query)
row = cur.fetchone()
for i in range(row[0]):
query = 'select a.color from p574036 a, p574036 b where b.OGC_FID='
query += '%d ' % (i+1)
query += 'and ST_TOUCHES(a.geometry,b.geometry)'
#print query
cur.execute(query)
m = cur.fetchall()
m.sort()
print m
if m[0][0]>1:
maxcolornum = 1
else:
maxcolornum = m[0][0]
for k in range(1,len(m)):
if maxcolornum + 1 < m[k][0]:
maxcolornum = maxcolornum + 1
break
else:
maxcolornum = m[k][0]
if maxcolornum == m[len(m)-1][0]:
maxcolornum = m[len(m)-1][0]+1
print maxcolornum
query = 'update p574036 set color='
query += '%d ' % maxcolornum
query += 'where p574036.OGC_FID='
query += '%d' % (i+1)
#print query
cur.execute(query)
con.commit()
cur.close()
con.close()
