QGISで実用性皆無なtky2jgdプラグインを作る その2

今回のお題

その１ではプラグイン上のボタンをクリックしたら
EPSG:4301なマルチポリゴンレイヤをEPSG:4612に座標変換できるプラグインを作りました。

しかしマルチポリゴンレイヤ限定では使い物になりませんし
平面直角座標系のレイヤにも対応していませんので
今回はその辺を強化して少し実用性を持たせてみたいと思います。

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

複数レイヤを一括変換するのはUI系の話になるので切り分けることにしました。
それでは頑張っていきましょう。

最終目標地点

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

作業開始

その１で作ったソースが無い方はリポジトリから落として来てね。
https://github.com/ozo360/tky2jgd/releases/tag/ddbc2ba

すべてのジオメトリタイプに対応

すべてに対応するとは言ってない。
QGISが扱えるベクトルレイヤのジオメトリタイプはQgsWkbTypes::geometryType
更にそのタイプとしてQgsWkbTypes::Typeがあります。

このプラグインの対象は、ポイント、ライン、ポリゴンで2次元データに限定します。
これで大半のデータをカバーできるはず。
CurveやTriangleのデータちゃんと変換できるのかは知らないので作った後でテストしてみます。
駄目なら除外するように変更しよう。

変換元のレイヤwkbTypeを取得して、新規レイヤのuriに埋め込むには以下のようにします。
例外処理はしませんが対象を限定するための分岐は入れておきます。

layer = self.iface.activeLayer()
wkbType = layer.wkbType()
# point line polygonを対象とする
if not QgsWkbTypes.geometryType(wkbType) in [
        QgsWkbTypes.PointGeometry, 
        QgsWkbTypes.LineGeometry, 
        QgsWkbTypes.PolygonGeometry]:
    # 対象外のgeometryType
    self.iface.messageBar().pushMessage("tky2jgd <b>{}</b>".format('activeLayer is out of scope geometryType'), Qgis.Warning)

# 2次元レイヤを対象とする
if QgsWkbTypes.hasZ(wkbType) or QgsWkbTypes.hasM(wkbType):
    # 対象外のwkbType
    self.iface.messageBar().pushMessage("tky2jgd <b>{}</b>".format('activeLayer is out of scope wkbType'), Qgis.Warning)

self.wkbType = QgsWkbTypes.displayString(wkbType)

layerUri = self.wkbType +"?crs=postgis:4612"
newlayer = QgsVectorLayer(layerUri, "exchange", "memory")

enumの文字列を取得する方法（displayString）が判らなくてとても時間がかかりましたが、なんとか対象レイヤのジオメトリタイプに対応したレイヤを作ることができました。

緯度経度、平面直角に関わらず変換

日本測地系が対象なので以下のような変換ができる必要があります。

系	変換元SRID	変換先SRID
	4301	4612
I	30161	2443
II	30162	2444
III	30163	2445
IV	30164	2446
・	・	・
・	・	・
・	・	・
XVIII	30178	2460
XIX	30179	2461

また、パラメータファイルはEPSG:4301における補正値を表していますので
平面直角は緯度経度に変換してから補正値移動をしないといけません。
よって変換先SRIDへの変換は、いずれのSRIDでもEPSG:4301からの変換となります。
それぞれ座標変換クラスを生成します。

layer = self.iface.activeLayer()
srid = layer.crs().postgisSrid()
if srid == 4301:
    self.isLonlat = True
    self.toSrid = 4612
elif (srid > 30160 and srid < 30180):
    self.isLonlat = False
    self.toSrid = 2442 + srid - 30160
else:
    # 対象外のSRID
    self.iface.messageBar().pushMessage("tky2jgd <b>{}</b>".format('activeLayer is not Tokyo Datum'), Qgis.Warning)

# 座標変換クラスを生成
crs4301 = QgsCoordinateReferenceSystem(4301, QgsCoordinateReferenceSystem.EpsgCrsId)
if not self.isLonlat:
    crsFrom = QgsCoordinateReferenceSystem(fromSrid, QgsCoordinateReferenceSystem.EpsgCrsId)
    self.trans4301 = QgsCoordinateTransform(crsFrom, crs4301, QgsProject.instance())

