waigani's diary

QGISを中心にFOSS4Gをいじくる

PostGISのST_AsMVTでバイナリベクトルタイル

これはベクトルタイル Advent Calendar 2017 12/22の記事です。

MIERUNEとベクトルタイル

OpenMapTiles

株式会社MIERUNEでは、MIERUNE地図としてタイル地図画像を配信しています。
f:id:waigani:20171220125421p:plain
何をベースに作っているかというとOpenMapTilesになります。
OpenMapTilesの作成手順に沿っているので、間でmbtilesを作成していて、バイナリベクトルタイルになっています。
諸事情によりバイナリベクトルタイルでの公開していないのですが、近々公開したいとは思っています。

改良したい点

MIERUNE地図として公開している範囲としては、日本+台湾になります。
本年度中に範囲をどんどん増やしていくつもりだったのですが、ストップしています。
なぜストップしているかというと、巨大なmbtilesを作るのが時間はかかるし、手間がかかるからです。
Klokan Technologiesのエンジニアにも聞いたことがあるのですが、「分散する以外にないガンバレ」との返答をいただきました。そうですよね、それしかないですよね。
postserveも出てきたので、mbtilesすっとばしてPostGISからの直接配信の方向にも力をいれていくのかなと思ったら、「あれKlokan Technologiesの人じゃないんだよね」とも言っていたので、やっぱり頑張ってmbtiles作るのかな。。。
おまけにスキーマ変更のコストも高いんです。Debugどうしてるの?と聞いたら「興味深いトピックスだよね(笑」で終わったので、エンジニアの地道な努力に支えられているのが現状です。

MIERUNEとしては、ベースマップに関するデータだけではなく、変更が頻繁にあるようなデータについてもベクトルタイルで配信していくことを想定しておきたいです。SQLだけでほぼスキーマ定義を変更することも目的に、PostGISからアクセスに応じて都度データを取得しての配信についても色々お試し中です。

ST_AsMVT

いろいろアプリケーションは出てきているのですが、もう1段低レベルなところから確認しましょう。
ST_AsMVTを使ってみます。
https://postgis.net/docs/ST_AsMVT.html

インストール

PostGIS ver2.4から、使えるようになっています。
ただし、libprotobuf-c が必要です。PostGIS ver2.4ですぐ使えるのかなとpgadminで関数一覧をみても、
f:id:waigani:20171220132921p:plain
st_asmvtgeomはあるのに、st_asmvtがない、となるかもしれません。
PostGISを、ソースからインストールしましょう。configure時に、protobu-c supportがされていることを確認しておきます。
f:id:waigani:20171220133541p:plain

tileから経緯度

tileのxyzから経緯度を求める関数は必要です。OpenStreetMap wikiの記事"Slippy_map_tilenames"に、plpgsqlでの例が載っていますので、使わせてもらいましょう。
http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#PostgreSQL

使い方例

  • テーブルjapanに、国土数値情報の行政区域が入っています
  • x:3638, y:1618, z:12で切り出します
  • st_makeenvelope()でタイル矩形を作る
  • st_intersection()でタイル矩形でjapanを切り出し
  • st_asmvtgeom()でタイル内座標(4096)に変換
  • st_asmvt()でベクトルタイル化

の例になります。

select st_asmvt(v, 'test', 4096, 'geom')
from (
select j.n03_001, j.n03_002, j.n03_003, j.n03_004, j.n03_007, st_asmvtgeom(st_intersection( j.geom, e.geom), e.geom,  4096, 0, true)  as geom
from
 japan j,
 (select st_makeenvelope(tile2lon(3638,12), tile2lat(1618,12), tile2lon(3639,12), tile2lat(1619,12), 4326) as geom) e
 where
  st_intersects(j.geom, e.geom)
) v

確認

上記SQLをtest.sqlというファイル名で保存しておいたとして、SQLの実行結果を一旦ファイルに書き出し。書き出したファイルをバイナリ化してから、geojsonに変換して中身を確認してみます。バイナリベクトルタイルからgeojsonへの変換は、vt2geojsonを使います。

psql -f test.sql -tAz0 "データベース名" > 12-3638-1618.txt
xxd -r -p 12-3638-1618.txt > 12-3638-1618.pbf
vt2geojson -x 3638 -y 1618 -z 12 --layer test 12-3638-1618.pbf > 12-3638-1618.geojson

geojsonをQGISに表示して、タイル番号と一緒に表示して確認してみます。
意図したタイルのデータを、切り出し出来ています。
f:id:waigani:20171220135757p:plain

補足(2017.12.28追記)

誤解を与えそうなので。mapbox vector tile specificationでは、タイルの外側にバッファを取ります。
上記例だと、きれいにタイルに合わせた絵にしようと思ってst_intersection()で切り取っているのですが、ここは本来要らないです。

  st_asmvtgeom(j.geom, e.geom,  4096, 256, true)

のようにして、4番目の引数がバッファサイズ、5番目の引数がclipするか、を指定してあげます。
ただ、st_asmvtgeomに渡してあげるj.geomは、このタイルに重なってくる範囲に限定してあげましょう。

ここで求人

MIERUNEで、一緒にベクトルタイルしたいエンジニアいませんか?
表立って求人は出していないのですが、ご興味のある方がいらっしゃればいつでもお声掛けください。
札幌に事務所はありますが、どこで勤務して頂いても構いません。

CS立体図関係で今年やったことと、やり残したこと

FOSS4G Advent Calendar 2017

qiita.com

今年はブログの更新もサボっていて、db analytics showcaseと共催させていただいたFOSS4G Hokkaidoの事を何も書いていないとか、8月のQGIS Conferenceの事を書いていないとか、やってないことだらけです。

Advent Calendarネタとしても、QGIS3.0の新機能紹介をしたかったのですが、試そう試そうと思っているうちに時間がなくなってしまったので、ダメですね...
年末年始に時間できたら、触ってどこかに書きます。

今年、中途半端でも成果として公開できたのがCS立体図なので、その振り返りだけ。

CS立体図関係で作ったもの

1月につくばで行われた、つくばFOSS4Gもくもく会・CS立体図勉強会に誘われて、ノコノコ出かけていったのが始まりでした。

兵庫県さんがまとめてくれていたQGISでのCS立体図作成マニュアルと、長野県 戸田様からの直々のご指導により、processingで処理をつなげればできそうなことは、その場で確認できました。その後も戸田様にアドバイスいただきながら、QGISプラグインとしてまとめて、公開まで持っていきました。
戸田様に直に会い、作り方を習えたのは良い機会でした。

G空間情報センターから作成したQGISプラグインへのリンク

G空間情報センターからgithubにリンクしています。
www.geospatial.jp

QGIS Python Plugins Repository

8月のQGIS Conferenceに参加した機会に、QGISプラグインリポジトリにも登録しました。
QGISのメニューからインストール可能になっています。
QGIS Python Plugins Repository
f:id:waigani:20171203222225p:plain

おまけツール

作成したQGISプラグインは、GeoTIFFからのCS立体図作成を前提にしています。
FOSS4G Tokyoの際に、x,y,z座標+α入ったテキストファイルからの作成としたいとの相談を受けました。FOSS4G Kansaiのハンズオンで使いたいとの要望でしたので、急ぎで取り掛かり。QGISプラグインには組み込まず、簡単なQGIS Processingスクリプトを作成しています。
テキストデータのうち、x,y,zがどのカラムに入っているかを指定して、GeoTIFFまで変換します。
github.com

ただ、なにぶん突貫工事だったので、テストも不十分な状態です。
どなたかデータを提供いただければ、色々試します。

積み残し

色々積み残してしまっているのですが、年末年始対応できるかな。もしくは来年に持ち越し。

QGIS Processing スクリプトへの変更

Githubにissueも上がっているのですが、CS立体図作成のようなデータ解析機能については、プラグインでの作成ではなくProcessingに入れるのがQGISの流れです。
そんなに大変な変更ではないのですが、さぼっている状態なので、そのうち時間を作って手を付けなくては。
github.com

QGIS 3.0への対応

これは早目に行わなくてはと思ってはいるのですが...
3.0対応のついでに、Processingスクリプトに移行しようかな。

おまけツールの改良とテスト

戸田様から、

  • 区切り文字の対応の確認
  • メッシュサイズの指定

をあげてもらっているのですが、こちらも手付かず。

以上

QGISで人口が半分になるように2分割してみた

人口が半分になるように都道府県を二分割してみた

togetter.com
というのがtwitterでまわっていましたね。
これをQGISでやるプラグインを作って、FOSS4G Advent Calendar
FOSS4G Advent Calendar 2016 - Qiita
のネタにしようと思っていました。ただ朝の時点で出来ていなかったので、他の投稿をしておきました。
UIを作るところで思いの外ハマったのでですが、なんとなくできたので晒しておきます。

QGISプラグイン

GitHubにおいておきました。
github.com

プラグインの動作

インストールして実行してもらうと、ベクトルレイヤ(ポリゴンのみ)を選択→レイヤ内の使用するカラムの選択を行えます。
合わせてスライダーで割合を指定します。
f:id:waigani:20161201224121p:plain
この画像の例だと、

  • h22ka01204レイヤを
  • JINKO カラムを使って
  • 人口上位の地域から加えていって30%になる地域:それ以外の地域で色分け(30%超えたところで止めてるのでだいたいね)

になります。

e-stat から平成22年の小地域をダウンロードしてきて試してみました。

  • 旭川市を、人口上位の地域から数えていって30%になるところまでを青:それ以外を赤で色分け

f:id:waigani:20161201224811p:plain

  • 札幌市を、人口上位の地域から数えていって50%になるところまでを青:それ以外を赤で色分け

f:id:waigani:20161201225042p:plain

処理の仕方

大事なところだけ。

  • comboboxとsliderから値を持ってきて
  • 指定のlayerの指定のカラムの値を全部足す
  • sliderの値の割合をかけて、目標の数字とする
  • 指定カラムの値をソート
  • ソートしたものを足していって、目標の数字を超えたところでストップ
  • この条件でQgsRuleBasedRendererV2を作成
  • layerにセットして再描画

としています。

QtDesinger with QGIS 2.14 custom widgets

どこでハマってたかというと。気づいてなかったのですがQGISをインストールすると、QtDesingerがついてくるんですね。
このQtDesingerにはQGISのcustom widgetsがついています。

例えばQgsFieldComboBoxだと、レイヤを渡すと勝手にFieldでComboBoxを作ってくれます。
ただ作ったUIを実行しようとすると、エラーが出てしまいました。このあたりを参考に修正しています。
osgeo-org.1560.x6.nabble.com

QGIS 3.0を後押ししよう

FOSS4G Advent Calendar 2016

今年もFOSS4G Advent Calendarを行います。
例年のような勧誘もしておりませんので、まだまだ参加に余裕がある状況です。
どのソフトに興味があるとか、この処理で困っているとか、些細なことで結構ですので、みなさまのご参加お待ちしております。
qiita.com

Let’s make a big funding push for QGIS 3.0!

QGISみなさん使ってますよね?
お仕事で各所回っていても、現場で広まっているのがわかるようになってきました。"QGIS"と提案して、すでに使ってる/検討しているというレベルでお話が通じるようになって来ています。

さて、そこでみなさんにお願い。

QGIS projectからQGIS 3.0リリースへ向けて、みなさまからの後押しのお願いがきております。
本来なら下記ブログを翻訳して各国で紹介してね、と言われているのですが、時間がなくて放置しててすいません。。。
blog.qgis.org

かいつまんで、

  • 2017年にQGIS 3.0のリリースが予定されています
  • 2.xシリーズで解決できなかった多くの問題を取り除くために、モダンな構造となります
  • 例えば、Qt5、Python3への移行、古いコードの見直し、多くのQGISにパフォーマンスの改善と安定性をもたらす調整を行います
  • その他のキーとなる問題としては.... (省略)

みなさまへのお願いとしては、

  • 個人ユーザの方には寄付の検討

Donations

  • 企業で働いている方は、スポンサーの検討

Sponsorship Program Overview

  • 社内にQGIS開発者がいる場合は、QGIS3.0開発への貢献に時間を割くことへの考慮を

などなどとなっていますので、ぜひご検討ください。

日本でも、QGIS3.0リリースの前には日本語環境での動作検証や、ドキュメントの翻訳など、さまざまな作業が必要になると思われます。この辺はOSGeo.JP内でも調整して追々。

みなさま、ご協力の程よろしくお願い致します。

合同会社MIERUNEを設立いたしました

このたび合同会社MIERUNEを設立いたしました。
www.mierune.co.jp

オープンなデータ、オープンなソフトウェアを駆使して、みなさまと一緒に課題の解決に取り組みます。
設立メンバーは地理空間情報分野で特に豊富な経験を持っておりますが、そこにこだわらず広く、可視化、データ解析、データ共有の技術を駆使したコンサルティング、システム開発を行っていきます。
事務所は札幌に置いておりますが、どこへでも行きます(もしくはネットで繋げます)ので、気軽に声をお掛けください。

まだ設立したばかりの会社で、必要最低限の設備でのスタートとなっております。
amazonで欲しいものリストも公開してみました。MIERUNEへの応援よろしくお願い致します。
www.amazon.co.jp

QGISプログラミング入門 2016Osaka編

QGISプログラミング入門

2013年に札幌と東京でハンズオンを行った際に資料を作成しました。
slideshareにあげてあるのですが、地味にいまでもアクセスあるんですよね。

大阪市大でワークショップの機会をいただきましたので、資料を最近のバージョンに合わせて修正いたしました。
内容は基本変更してないです。参考資料を紹介で済ましているデバックの部分や、プロセッシングの補足を充実させたいところなのですが、またの機会に気が向いたら。

GNOSIS Cartographer(まだまだα版)を使ってみたよ

GNOSIS

対応していただきありがとうございました

ecere.ca

GNOSISを試してみたいと思い、FOSS4G Seoulの後から開発者のJeromeさんと連絡を取ってみました。私のラップトップPCが貧弱なためインストールできず、ご迷惑をかけております。試してみたいというだけの問い合わせに、二か月くらいかかって対応していただくという手間を取らせてしまい、本当に申し訳ありませんでした。Jeromeさんとってもいい人です。

まだまだ開発中ということでα版ではありますが、GNOSIS Cartographerを動かしてみます。

まずはサンプルデータを表示

サンプルデータとして"12 Globe Gores by Giovanni Maria Cassini, Published in Rome, 1790"をいただいたので表示してみます。
Globo terrestre / Cassini, Giovanni Maria / 1790

つづいて同じくサンプルとしていただいた、SRTMを表示してみます。
DEMとして扱って陰影・段彩を設定できます。

基盤地図情報を表示

shapefile,GeoTIFFはそのまま読み込めるということなので、基盤地図情報から札幌市周辺をダウンロードして表示してみました。5mDEM(容量約200MBくらい)と基盤地図情報から水域・水涯線・海岸線・軌道の中心線・道路構成線・道路縁を入れています。


引き続き

まずは触ってみたところまででした。
タイル地図読み込めるのかとか、地理院地図標高タイル読み込める、とかとか聞いてみようかな。
SDK使うと何できるかも知りたいかな。