はじめに

新旧測地系の座標変換をしたいことありますよね。
ArcGISを使えばさくっと高速に処理できますし、
各種FOSS4Gでも実用的な変換ができるそうです。
すごーい。

そんなわけで@tohka383さんからQGISネイティブでパラメータファイルを使った変換ができると
教えてもらったうえで、あえて自作tky2jgdプラグインを作ってみようと思います。

しかしながらこのプラグインは作る前から以下の要素において
他の製品に劣るものになるので実用性は無いと思います。
・速度
・精度
・汎用性

要するに何が言いたいかというと
座標変換がしたいわけじゃなく、座標変換プラグインが作りたい。

最終的な目標地点

・順変換（日本測地系→世界測地系）できること
・緯度経度、平面直角に関わらず変換できること
・すべてのジオメトリタイプの変換ができること
・属性値を引き継いだ変換後レイヤができること
・パラメータファイルの補正値を補間して精度をあげること ※難易度:中
・複数のレイヤを選択して一度に変換できること
・逆変換（世界測地系→日本測地系）もできること ※難易度:高

今回のお題

・parファイルデータの読み込みを実装
・日本測地系（EPSG:4301）から世界測地系（EPSG:4612）に変換する
・対象レイヤのジオメトリタイプはマルチポリゴンのみ
・変換結果はジオメトリのみ（属性無し）
・例外処理はしない

作業開始

それではさっそく作ってまいりましょう。
まずはプラグインのベース部分をQGIS Plugin Builderを使って作ります。
こいつの使い方は省略しますので詳しく知りたい方は以下を参照してください。
QGIS3でpythonプラグインを作ってみたその1 ベース作成

ベース部分ができたら実装に取り掛かります。

parファイルを読み込み

まずはparファイルを読み込んでメッシュコードから補正値を取得できる辞書(dict)を作ります。
肝心のparファイルは以下からダウンロードします。
https://www.gsi.go.jp/sokuchikijun/tky2jgd_download.html
ダウンロードしたら解凍してプラグインフォルダに入れておきましょう。

parファイルをエディタで開くとこんな感じ。

JGD2000-TokyoDatum Ver.2.1.1
MeshCode   dB(sec)   dL(sec)
46303582  12.79799  -8.13354
46303583  12.79879  -8.13749
46303584  12.79959  -8.14144
46303592  12.79467  -8.13426
46303593  12.79544  -8.13819
46303594  12.79627  -8.14216

先頭の2行がヘッダ情報、3行目から固定長データなのが判りますね。

内容	型	サイズ
メッシュコード	int	8
補正値:dB	float	10
補正値:dL	float	10

読み込み関数は以下のような感じにしました。

    # parファイル読み込み
    def loadPar(self):
        self.par = {}
        parfile = 'TKY2JGD.par'
        with open(os.path.join(os.path.abspath(os.path.dirname(__file__)), parfile)) as f:
            # ヘッダの2行飛ばす
            skip = 2
            for i in range(skip):
                next(f)

            # 1行づつ読み込んでdistに格納
            while True:
                line = f.readline()
                if not line:
                    # EOF
                    break
                # 値は秒なのでここで割っておく
                self.par[int(line[0:8])] = (float(line[8:18]) / 3600, float(line[18:28]) / 3600)

インスタンス変数parに変換パラメータ辞書が完成しました。

例えばself.par[46303582]にアクセスすると(0.0035549972222222222, -0.0022593166666666667)というタプルが返ります。

座標変換を実装

