オープンデータとオープンソースGisを用いたweb上でのインタラクティブ可視化手法について

オープンデータとオープンソースGIS
を用いたWEB上でのインタラクティブ
可視化手法について
株式会社ノーザンシステムサービス
和山亮介
FOSS4G 2014 Tokyo

自己紹介
岩手にある大草原の小さな会社
(株)ノーザンシステムサービス
で
オープンソースとオープンデータとビッグデータ
の御三家と戯れて、たまにソーシャルゲーム作っ
たりしてます。
趣味は旅行と空撮。
最近はドローン空撮にはまってます。
Twitter：@wayama_ryousuke
mixiゲーム

これからお話するのは
GISとかオープンデータにほとんど触れた事がない
人が、WEB上での地図とグラフを用いた可視化シス
テムを構築するときにぶち当たった問題を試行錯誤
しながら解決する時の備忘録的なお話になります。
GISに触れたことがない人がどの様に模索するのか
を参考にしていただければ幸いです。
意訳）間違ったこと言ったら後で突っ込んでください(´・ω・｀)

全てはお客様の
一言から始まりました。

お客様「オープンデータとビッグ
データを組み合わせて集計して
地図で出したい。グラフも連動
で。あとIE7で動くように」
＿人人人人人人人人＿
＞ IE7で動くように＜
￣^Y^Y^Y^Y^Y^Y^Y＾
Y￣

地図？
それGoogleMapsでやれば簡単
じゃね？
グラフ？
それGoogleChartsでやれば簡単
じゃね？

