$30 off During Our Annual Pro Sale. View Details »

[FOSS4G 2014 Tokyo] オープンデータとオープンソースGISを用いたWEB上で...

[FOSS4G 2014 Tokyo] オープンデータとオープンソースGISを用いたWEB上でのインタラクティブ可視化手法について

★ この資料は FOSS4G 2014 Tokyo の登壇発表資料です ★

全てはお客様の一言から始まりました。
お客様「オープンデータとビッグデータを組み合わせて集計して地図で出したい。グラフも連動で。あとIE7で動くように

この発表は、GISとかオープンデータにほとんど触れた事がない人が、WEB上での地図とグラフを用いた可視化システムを構築するときにぶち当たった問題を試行錯誤しながら解決していった、闘いの記録である...。

発表者: 和山 亮介(株式会社ノーザンシステムサービス 研究開発部)

More Decks by ノーザンシステムサービス | Northern System Services

Other Decks in Technology

Transcript

  1. そんなふうに考えていた時期が 俺にもありました。 .i;;;;;;;;;;;;;;;;;;;;i ;,,_,,;;;;;;;;,,..;_ 、 ,.__..,;_,,,;;;,,;;;;;;;;__,,___ ....;;;;;;;;;;;;;;;;; :!;;;;;;;;;;;;;;;;;;フ;'"゙゙゙゙゙゙゙゙゙~~ ̄゛ `!i、  ̄´

     ̄ .`''‐ i;;;;;;;;;;;;;;;; . l;;;;;;;;;彡;;;;;ゝ .if'=====ー゙ :: ,.========r ゙i;|.l;;;;;;;;;; ヽ;;ノ;;;;;;;;;;;;;;; .`''`-ヽ--''゙゙゙ ;; '゙ゝヽ-ノ-‐'゙´ ;.i;;;;;;i.フ;;;;;;l′ ゙l;;;;;;;;;;;;`、 ,! : ,、|ll/ ;;;;;;;;r" i;;;;;;;;;;;;;;l、 / ;:;: :.゛゛:l ;;;;;;;;./ '';;;;;;;;;;;;;i i;;;;;;.;:;:;;;;;;;: .i;;;;;;;;;;;;;;; /;;;ゝ ./´:::: ´;:;:;;;;;; ,!;;;;;;;;;;;;;;;;;i ミ;;;;'! .ヽ;;、_;.::__::::;;:: /;;;;;;;;;;;;;;;;;;;;;! : 、.-、/;;;|l `;:;:;: ,.ノ.::: :lく;;;;;;;;l゙゙′ ''";;;;;;;;;;;;;'l、 .,, :::::::::___ ・ ,'" :::::: .ヽ;;;;;丶; ;:;;;;;;;;;;;;;;;;;;ヽ, ゙゙--= ゙̄~゛`''>,,._,..,r;" ,,l゙ :::::::::::: i;;;;;; .`";'"゙;;;;;;;;;;;;;i、 ヽ_ ゙̄ ̄゛_、 __r::::::::::::::::: ヽ ;;.'ミ;;;;;;;;;:'.、 ゙゙゙゙"'''"~ ‐"":::::::::::::::::::: i
  2. といいつつ現状分析してみる。 Geojsonのポリゴン数が多いのか? ↓ ポリゴン軽くすればよいのか? ↓ PostGISのST_Simplifyとか使ってみる ↓ 隣接地物とくっつかない。隙間が開いて見た目がやば い ↓

    ならばトポロジーか? ↓ 属性値が全部消えて地物の判別がつかない! ↓ 簡略化前の地物と比較して領域が含まれているならそ の地物IDを割り振ればよいか? ↓ きちんと属性値がとれた! ポリゴンは小さくなるが スカスカでヤバイ! ぴったりくっついていい 感じ!
  3. 【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などで確認できる 参考1:トポロジーの作り方(yellow_73さんのページを参考にしました) http://d.hatena.ne.jp/yellow_73/20120323
  4. ①岩手県と接しているトポジオメトリを取得。(ST_Intersects) ②接している部分の面積を求める。(ST_Area) ③面積が1以上のトポジオメトリが岩手県に該当する。 面積: 1.483030 面積: 0.015550 面積: 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 参考2:トポロジーへの属性値の割り当て方法
  5. 簡略化&トポロジ化で多少は改善したがまだ重い! ↓ 札幌市など小地域が極端に多いところで選択するとフリーズする ↓ オブジェクト数が多い=Geojsonで地物ごと描画しているからか? ↓ 地物ではなく画像ならよいかも? ↓ 画像でインタラクティブな選択はできないから、選択している範囲を渡して、画像作れば よいか?

    ↓ 選択画像はリアルタイムで。背景画像はtileで行ってみよう ↓ PostGISでラスタ画像作成できるのでやってみたが枠線の調整等細かい調整ができない& Tile画像が生成できないため断念 ↓ qgisのpythonインターフェースもやってみたが少し調べた感じだと地物単位の画像生成 ができないため断念 ↓ Webkitを使ったwkhtmltoimageでサーバ上でレンダリングさせる方法も試したが、そも そもchromeでも10秒位かかっていたので、落ちないけど遅いので断念 ↓ 地図画像のリアルタイム生成+Tile画像生成で よいのがないか探してみると…
  6. Mapnikとは? ・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 http://mapnik.org/
  7. 例2:介護事業所からの到達圏域と65歳以上人口のメッシュを 重ねることで、介護事業所のカバレッジを分析したいといった 場合、普通はArcGIS、もしくはQGIS or Grass + pgRouting を使ったりするが、mapnik +pgRoutingだとWEB上でインタ ラクティブな操作が可能に。

    ・用意するデータは国土数値データの道路情報(国道と高速道 路しかない為、正確に分析したい場合は国土地理院の基盤地図 情報が必要)と介護事業所の緯度経度データ(介護サービス情 報公表システムより取得) ・馬鹿正直に計算するとかなり遅くなり現実的ではなくなるの で介護事業所から50分圏内までを計算範囲にしたり、ノードを 最適化したり、ノード間での範囲内のみにしたりと高速化した (Fig.1-2)
  8. 事業所 事業所最短ノード 【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つにまとめて、計算するノードを数を減らす。
  9. 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地点を囲む矩形の 範囲を作成し、その範囲内で ルート検索を行う。
  10. -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 カラーテーブルファイル Fig. 3 GDALでヒートマップを作成 STEP1: gdal_gridでデータ補間 事業所⇔ノード間の最短距離(drivetime_cost)を補間し、tif画像を作成 STEP2: gdalwarpでtif画像を切り取る STEP1で作成したtif画像を都道府県の形に切り取る STEP3: gdaldemで画像に色付けをする STEP2で切り取ったtif画像をカラーテーブルファイルを使用して色付けをする 7分位で定義
  11. 【 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 Fig. 4 500mメッシュ単位の高齢者人口を求めたい 高齢者の分布をメッシュで表したいが、メッシュ単位ではデータ をもっていない。 →面積按分にて算出。 ここの部分の高齢者人口を求めたい 面積:1420k㎡ 65歳以上の人口:18人 重なっている部分の面積: 154k㎡ ≪面積按分の求め方≫ 65歳以上の人口= (154k㎡ × 18人) ÷ 1420k ㎡ ≒ 1.9人
  12. オープンデータ&オープンソース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)で使われている。