crsTo = QgsCoordinateReferenceSystem(self.toSrid, QgsCoordinateReferenceSystem.EpsgCrsId)
self.transTo = QgsCoordinateTransform(crs4301, crsTo, QgsProject.instance())

座標変換クラスさえできていれば変換処理時に平面直角なら緯度経度に変換する前処理入れるだけですね。

属性値を引き継ぐ

属性値を引き継ぐためには新レイヤにカラムを追加する必要があります。
元レイヤのカラムを新レイヤにコピーします。

# newLayerを新レイヤとした場合
layer = self.iface.activeLayer()
newLayer.dataProvider().addAttributes(layer.fields())

地物を作る時に、元の属性値を引き継ぎます。

inFeat = QgsFeature()
feats = layer.getFeatures()
while feats.nextFeature( inFeat ):
    attributes = inFeat.attributes()
    geometry = QgsGeometry( inFeat.geometry() )
    # 本当はここで座標変換する
    feature = QgsFeature(fields)
    feature.setAttributes(attributes)
    feature.setGeometry(geometry)

    newLayer.addFeature(feature)

プラグインに埋め込んでいく

各処理ができましたのでプラグインに埋め込んでいきます。

###################################
# 処理開始ボタンクリック
###################################
def buttonClicked(self):
    self.layer = self.iface.activeLayer()
    if not self.isScope():
        return

    self.wkbType = QgsWkbTypes.displayString(self.layer.wkbType())

    # 座標変換パラメータ辞書を生成
    self.loadPar()
    # 座標変換クラスを生成
    self.setCrsTrans()

    # 変換処理
    self.execTrans()
    QMessageBox.information(self.dlg, 'tky2jgd', 'finished')

###################################
# 変換パラメータファイルを読み込み
###################################
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)

###################################
# 変換対象レイヤかどうかを判定
# （ついでに変換前後のsridと緯度経度フラグを保持する）
###################################
def isScope(self):
    try:
        if not isinstance(self.layer, QgsVectorLayer):
            raise ValueError("activeLayer is not QgsVectorLayer")

        self.srid = self.layer.crs().postgisSrid()
        if self.srid == 4301:
            self.isLonlat = True
            self.toSrid = 4612
        elif (self.srid > 30160 and self.srid < 30180):
            self.isLonlat = False
            self.toSrid = 2442 + self.srid - 30160
        else:
            # 対象外のSRID
            raise ValueError("activeLayer is not Tokyo Datum")

        wkbType = self.layer.wkbType()
        # point line polygonを対象とする
        if not QgsWkbTypes.geometryType(wkbType) in [
                QgsWkbTypes.PointGeometry, 
                QgsWkbTypes.LineGeometry, 
                QgsWkbTypes.PolygonGeometry]:
            # 対象外のgeometryType
            raise ValueError("activeLayer is out of scope geometryType")

        # 2次元レイヤを対象とする
        if QgsWkbTypes.hasZ(wkbType) or QgsWkbTypes.hasM(wkbType):
            # 対象外のwkbType
            raise ValueError("activeLayer is out of scope wkbType")

    except ValueError as e:
        # 対象外
        self.iface.messageBar().pushMessage("tky2jgd <b>{}</b>".format(e), Qgis.Warning)
        return False
    else:
        # 対象
        return True

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

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

    inFeat = QgsFeature()
    feats = self.layer.getFeatures()
    while feats.nextFeature( inFeat ):
        attributes = inFeat.attributes()
        beforeGeom = QgsGeometry( inFeat.geometry() )
        if not self.isLonlat:
            beforeGeom.transform(self.trans4301)
        afterGeom = self.moveCorrection(beforeGeom)
        if not self.isLonlat:
            afterGeom.transform(self.transTo)

        feature = QgsFeature(fields)
        feature.setAttributes(attributes)
        feature.setGeometry(afterGeom)

        afterLayer.addFeature(feature)

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

###################################
# 変換後のメモリレイヤを作成
###################################
def createAfterLayer(self):
    layerUri = self.wkbType +"?crs=postgis:" + str(self.toSrid)
    afterLayer = QgsVectorLayer(layerUri, "exchange", "memory")
    afterLayer.dataProvider().addAttributes(self.layer.fields())
    return afterLayer