そんなふうに考えていた時期が
俺にもありました。
.i;;;;;;;;;;;;;;;;;;;;i ;,,_,,;;;;;;;;,,..;_ ､,.__..,;_,,,;;;
,,;;;;;;;;__,,___ ....;;;;;;;;;;;;;;;;;
:!;;;;;;;;;;;;;;;;;;ﾌ;'"ﾞﾞﾞﾞﾞﾞﾞﾞﾞ~~￣゛｀!i、￣´ ￣.
｀''‐ i;;;;;;;;;;;;;;;;
. l;;;;;;;;;彡;;;;;ゝ.if'＝=＝==ｰﾞ:: ,.=＝
======r ﾞi;|.l;;;;;;;;;;
ヽ;;ノ;;;;;;;;;;;;;;; .｀''｀-ヽ--''ﾞﾞﾞ;; 'ﾞゝヽ-ﾉ-‐'ﾞ
´ ;.i;;;;;;i.ﾌ;;;;;;l′
ﾞ
l;;;;;;;;;;;;`、,! : ,､
|ll/ ;;;;;;;;r"
i;;;;;;;;;;;;;;l、
/ ;:;: :.゛゛:l ;;;;;;;;./
'';;;;;;;;;;;;;i
i;;;;;;.;:;:;;;;;;;: .i;;;;;;;;;;;;;;;
/;;;ゝ./´::::
´;:;:;;;;;; ,!;;;;;;;;;;;;;;;;;i
ﾐ;;;;'! .ヽ;;､_;.::__::::;;::
/;;;;;;;;;;;;;;;;;;;;;!

Google Mapsは色々制
約がある。
・25000以上のリクエストは有料
・地図は個人利用の範囲でないと印刷できな
い→資料として配布できない
・無料版では非公開サイトでは使えない
・日本の領土に関わる地名表記が日本の主張
通りではない場合がある→官公庁関係の仕事
ではNG
・そもそも市区町村などの境界データがない

よろしい
ならば
オープンデータだ

といいつつ実はほとんどオープン
データに関してはわかってない！
オープンソースもよく使うけどGIS周り
はさっぱり手つかず。
まずは調査してみよう！

オープンデータはどこになにがあるのか？
こちらのまとめサイトが大変参考になった。
作ってくれた人に感謝！
オープンデータ＆公開データのまとめ
http://matome.naver.jp/odai/2135883023685089601?&page=1

オープンデータは何が使えそうか？
政府統計の総合窓口「e-Stat」
http://www.e-stat.go.jp/
国土数値情報ダウンロードサービス
http://nlftp.mlit.go.jp/ksj/
・地図の元データは国勢調査の小地域shpデータや国土数値地図が使えそう。
・市区町村によってはかなり詳細なオープンデータ(街灯のデータ等)もあったが、
かなりばらつきがある。データの内容も様々。国で統計しているデータが一番扱い
やすい。e-statのデータがオススメ。

Qgisにて少地域の人口を5分位で表示基盤地図情報の提供地域（縮尺2500分1相当以上の概ねの範囲）
http://portal.cyberjapan.jp/site/mapuse4/index.html?tab=1&dir=teikyo#tab=1&di
r=teikyo&zoom=5&lat=38.53615&lon=136.36727&layers=BTTTT
・国勢調査のshpデータには全人口は入っているが、男女別五歳階級人口はないため、
csvデータとkey_codeで紐付ける必要があった。
・都道府県、市区町村の行政区域データは小地域のデータをマージして作成した。
・国土地理院の基盤地図情報も考慮したが提供している範囲が国土の全領域ではなかった
ため、断念。
・少地域は町丁・字単位なので市区町村の中で地域間分析（南区と東区での高齢者人口の
割合を知りたいなど）を行う場合の最小単位として活用できる。

オープンソースは何が使えそうか？（1/2）
・クライアントはOpenLayerでいけそう。D3も候補に挙がったが
IE7以降という要件があり断念。
D3はSVG必須＝IE9以降でないと使えないため。ただしFlashが使
えるならsvgwebなどFlash経由のSVG代替ライブラリで行けるぽい。
しかし、遅かったり機能が十分に実装されていないため断念。
・サーバはデータストレージとしてPostGISは確定。
・クライアントにどうやって地図画像を生成するか？
MapServer？GeoServer？
・試した結果、MapServer系はインタラクティブには向かないこと
が分かった。もっさり感あり。あと定義が多すぎて面倒そう。

オープンソースは何が使えそうか？（2/2）
・PostGISから直接GeoJSON形式（JSON形式で書かれたjavascript
用ジオメトリ規格）で出せるぽい。
・GeoJSONは属性値を渡せるため、Openlayersでインタラクティ
ブに操作する場合、相性が良い。
・TopoJSONはトポロジ構造のデータなためデータ容量が少なくて
よいがOpenlayersが対応してないため断念。ただし開発者バージョ
ンのv3では対応している。が、IE7対応ではないため断念。
・UTFGridはtile画像に合わせた2次元行列でピクセル単位での属性
値を保持することができる。データ量が少なく、インタラクティブ
にマウスカーソル上の属性値を表示したりなど期待したが、tile画像
に合わさる＝ズーム値が整数となり、今回はズーム値が小数点の倍
率なため断念。

地図はOpenlayers＋GeoJSONでよさそう
グラフはどうか？
調べるとわんさか出てくる…
D3(svg)
GoogleCharts(apiでflash)
jqChart(canvas)
OpenFlashChart(falsh)
FusionCharts(flash) ・・・etc

・絞り込み条件として
インタラクティブである
アニメーションできる
ie7で動く
noプラグイン
グラフの種類が豊富
矩形選択ができる
…
を追加した結果…
該当件数は1件です

＿人人人人人人人人人＿
＞その名はamCharts＜
￣Y^Y^Y^Y^Y^Y^Y^Y^Y
￣
http://www.amcharts.com/inspiration/motion-chart/

http://www.amcharts.com/demos
/
amChartsとは？
・指定できるパラメータやグラフの種類が多く、イベント関連もかなり豊富で、インタラク
ティブな操作にうってつけなjavascriptグラフライブラリ。
・ほとんど要求仕様のグラフはあるが箱ひげ図など特殊なグラフが対応していない。
→が、何とか既存の組み合わせでできることが判明。
・グラフと地図の連動もいけるようだ。
・無償版と有償版があり、無償版はロゴが入る以外は普通に使える。有償版はロゴが消え、
サポートのインシデントがつく。
・有料($140)だがサポート有でいろいろ開発者に要望を出せる。バグとか要望、質問を実際
に聞いてみたらすぐに返事が来た。ただしもちろん英語。

http://dev.openlayers.org/examples/zoom.html http://www.amcharts.com/dem
os/
http://s.kitazaki.name/docs/geojs
on-spec-ja.html
Openlayers＋
GeoJSON＋
amChartsで本開発に
入ることが決定！

グラフ＆地図
連携デモ

開発中はいくつか
小さい問題はあったが
ほぼ問題なく進んだ…

…ように見えたが、
テスト工程で
大問題が発覚！

IE7で
重すぎて
IEとOSを道連れに落ちる
自爆テロ案件が発生

問題の箇所は地物を選択して保存する機能。
原因はIE8以前はSVGに対応していなく、
VMLで地図やグラフを描画していたが、オブジェクト数が多
くなると比例して急激に重くなることが判明。
リリースまであと2週間…
どないしよ…

といいつつ現状分析してみる。
Geojsonのポリゴン数が多いのか？
↓
ポリゴン軽くすればよいのか？
↓
PostGISのST_Simplifyとか使ってみる
↓
隣接地物とくっつかない。隙間が開いて見た目がやば
い
↓
ならばトポロジーか？
↓
属性値が全部消えて地物の判別がつかない！
↓
簡略化前の地物と比較して領域が含まれているならそ
の地物IDを割り振ればよいか？
↓
きちんと属性値がとれた！
ポリゴンは小さくなるが
スカスカでヤバイ！
ぴったりくっついていい
感じ！

参考１：トポロジーの作り方（yellow_73さんのページを参考にしました）
http://d.hatena.ne.jp/yellow_73/20120323
【STEP1】
トポロジ化したいジオメトリを単純化する(MultiPolygon→Polygon)
①pgis_pref_dumpテーブル作成
CREATE TABLE pgis_pref_dump(
gid SERIAL PRIMARY KEY,
ken varchar(2),
geom GEOMETRY(POLYGON,4326),
area DOUBLE PRECISION
);
②単純なポリゴン化したものをpgis_pref_dumpテーブルに入れる
INSERT INTO pgis_pref_dump (ken,geom)
SELECT ken,(ST_Dump(geom)).geom FROM pgis_pref WHERE cache_flag='2'
③面積を入れる
UPDATE pgis_pref_dump
SET area = st_area(geom ::geography)
【STEP2】
トポロジ1(ポリゴン)を作成する
①pgis_pref_topo1を作成
SELECT topology.CreateTopology('pgis_pref_topo1',4326);
→スキーマにpgis_pref_topo1が追加される
②pgis_pref_dumpのポリゴンデータをpgis_pref_topo1にいれる
SELECT TopoGeo_AddPolygon('pgis_pref_topo1',geom,0)
FROM pgis_pref_dump
WHERE area >= 100000000;
【STEP3】
トポロジ2(簡略化したラインストリング)を作成する
①pgis_pref_topo2を作成
SELECT topology.CreateTopology('pgis_pref_topo2',4326);
→スキーマにpgis_pref_topo2が追加される
②pgis_pref_topo1の簡素化したラインストリングデータをpgis_pref_topo2にいれる
SELECT TopoGeo_AddLineString('pgis_pref_topo2', ST_SimplifyPreserveTopology(geom,
1)) FROM pgis_pref_topo1.edge_data;
【STEP4】
トポロジ2をポリゴン化する
SELECT topology.Polygonize('pgis_pref_topo2');
【STEP5】
トポジオメトリの作成
①topoTempテーブルを作成
CREATE TABLE topotemp (
gid SERIAL PRIMARY KEY,
face_id INT
);
②トポジオメトリカラムを作成
SELECT AddTopoGeometryColumn('pgis_pref_topo2', 'public', 'topotemp',
'topogeom', 'POLYGON');
③topoTempテーブルにトポジオメトリを追加
INSERT INTO topotemp(topogeom)
SELECT topology.CreateTopoGeom('pgis_pref_topo2',3,1,Q.faces)
FROM (
SELECT ARRAY[ARRAY[f.face_id,3 ] ] as faces FROM
pgis_pref_topo2.face AS f
WHERE NOT mbr IS NULL
) AS Q;
【STEP6】
トポジオメトリをジオメトリにしたテーブルを作成
CREATE TABLE topoPref AS
SELECT gid, topogeom::geometry FROM topoTemp;
→これでqgisなどで確認できる

参考２：トポロジーへの属性値の割り当て方法
①岩手県と接しているトポジオメトリを取得。(ST_Intersects)
②接している部分の面積を求める。(ST_Area)
③面積が1以上のトポジオメトリが岩手県に該当する。
面積：
0.015550
面積：
1.483030
面積：
0.00523
面積：
0.01423
【対象のジオメトリと接する面積が1以上のトポジオメトリを取得する】
SELECT f2.gid, f2.topogeom
FROM
(SELECT geom FROM pgis_pref WHERE ken='03' AND cache_flag='2') AS
f1,
(SELECT gid, topogeom FROM topopref) AS f2
WHERE ST_Intersects(f1.geom, f2.topogeom)
AND ST_Area(ST_Intersection(f1.geom, f2.topogeom)) >= 1

簡略化＆トポロジ化で多少は改善したがまだ重い！
↓
札幌市など小地域が極端に多いところで選択するとフリーズする
↓
オブジェクト数が多い＝Geojsonで地物ごと描画しているからか?
↓
地物ではなく画像ならよいかも？
↓
画像でインタラクティブな選択はできないから、選択している範囲を渡して、画像作れば
よいか?
↓
選択画像はリアルタイムで。背景画像はtileで行ってみよう
↓
PostGISでラスタ画像作成できるのでやってみたが枠線の調整等細かい調整ができない＆
Tile画像が生成できないため断念
↓
qgisのpythonインターフェースもやってみたが少し調べた感じだと地物単位の画像生成
ができないため断念
↓
Webkitを使ったwkhtmltoimageでサーバ上でレンダリングさせる方法も試したが、そも
そもchromeでも１０秒位かかっていたので、落ちないけど遅いので断念
↓
地図画像のリアルタイム生成＋Tile画像生成で
よいのがないか探してみると…

リアルタイムで地図画像を生成できる。
Tile画像も作れる。
そう、mapnik ならね
http://ja.wikipedia.org/wiki/%E3%82%B9%E3%83%86%E3%82%A3%E3%83%BC%E3%83%96
%E3%83%BB%E3%82%B8%E3%83%A7%E3%83%96%E3%82%BA#mediaviewer/File:Steve_J
obs_WWDC07.jpg

Mapnikとは？
http://mapnik.org/
・mapnik は地図描画ツールキットでOpenStreetMapのtileレンダリン
グエンジンとかで有名。
・tile生成はmapnikでいけそう！
・リアルタイムで選択範囲の地物ID＋地物単位の画像を作成もmapnik
＋PostGISでいけそう
・さらに複数のレイヤーを合成したり、表示Styleをかなり弄れるため
OpenLayersでできないor クライアントのwebブラウザ(IE7とか)で非
常に時間がかかるような処理をmapnikに投げてやることで処理時間が
大幅に改善できることが判明。
・mapnik先生のお陰で何とかリリースに間に合ったよ＼(^o^)／
http://upload.wikimedia.org/wikipedia/commons/thumb/c/c8/Ljubljana-OpenStreetMap-Mapnik-
100k.svg/741px-Ljubljana-OpenStreetMap-Mapnik-100k.svg.png

mapnik時代の幕開けである
http://mapnik.org/maps/

「それmapnikでできるよ」
が社内流行語に(嘘)
・本来であればQGISやGrassなどデスクトップGISで行う作業もmapnik＋
PostGISでできるようになった。
例1：地図生成の際に地名を表示する際に通常は重心に設置するが、場所に
よっては密集したり、飛び地などがあり、手動で調整する必要がある。そこで
PostGIS＋openlayersを組み合わせて直感的にラベル位置を修正し、画像地図
をmapnikで自動生成するツールを作成した。
ラベルのポイントデータ座標の変更を直感的に行うツール
密集していて分かりづらいドラッグでポイントすっきりして分かりやすくなった
データを移動する
離れた場所は線で表現する
(drawControlで線を描画)

例2：介護事業所からの到達圏域と65歳以上人口のメッシュを
重ねることで、介護事業所のカバレッジを分析したいといった
場合、普通はArcGIS、もしくはQGIS or Grass + pgRouting
を使ったりするが、mapnik +pgRoutingだとWEB上でインタ
ラクティブな操作が可能に。
・用意するデータは国土数値データの道路情報（国道と高速道
路しかない為、正確に分析したい場合は国土地理院の基盤地図
情報が必要）と介護事業所の緯度経度データ（介護サービス情
報公表システムより取得）
・馬鹿正直に計算するとかなり遅くなり現実的ではなくなるの
で介護事業所から50分圏内までを計算範囲にしたり、ノードを
最適化したり、ノード間での範囲内のみにしたりと高速化した
(Fig.1-2)

【SQL】特定のジオメトリから指定した距離内にあるジオメトリを取得
例： IDが24037のジオメトリから1km圏内のジオメトリを取得
SELECT sample.id, sample.geom
FROM
sample,
(select geom from sample where id='24037') as target
WHERE ST_DWithin(sample.geom, target.geom, 0.0089831666666667)
事業所
事業所最短ノード
ID：
24037
Fig. 1 道路ノード最適化
経路検索は、事業所に一番近いノードを使用して計算を行う。
事業所の最短ノード数が多いと計算時間が遅くなる。
→密集している地域は1つにまとめて、計算するノードを数を減らす。

Fig. 2 ざっくりとした枝刈り
pgRoutingの最短経路を求める関数「pgr_dijkstra」は全ての
ノードを探索するため、計算時間が遅くなる。
→探索する範囲を狭めて計算時間を短縮する。
start ID:4371
end ID:9556
【SQL】2つの地点を囲むバウンディングボックス(矩形)を作成
SELECT ST_SetSRID(ST_Expand(ST_Extent(geom),0.1),4326) AS bbox
FROM node
WHERE id IN (4371,9556)
start地点、end地点を囲む矩形
の範囲を作成し、その範囲内で
ルート検索を行う。

・高速化のおかげで
１７サービスの事業所（9997） x ノード数（26105）
で416時間かかっていた計算が15時間程度になり約２６倍高速化した。
・しかし検証用(広島県の隣接のみ)なのでノード数と事業所数は少ない。
詳細な検索をしたい場合はノード数と事業所数が膨大になる可能性が
ある。
・到達圏域はある意味、交通グラフの最短経路問題なのでグラフ最適
化のアルゴリズムが使えるか？
2-Hop LabelingやHub Labelingとか要検討。HadoopやGPGPUなど
で並列化しやすそう。
調査中！
・この到達圏域データをGDAL(ラスターデータフォーマット用ライブ
ラリ)に食わせてごにょごにょ(Fig.3)すると到達圏域ヒートマップtile
画像が出来上がる。

Fig. 3 GDALでヒートマップを作成
ＳＴＥＰ１： gdal_gridでデータ補間
事業所⇔ノード間の最短距離(drivetime_cost)を補間し、tif画像を作成
ＳＴＥＰ２： gdalwarpでtif画像を切り取る
ＳＴＥＰ1で作成したtif画像を都道府県の形に切り取る
ＳＴＥＰ３： gdaldemで画像に色付けをする
ＳＴＥＰ2で切り取ったtif画像をカラーテーブルファイルを使用して色付けをする
カラーテーブルファイル
-999 0 0 0 0
-1 0 0 0 0
5.77664 0 193 0 255
11.43873 145 224 145 255
17.10082 255 242 138 255
22.76291 253 141 141 255
28.425 221 60 79 255
7分位で定義

・お次は65歳以上人口の500mメッシュ作成。
・本当はズバリこのデータがあればよいけど、
無いので国勢調査の小地域毎の男女5歳階級
別人口のcsvデータから500mメッシュに面積
按分して人口分布の7分位コロプレスtile画像
を作成(Fig.4)

Fig. 4 500mメッシュ単位の高齢者人口を求めたい
高齢者の分布をメッシュで表したいが、メッシュ単位ではデータ
をもっていない。
→面積按分にて算出。
【SQL 】
SELECT mesh.id,
sum(mergearea.value * ST_Area(ST_Intersection(mesh.geom, mergearea.geom))
/ ST_Area(mergearea.geom))
FROM
(SELECT geom,
CASE when senior_population in ('-','X') then 0
when senior_population is null then 0
else cast(senior_population as integer)
END AS value
FROM pgis_smallarea_merge
)AS mergearea,
pgis_kokudo_mesh AS mesh
WHERE ST_Intersects(mesh.geom, mergearea.geom)
GROUP BY mesh.id
面積：1420k㎡
65歳以上の人口：18人
ここの部分の高齢者人口を求めたい
重なっている部分の面積：
154k㎡
≪面積按分の求め方≫
65歳以上の人口=
(154k㎡ × 18人) ÷ 1420k
㎡ ≒ 1.9人

この二つの画像を
重ねることで
分析できる！？
＋ =?

到達圏域カバレッジ
（OpenLayers版）
デモ

透明度で重ね合わせる
とよく分からない
(´・ω・｀)

・レイヤーの合成具合を変えることで見やすいようにできるか?
・OpenLayersではレイヤーの合成は透明度以外出来ない。
・ここまで来たなら答えは簡単。
mapnik名言集：
OpenLayers で
出来ないなら
mapnikで
やればいいじゃない。

到達圏域カバレッジ
（mapnik版）デモ

・mapnikではPhotoshopなみに
様々な合成手法が可能。
・合成手法を変えたり、混じり具
合を変えたりなんかもインタラク
ティブに調整可能。
・速度もかなり早い。

オープンデータ＆オープンソースGISを
使ってみた所感(1/2)
オープンデータは・・・
・行政関係のお仕事では北方四島、尖閣諸島、竹島は人がいなくても
表示しないと色々まずい。
・国勢調査のshpデータには竹島や尖閣諸島はあるが北方四島はない
ため、国土数値地図とマージした。
・行政関係は特に出典・原典を気にするのでデータの出典・原典は大
事。民間のデータはほぼダメ。データの作成日時も大事。出典・原典
と作成日時は画面上に明示する必要がある。
・実際にデータを見てみると合併周りが結構大変。
・合併は最新の合併状況に合わせて過去データを合算して集計するこ
とで対処した。
・国勢調査の中には水上生活者という分類もあり、湖や港なども領域
としてあり、手作業で消さなければならなかった。

オープンデータ＆オープンソースGISを
使ってみた所感(2/2)
オープンソースGISは・・・
・英語のサイトが主である為、ニュアンスで理解しなければならず大
変だった。
→Google翻訳をフル活用。Google先生に足を向けて寝れない・・・
・ドキュメントが少ない(特にmapnik)
→Stack Overflow等、Q&Aサイトで検索したりした。
・状況に応じて適切に空間参照系を変換しなければならず苦労した。
→OpenLayersで選択した範囲の画像をmapnikで作成する場合、
OpenlayersがEPSG: 900913（Googleメルカトル）なので、
OpenLayersの選択範囲をいったんEPSG: 4326 （WGS84）の座標
に変換し、PostGIS（EPSG: 4326）で計算した結果をまたEPSG:
900913に変換して画像出力。
EPSG: 900913 はTMS (Tile Map Service)で使われている。

まとめ
・いままでweb上で地図を見せる場合は既に作成済みのtile画像や
簡単なジオメトリを作成したりプロットしたりすることが多かっ
た。
・クライアントやWEBブラウザの性能の向上によりデスクトップ
用GISソフト並みのインタラクティブな操作や分析が可能になって
きていることが分かった。
・またレガシー環境でもクライアント及びサーバの協働によりレ
スポンスを確保することで、より大勢の人がインタラクティブに
オープンデータに触れることで、予想もできないような知見が生
まれるのではないかと期待している。

ご清聴ありがとうございました。

オープンデータとオープンソースGisを用いたweb上でのインタラクティブ可視化手法について

More Related Content

オープンデータとオープンソースGisを用いたweb上でのインタラクティブ可視化手法について