変換処理の流れはざっくりこんな感じ。
メッシュコードの算出は@mormorさんのコードをお借りします。
緯度経度からメッシュコードに変換するPython script

    # 座標変換処理
    def execTrans(self, layer):
        # 変換後レイヤを生成
        afterLayer = self.createAfterLayer()

        # 編集開始
        afterLayer.startEditing()
        fields = afterLayer.fields()

        inFeat = QgsFeature()
        feats = layer.getFeatures()
        while feats.nextFeature( inFeat ):
            beforeGeom = QgsGeometry( inFeat.geometry() )
            afterGeom = self.moveCorrection(beforeGeom)
            # バグ修正 <冗長な変換を削除>
            # afterGeom.transform(self.crsTrans)

            feature = QgsFeature(fields)
            feature.setGeometry(afterGeom)
            afterLayer.addFeature(feature)

        # 編集終了
        afterLayer.commitChanges()
        QgsProject.instance().addMapLayer(afterLayer)

    def createAfterLayer(self):
        # EPSG:4612のメモリレイヤを作成（マルチポリゴン）
        layerUri = "multipolygon?crs=postgis:" + "4612"
        layer = QgsVectorLayer(layerUri, "exchange", "memory")
        return layer

    def moveCorrection(self, geom):
        for i, v in enumerate(geom.vertices()):
            meshcode = self.Coordinate2MeshCode(v.y(), v.x())
            correction = self.par[meshcode]
            geom.moveVertex(v.x() + correction[1], v.y() + correction[0], i)
        return geom

    """ 緯度経度変換時にQgsCoordinateTransformは不要
    def setCrsTrans(self):
        fromCrs = QgsCoordinateReferenceSystem(4301, QgsCoordinateReferenceSystem.EpsgCrsId)
        toCrs = QgsCoordinateReferenceSystem(4612, QgsCoordinateReferenceSystem.EpsgCrsId)
        self.crsTrans = QgsCoordinateTransform(fromCrs, toCrs, QgsProject.instance())
    """

    def Coordinate2MeshCode(self, dLat, dLng ):
        # cf: http://white-bear.info/archives/1400
        # Make sure the input values are decimal 
        iMeshCode_1stMesh_Part_p = dLat *60 // 40
        iMeshCode_1stMesh_Part_u = ( dLng - 100 ) // 1
        iMeshCode_2ndMesh_Part_q = dLat *60 % 40 // 5
        iMeshCode_2ndMesh_Part_v = ( ( dLng - 100 ) % 1 ) * 60 // 7.5
        iMeshCode_3rdMesh_Part_r = dLat *60 % 40 % 5 * 60 // 30
        iMeshCode_3rdMesh_Part_w = ( ( dLng - 100 ) % 1 ) * 60 % 7.5 * 60 // 45
        iMeshCode = iMeshCode_1stMesh_Part_p * 1000000 + iMeshCode_1stMesh_Part_u * 10000 + iMeshCode_2ndMesh_Part_q * 1000 + iMeshCode_2ndMesh_Part_v * 100 + iMeshCode_3rdMesh_Part_r * 10 + iMeshCode_3rdMesh_Part_w
        return iMeshCode

ボタンを押したら実行

uiファイルにボタンオブジェクトを配置して、
それをクリックしたらアクティブレイヤを対象に変換を実行します。

    def initGui(self):
    ...
        self.dlg.btn1.clicked.connect(self.buttonClicked)
    ...

    def buttonClicked(self):
        layer = self.iface.activeLayer()
        crs = layer.crs().postgisSrid()
        if layer.crs().postgisSrid() != 4301:
            return
        self.loadPar()
        # self.setCrsTrans()
        self.execTrans(layer)
        QMessageBox.information(self.dlg, 'tky2jgd', 'finished')

検証

一通りお題を達成したつもりですので本当に出来たのか検証します。
適当なEPSG:4301のレイヤをアクティブにしてボタンを押すと、なんかいい感じにずれたレイヤが増えました。

頂点座標を比較してみましょう。
ポリゴンの頂点を抽出して、任意の点の座標を取得します。

レイヤ	X	Y
変換前	136.863130	35.204678
変換後	136.860179	35.207900

ではWeb版TKY2JGD - 国土地理院の変換結果と比較してみます。

レイヤ	X	Y
プラグイン	136.860179	35.207900
Web版	136.860179203	35.207899483

補完処理をしていませんので誤差が出ていますが大体同じになりました。
（QGISの属性表示上の桁丸めもあります）
やったね。

ここまでのソースはこちらになります。（2020/01/14 不具合修正版）
https://github.com/ozo360/tky2jgd/releases/tag/ddbc2ba

次回の予定

次回はこんなことを実装しよう。
・緯度経度、平面直角に関わらず変換できること
・すべてのジオメトリタイプの変換ができること
・属性値を引き継いだ変換後レイヤができること
・複数のレイヤを選択して一度に変換できること

それでは皆さん、良いクリスマスを。

本記事のライセンス

この記事はクリエイティブ・コモンズ表示 4.0 国際ライセンスの下に提供されています。