###################################
# 地物の頂点移動
###################################
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

###################################
# 座標変換クラスを生成
###################################
def setCrsTrans(self):
    crs4301 = QgsCoordinateReferenceSystem(4301, QgsCoordinateReferenceSystem.EpsgCrsId)
    if not self.isLonlat:
        crsFrom = QgsCoordinateReferenceSystem(self.srid, QgsCoordinateReferenceSystem.EpsgCrsId)
        self.trans4301 = QgsCoordinateTransform(crsFrom, crs4301, QgsProject.instance())

    crsTo = QgsCoordinateReferenceSystem(self.toSrid, QgsCoordinateReferenceSystem.EpsgCrsId)
    self.transTo = QgsCoordinateTransform(crs4301, crsTo, QgsProject.instance())

ちゃんとクラス設計から始めないから、なんだかごちゃっとしてしまいましたね。
この辺りは最後にクラスに切り出して整理するとして、まずは動くことを優先していきます。
isScope関数が死にたくなるほどダサいのも目をつぶる。

検証

地味にサンプルデータを作るのが面倒です。

自前で作ったサンプルデータで平面直角でも変換できることを検証してみます。
１．スクラッチレイヤで日本測地系第16系 EPSG:30176のラインレイヤを作成。
２．OSM読み込んで西表島に落書き。
３．tky2jgdプラグインで変換。

変換後のレイヤのSRIDがちゃんと2458になってます。
© OpenStreetMap contributors

４．ポイント化してWeb版TKY2JGD - 国土地理院を使って確認。

レイヤ	X	Y
変換前	-21165	-177520
変換後	-21313	-177066
Web版TKY2JGD	-21315.4052	-177083.9327

Y方向の誤差がすごいけど多分計算間違ってたらこんな誤差で済まないはずなので多分大丈夫じゃないかな？

引き続き東京でもやってみよう。
併せてジオメトリタイプの検証も行う。
１．スクラッチレイヤで日本測地系第9系 EPSG:30169で今度はマルチポイントレイヤを作成。
２．皇居付近にマルチポイントを書いて。
３．tky2jgdプラグインで変換。

マルチポイントのレイヤも変換出来て、今度のレイヤはSRIDが2451になってます。
© OpenStreetMap contributors

レイヤ	X	Y
変換前	-6984	-34730
変換後	-7276	-34371
Web版TKY2JGD	-7277.1503	-34374.4408

なんか嫌な予感がしないでもないが一旦目を瞑って次へ進む。

今度は属性が引き継がれているかどうかを検証
メモリレイヤに文字型、整数型、浮動小数点のカラムを追加してから地物を追加。
属性値を与えておこう。
変換後のレイヤの属性を表示するとちゃんと保持されている。

これで今回の課題としていたすべての機能は実装できた。

ここまでのソースはこちら
※ただし盛大にバグっている
https://github.com/ozo360/tky2jgd/releases/tag/sono2

嫌な予感の正体

変換後のレイヤをShapefileに保存してみる。
もちろん指定するCRSは変換後のレイヤのものとする。
保存したShapefileを再度QGISに読み込んでみるとズレる。
オンザフライで元レイヤに近い場所に乗らないといけないのに変換後の場所に乗る。
座標を調べるとめっちゃズレてる。
旧座標を新座標に変換した値を更に新座標に移動したような値。

プロジェクトへの追加方法が間違ってるのかもしれないし、それ以前に根本的なミスを犯しているのかもしれない。
本来であれば不具合を解消してからポストすべきだと思うけどまあこういう試行錯誤があってもいいじゃないですか。
というわけで冬眠に入るのでしばらく更新できないが次回はこの不具合を調べる。

追記 2020/01/14

緯度経度の変換でパラメータファイルで移動後、更にQgsCoordinateTransformで旧座標から新座標に変換をかけているのは間違いだったので「その１」も含めて修正。
これにより緯度経度の変換は正常になった。
平面直角の変換時にパラメータファイルの移動量の倍移動してしまっていることが判った。
こちらの原因はまだ見つけられていない。

本記事のライセンス

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