QGIS ボロノイ図でバッファを指定すると発生する不具合を修正してみる

目標

みんな大好きボロノイポリゴン。
GISのチュートリアルなどで商圏分析をやってみましょうとか言って必ずといってもいいほど使われるボロノイポリゴンですが、QGISのジオメトリツールにあるボロノイポリゴンはバッファ領域を設定すると正しいボロノイ図ができません。これでは困るんで治しましょう。

どんな不具合か

バッファをかけていないものとバッファを20%かけて作ったボロノイポリゴンを重ねてみます。

FotoJet.png

矢印の部分で大きくボロノイポリゴンに差分が出てしまっています。
バッファを変えただけでボロノイの線分がずれるなんておかしいですね。

そもそもジオメトリツールって修正できるの?

はい、できます。
ジオメトリツールとは、元々fToolsと呼ばれていた今はProcessingなコアプラグインのことです。Pythonで書かれているのでコンパイルも必要ありません。
普通のプラグイン置き場とは違うところにありますので判りにくいですがQGISのインストールフォルダ内にありますのでエクスプローラーで開きましょう。

C:\Program Files\QGIS 3.0\apps\qgis\python\plugins\processing

この中のファイルを変更すれば修正ができます。

不具合の原因

ボロノイのバッファ領域を指定するとおかしくなるわけですから、バッファを設定している場所を探します。

processing\algs\qgis\VoronoiPolygons.py

こいつが怪しいですね。

VoronoiPolygons.py
        extraX = extent.height() * (buf / 100.0)
        extraY = extent.width() * (buf / 100.0)

90行目にバッファを指定しているか箇所を見つけました。
clip_voronoiにextraXとextraYを渡して処理していますが、おそらくその処理の中で問題が起きていそうです。
198行目から255行目までのどこかでextraX(exX)とextraY(exY)の考慮漏れがあるのではないかと思いますがぱっとは判らないですね。今後調査してみますが、とりあえずはここまでにしておきましょう。

VoronoiPolygons.py
        for edge in edges:
            if edge[1] >= 0 and edge[2] >= 0:
                # Two vertices
                [x1, y1, x2, y2] = clip_line(
                    c.vertices[edge[1]][0],
                    c.vertices[edge[1]][1],
                    c.vertices[edge[2]][0],
                    c.vertices[edge[2]][1],
                    width,
                    height,
                    exX,
                    exY,
                )
            elif edge[1] >= 0:
                # Only one vertex
                if c.lines[edge[0]][1] == 0:
                    # Vertical line
                    xtemp = c.lines[edge[0]][2] / c.lines[edge[0]][0]
                    if c.vertices[edge[1]][1] > (height + exY) / 2:
                        ytemp = height + exY
                    else:
                        ytemp = 0 - exX
                else:
                    xtemp = width + exX
                    ytemp = (c.lines[edge[0]][2] - (width + exX) *
                             c.lines[edge[0]][0]) / c.lines[edge[0]][1]
                [x1, y1, x2, y2] = clip_line(
                    c.vertices[edge[1]][0],
                    c.vertices[edge[1]][1],
                    xtemp,
                    ytemp,
                    width,
                    height,
                    exX,
                    exY,
                )
            elif edge[2] >= 0:
                # Only one vertex
                if c.lines[edge[0]][1] == 0:
                    # Vertical line
                    xtemp = c.lines[edge[0]][2] / c.lines[edge[0]][0]
                    if c.vertices[edge[2]][1] > (height + exY) / 2:
                        ytemp = height + exY
                    else:
                        ytemp = 0.0 - exY
                else:
                    xtemp = 0.0 - exX
                    ytemp = c.lines[edge[0]][2] / c.lines[edge[0]][1]
                [x1, y1, x2, y2] = clip_line(
                    xtemp,
                    ytemp,
                    c.vertices[edge[2]][0],
                    c.vertices[edge[2]][1],
                    width,
                    height,
                    exX,
                    exY,
                )

修正

extraXとextraYの考慮漏れを回避すればいいわけですから、両方の値を0として渡してみましょう。当然ですがそれだと出力範囲が広がらないのでエクステントの領域を広げて渡してやったらどうでしょうか。
試してみましょう。

VoronoiPolygons.py
        # extraX = extent.height() * (buf / 100.0)
        # extraY = extent.width() * (buf / 100.0)
        extraX = 0
        extraY = 0
        extent.setXMinimum(extent.xMinimum() - (extent.height() * (buf / 100.0)))
        extent.setYMinimum(extent.yMinimum() - (extent.width() * (buf / 100.0)))
        extent.setXMaximum(extent.xMaximum() + (extent.height() * (buf / 100.0)))
        extent.setYMaximum(extent.yMaximum() + (extent.width() * (buf / 100.0)))

QGISを再起動してテストしてみます。

修正後.png

思った通りの出力になりました!

検証

Blogopolisから学ぶ計算幾何 第11回 ボロノイ図の作成(前編)を参考に検証していきます。

ボロノイ辺の性質
aとbを隔てるボロノイ辺は,aとbからちょうど距離の等しい点を集めた直線になることが分かります。すなわち,これはaとbの垂直二等分線(perpendicular bisector)ということになります。

では各点の中間地点にポイントを発生させてみましょう。
検証.png
ばっちりボロノイ辺の性質を再現できています。
浮いている点も接していない線分の延長線上にありそうです(適当)。

もしかしたら特定の状況では不具合が出るかもしれませんがこれで修正完了とします。

ちなみに編集する前にバックアップは取りましょう。また、書き込み不可になってますので管理者権限で編集する必要があります。

最後に

ボロノイポリゴンの不具合を修正してみました。
なんだかんだでfTools時代からずーーーとバグってる機能なので公式で今更対応されるということは期待薄かもしれません。
ソースを読んで想像するには、ボロノイポリゴン機能の当初実装時にはバッファ機能がなくて、後からバッファ機能を追加したときに理論を把握できずにごちゃっと実装しちゃったんだろうなってことです。そのくらい冗長なコードですよね。実際この記事で書いた修正方法で良いならメソッド自体随分シンプルに書き直せると思うんですよね。
まあ、理論を把握できていないのは私もなんで人のことは言えませんけどw

場当たり的な対応かもしれませんが少なくとも元の状況よりは良いのではないでしょうか。

こちらの記事でも不具合っぽい書かれ方していますので不具合で間違いないとは思うのですが、ボロノイポリゴンのアリゴリズム自体を把握していないので詳しい方ご指摘よろしくお願いします!

追記:2018/4/11

https://issues.qgis.org/issues/8002
こちらのissueでは生成されるボロノイポリゴンが重なる不具合が報告されています。(なんと5年前から!)
そして3年前に投稿されたパッチ(voronoi.patch)を適用したところ、この記事で指摘している不具合も正常になることが確認できました。
ただ、パッチ投稿時と現行ソースが変わってしまっていることや、extentをバッファサイズに広げることでextraXとextraYを引き回す処理は冗長になっているため、それらを整理したファイルをQGISリポジトリのforkにコミットしました。
以下から改訂版をダウンロードできますので興味がある方は使用してみてください。
https://github.com/ozo360/QGIS/raw/master/python/plugins/processing/algs/qgis/VoronoiPolygons.py

しかし何故パッチファイルまで投稿されているのにずっと放置されているのでしょうね?

以上です。

Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account log in.