Encoded polyline を生成するQGISのカスタム関数 #GIS

QGISにはEncoded polyline を生成する機能がないので、カスタム式を作ってみた。

Encoded polyline とは

エンコードポリラインアルゴリズム形式

Googleが策定したGIS上のポリラインをbase64形式の文字列で表現するアルゴリズム。ラインを圧縮した文字列で表現できる。
理屈としては、始点の緯度経度をまず文字列化し、そこからの差分を連結していく。差分を取っているので日本やアメリカなど経度の桁が大きい地域とそうじゃない地域で、同じ圧縮率となる。
詳しくは、下記ドキュメント参照、
エンコードポリラインアルゴリズム形式

Python の Encoded Polyline アルゴリズム

QGISのカスタム式はPythonで記載するが、Encoded Polyline のアルゴリズムそのものは下記のライブラリの式をそのまま使用させていただく。
polyline (PyPI)
なお、このライブラリの元ネタは Mapbox Polyine だそうだ。
@mapbox/polyline

QGISのカスタム関数

実際のソースの説明の前に、関数の作り方と実行方法を簡単に説明しておく

参考資料

カスタム式の作り方は下記参照
QGISユーザガイド - 13. 式でレベルアップ - 13.1. 式
式の中で使用するオブジェクトなどの説明は下記参照
PyQGIS 開発者用 Cookbook

作り方

QGISの式文字列ビルダーで、"関数エディタ"タブを開くとカスタム関数を作成できるようになる。
"+"ボタンで関数のファイルを追加可能。

実行方法

出来上がった関数は組み込み関数と同じように使用できる。
もちろん、プレビューも表示される。

encodedPolyline 関数

今回作成した関数 encodedPolyline の仕様は単純。
地物からEncoded polyline化した文字列を生成し、戻り値とするだけ。引数はなし。
対象となる地物はラインのみで、なおかつシングルパートとする。ライン以外やマルチパートの地物は固定のエラー文字列を戻すことにする。

完成形

まず、完成ソースを確認する。

import io
import itertools
import math
from qgis.core import *
from qgis.gui import *
from typing import List, Tuple

def _pcitr(iterable):
    return zip(iterable, itertools.islice(iterable, 1, None))

def _py2_round(x):
    # The polyline algorithm uses Python 2's way of rounding
    return int(math.copysign(math.floor(math.fabs(x) + 0.5), x))

def _write(output, curr_value, prev_value, factor):
    curr_value = _py2_round(curr_value * factor)
    prev_value = _py2_round(prev_value * factor)
    coord = curr_value - prev_value
    coord <<= 1
    coord = coord if coord >= 0 else ~coord

    while coord >= 0x20:
        output.write(chr((0x20 | (coord & 0x1f)) + 63))
        coord >>= 5

    output.write(chr(coord + 63))

@qgsfunction(group='Custom', usesgeometry=True)
def encodePolyline(feature, parent):
    """
    Encode poliline to encoded polyline string.
    Only accepting single part polyline.
    <h2>Example usage:</h2>
    <ul>
      <li>encodePolyline() -> "`~oia@..."</li>
    </ul>
    """
    geom = feature.geometry()
    geomSingleType = QgsWkbTypes.isSingleType(geom.wkbType())
    res = ""
    if geom.type() == QgsWkbTypes.LineGeometry:
        if geomSingleType:
            l = []
            for pnt in geom.asPolyline():
                l.append((pnt.x(), pnt.y()))
            output, factor = io.StringIO(), int(10 ** 5)
            
            _write(output, l[0][1], 0, factor)
            _write(output, l[0][0], 0, factor)
            
            for prev, curr in _pcitr(l):
                _write(output, curr[1], prev[1], factor)
                _write(output, curr[0], prev[0], factor)
                
            res = output.getvalue()
        else:
            res = "Not Single part"
    
    else:
        res = "Not polyline"
    
    return res

説明

内部関数の _pcitr, _py2_round, _write は前述の polyline ライブラリーをそのままコピーしただけなので、解説は省略する。

シグニチャ

@qgsfunction(group='Custom', usesgeometry=True)
def encodePolyline(feature, parent):
    """
    Encode poliline to encoded polyline string.
    Only accepting single part polyline.
    <h2>Example usage:</h2>
    <ul>
      <li>encodePolyline() -> "`~oia@..."</li>
    </ul>
    """

groupは関数のメニューに表示するグループ今回は"Custom"に表示している。
usesgeometryを true にしておかないと、地物を処理できない。
コメントをHTMLで記載しておくと、関数メニューに説明として表示される。

QgsWkbTypes

geom = feature.geometry()
geomSingleType = QgsWkbTypes.isSingleType(geom.wkbType())
geom.type() == QgsWkbTypes.LineGeometry

引数のfeatureは、対象の地物を表し、geometry()メソッドでジオメトリ情報をしゅとくできる。さらにwkbType()メソッドで地物の種別¹が取得できるので、ラインを表すQgsWkbTypes.LineGeometryのみに絞り込む。
また、QgsWkbTypes.isSingleType()メソッドで地物がシングルパートかどうかの確認を行うことができる。

asPolyline()

geometry()で取得したジオメトリに対し、asPolyline()を使うと QgsPolylineXY² オブジェクトが取得できる。これは、座標を表す QgsPointXY³ のリストである。

x(), y()

QgsPointXYからは、x()とy()でX座標とY座標を取得できる。QGISの式で使用する$X、$Yと同じ値である。
あとは、 polyline ライブラリに順に値を渡していき、最終的に取得した文字列を戻り値とすれば完成
一点注意が必要だが、QGISでは X が経度(longitude/long/lng)で Y が緯度(latitude/lat)である。各種ライブラリでは、lat/lng の順序であることが多いので注意が必要。

_write(output, l[0][1], 0, factor)
_write(output, l[0][0], 0, factor)
            
for prev, curr in _pcitr(l):
    _write(output, curr[1], prev[1], factor)
    _write(output, curr[0], prev[0], factor)

上記の部分で、タプルの0番目(X)と1番目(Y)を入れ替えているのは、 Encoded Polyline の順序が lat/lng だからである。
また、factor = int(10 ** 5)の部分だが、Encoded Polyline の精度が緯度/経度小数点5桁であるため、その数字に丸めるために_write関数に引き渡している。

生成されたEncoded Polyline の確認

Polyline decoder utility

生成された文字列を確認するときは、Google の Polyline decoder utility が便利。
文字列を張り付けると、地図上にラインが描画される。
Polyline decoder utility

例

下記の文字列は日本一短い路線・芝山鉄道の路線を表した文字列⁴
ukiyEghzwYhGcEtRgMrE[pG}ApGgDvDqExE_BvH}EhMuI~D{DdAmCxBgHjBgG
この文字列を、 Polyline decoder utility に入力すると、下記のように路線が表示される。

Encoded polyline を生成するQGISのカスタム関数