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

今回のお題

その２ではプラグイン上のボタンをクリックしたら旧座標の2次元のポイント、ライン、ポリゴンのレイヤが属性値を引き継ぎつつ新座標に変換されるところまで作りました。
しかし、平面直角座標系の変換時に正しく変換されないバグが発生しています。
今回はこちらのバグの原因を調べて修正していきます。

最終目標地点

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

デバッグ

どのタイミングで座標がおかしくなっているのか調べるために座標変換の前後で値をダンプしてみる。

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

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

        inFeat = QgsFeature()
        feats = self.layer.getFeatures()
        cnt = 0
        while feats.nextFeature( inFeat ):
            attributes = inFeat.attributes()
            beforeGeom = QgsGeometry( inFeat.geometry() )
            beforewkt1 = beforeGeom.asWkt() #debug
            if not self.isLonlat:
                beforeGeom.transform(self.trans4301)
            beforewkt2 = beforeGeom.asWkt() #debug
            afterGeom = self.moveCorrection(beforeGeom)
            afterwkt1 = afterGeom.asWkt() #debug
            if not self.isLonlat:
                afterGeom.transform(self.transTo)
            afterwkt2 = afterGeom.asWkt() #debug

            ###############
            #debug
            ###############
            with open('d:/exchange.log', 'a') as f:
                print("-------------旧座標-------------", file=f)
                print(beforewkt1, file=f)
                print("-------------旧座標 → 4301-------------", file=f)
                print(beforewkt2, file=f)
                print("-------------4301 move-------------", file=f)
                print(afterwkt1, file=f)
                print("-------------4301 → 新座標-------------", file=f)
                print(afterwkt2, file=f)

メッシュコードと移動量も吐いておこう。

    def moveCorrection(self, geom):
        for i, v in enumerate(geom.vertices()):
            meshcode = self.Coordinate2MeshCode(v.y(), v.x())
            correction = self.par[meshcode]
            ###############
            #debug
            ###############
            with open('d:/exchange.log', 'a') as f:
                print("####################\n meshcode:" + str(meshcode) + " : " + str(correction), file=f)

            geom.moveVertex(v.x() + correction[1], v.y() + correction[0], i)
        return geom

まずは正常に動いていると思われる緯度経度のサンプルデータを使って検証
シンプルに数ポイントのレイヤにしておく。

実行結果

####################
 meshcode:52365699.0 : (0.0032275777777777776, -0.002950277777777778)
-------------旧座標-------------
Point (136.87317102877543107 35.16588581584142048)
-------------旧座標 → 4301-------------
Point (136.87317102877543107 35.16588581584142048)
-------------4301 move-------------
Point (136.87022075099764606 35.16911339361919886)
-------------4301 → 新座標-------------
Point (136.87022075099764606 35.16911339361919886)
####################
 meshcode:52365750.0 : (0.003231211111111111, -0.002950475)
-------------旧座標-------------
Point (136.87685587271937493 35.13286893949518941)
-------------旧座標 → 4301-------------
Point (136.87685587271937493 35.13286893949518941)
-------------4301 move-------------
Point (136.87390539771936915 35.13610015060630332)
-------------4301 → 新座標-------------
Point (136.87390539771936915 35.13610015060630332)
####################
 meshcode:52365782.0 : (0.0032296583333333334, -0.00295335)
-------------旧座標-------------
Point (136.90695057909175603 35.15078176167821056)
-------------旧座標 → 4301-------------
Point (136.90695057909175603 35.15078176167821056)
-------------4301 move-------------
Point (136.90399722909174329 35.15401142001154255)
-------------4301 → 新座標-------------
Point (136.90399722909174329 35.15401142001154255)

「旧座標 → 4301」と「4301 → 新座標」は処理していないので無視して旧座標と移動後の値を確認する。

136.87317102877543107 - 136.87022075099764606 ≒ 0.002950277778
35.16588581584142048 - 35.16911339361919886 ≒ -0.0032275777777

パラメータファイルの移動量が正しく反映されている。

Web版TKY2JGDでも確認

良い感じだと思います。

では今度はこのサンプルデータを平面直角座標系にして検証します。

####################
 meshcode:52365699.0 : (0.0032275777777777776, -0.002950277777777778)
-------------旧座標-------------
Point (-26732.67109687992706313 -92487.54723248640948441)
-------------旧座標 → 4301-------------
Point (136.87317102877543107 35.16588581584142048)
-------------4301 move-------------
Point (136.87022075099764606 35.16911339361919886)
-------------4301 → 新座標-------------
Point (-27272.61809458389689098 -91780.23878862317360472)
####################
 meshcode:52365750.0 : (0.003231211111111111, -0.002950475)
-------------旧座標-------------
Point (-26407.70557931492294301 -96150.79488501154992264)
-------------旧座標 → 4301-------------
Point (136.87685587271937493 35.13286893949519651)
-------------4301 move-------------
Point (136.87390539771936915 35.13610015060631042)
-------------4301 → 新座標-------------
Point (-26947.79513530126496335 -95443.0704856463271426)
####################
 meshcode:52365782.0 : (0.0032296583333333334, -0.00295335)
-------------旧座標-------------
Point (-23660.26983992193709128 -94171.45945985670550726)
-------------旧座標 → 4301-------------
Point (136.90695057909175603 35.15078176167821766)
-------------4301 move-------------
Point (136.90399722909174329 35.15401142001154966)
-------------4301 → 新座標-------------
Point (-24200.68891517832889804 -93464.00120055227307603)

先ほどの結果と比較してメッシュコードの取得、パラメータファイルの移動量も同じです。
新座標の変換結果をWeb版TKY2JGDで確認

1点目ですが

処理	X	Y
プラグイン	-27272.61809458389689098	-91780.23878862317360472
Web版TKY2JGD	-27003.6517	-92138.6062

あちゃー、全然違ってしまっていますね。
緯度経度の結果と比較すると移動までは正常で新座標への変換でおかしくなっていることが判りました。

何が悪いのか見直してみます。
　！！！
最後の変換は4301ではなく4612から変換しないといけないことに気が付きました。
そりゃズレるわ。

修正

self.transToの変換前CRSを4612に変更するようソースを修正します。

    def setCrsTrans(self):
        crs4301 = QgsCoordinateReferenceSystem(4301, QgsCoordinateReferenceSystem.EpsgCrsId)
        crs4612 = QgsCoordinateReferenceSystem(4612, 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(crs4612, crsTo, QgsProject.instance())

再度実行

####################
 meshcode:52365699.0 : (0.0032275777777777776, -0.002950277777777778)
-------------旧座標-------------
Point (-26732.67109687992706313 -92487.54723248640948441)
-------------旧座標 → 4301-------------
Point (136.87317102877543107 35.16588581584142048)
-------------4301 move-------------
Point (136.87022075099764606 35.16911339361919886)
-------------4612 → 新座標-------------
Point (-27003.54940856946996064 -92138.51321844129415695)
####################
 meshcode:52365750.0 : (0.003231211111111111, -0.002950475)
-------------旧座標-------------
Point (-26407.70557931492294301 -96150.79488501154992264)
-------------旧座標 → 4301-------------
Point (136.87685587271937493 35.13286893949519651)
-------------4301 move-------------
Point (136.87390539771936915 35.13610015060631042)
-------------4612 → 新座標-------------
Point (-26678.68443605181892053 -95801.75775258826615755)
####################
 meshcode:52365782.0 : (0.0032296583333333334, -0.00295335)
-------------旧座標-------------
Point (-23660.26983992193709128 -94171.45945985670550726)
-------------旧座標 → 4301-------------
Point (136.90695057909175603 35.15078176167821766)
-------------4301 move-------------
Point (136.90399722909174329 35.15401142001154966)
-------------4612 → 新座標-------------
Point (-23931.23210011895207572 -93822.46488364967808593)

今度の結果は

処理	X	Y
プラグイン	-27003.54940856946996064	-92138.51321844129415695
Web版TKY2JGD	-27003.6517	-92138.6062

変換後の座標が正しい位置になりました。

続デバッグ

続いてポイント3点繋いでラインにした地物を変換テストします。

####################
 meshcode:52365699.0 : (0.0032275777777777776, -0.002950277777777778)
####################
 meshcode:52365750.0 : (0.003231211111111111, -0.002950475)
####################
 meshcode:52365782.0 : (0.0032296583333333334, -0.00295335)
-------------旧座標-------------
LineString (-26732.67109687992706313 -92487.54723248640948441, -26407.70557931492294301 -96150.79488501154992264, -23660.26983992193709128 -94171.45945985670550726)
-------------旧座標 → 4301-------------
LineString (136.87317102877543107 35.16588581584142048, 136.87685587271937493 35.13286893949519651, 136.90695057909175603 35.15078176167821766)
-------------4301 move-------------
LineString (136.87022075099764606 35.16911339361919886, 136.87390539771936915 35.13610015060631042, 136.90399722909174329 35.15401142001154966)
-------------4612 → 新座標-------------
LineString (-27003.54940856946996064 -92138.51321844129415695, -26678.68443605181892053 -95801.75775258826615755, -23931.23210011895207572 -93822.46488364967808593)

ポイントの変換結果と比較して、問題なく変換できていることが確認できました。

今度はポイント3点繋いでポリゴンにした地物を変換テストします。

####################
 meshcode:52365699.0 : (0.0032275777777777776, -0.002950277777777778)
####################
 meshcode:52365750.0 : (0.003231211111111111, -0.002950475)
####################
 meshcode:52365782.0 : (0.0032296583333333334, -0.00295335)
####################
 meshcode:52366609.0 : (0.003226383333333333, -0.00295045)
-------------旧座標-------------
Polygon ((-26732.67109687992706313 -92487.54723248640948441, -26407.70557931492294301 -96150.79488501154992264, -23660.26983992193709128 -94171.45945985670550726, -26732.67109687992706313 -92487.54723248640948441))
-------------旧座標 → 4301-------------
Polygon ((136.87317102877543107 35.16588581584142048, 136.87685587271937493 35.13286893949519651, 136.90695057909175603 35.15078176167821766, 136.87317102877543107 35.16588581584142048))
-------------4301 move-------------
Polygon ((136.86727030099766012 35.17233977695253344, 136.87390539771936915 35.13610015060631042, 136.90399722909174329 35.15401142001154966, 136.86727030099766012 35.17233977695253344))
-------------4612 → 新座標-------------
Polygon ((-27271.2323161015738151 -91779.79549598139419686, -26678.68443605181892053 -95801.75775258826615755, -23931.23210011895207572 -93822.46488364967808593, -27271.2323161015738151 -91779.79549598139419686))

1点目と最終点がおかしいです。
ポイントの結果と比較してみると移動後の段階でおかしくなっていることが判ります。
値の差を確認するとパラメータファイルの移動量分ずれているようです。

うーん、リングの終点に移動をかけてはいけないのかな？

テストデータのポリゴンにリングを追加してテストしてみます。

####################
 meshcode:52365699.0 : (0.0032275777777777776, -0.002950277777777778)
####################
 meshcode:52365782.0 : (0.0032296583333333334, -0.00295335)
####################
 meshcode:52365750.0 : (0.003231211111111111, -0.002950475)
####################
 meshcode:52366609.0 : (0.003226383333333333, -0.00295045)
####################
 meshcode:52365780.0 : (0.0032288083333333337, -0.0029512249999999996)
####################
 meshcode:52365770.0 : (0.0032296333333333336, -0.0029509749999999998)
####################
 meshcode:52365771.0 : (0.0032301111111111107, -0.002952058333333333)
####################
 meshcode:52365780.0 : (0.0032288083333333337, -0.0029512249999999996)
-------------旧座標-------------
MultiPolygon (((-26732.67109687992706313 -92487.54723248640948441, -23660.26983992193709128 -94171.45945985670550726, -26407.70557931492294301 -96150.79488501154992264, -26732.67109687992706313 -92487.54723248640948441),(-26087.2897302932869934 -93786.75537118859938346, -26087.2897302932869934 -94614.34977696483838372, -25276.93687463738024235 -94441.93427576145040803, -26087.2897302932869934 -93786.75537118859938346)))
-------------旧座標 → 4301-------------
MultiPolygon (((136.87317102877543107 35.16588581584142048, 136.90695057909175603 35.15078176167821766, 136.87685587271937493 35.13286893949519651, 136.87317102877543107 35.16588581584142048),(136.88029758847102357 35.15418997263552825, 136.8803237277067808 35.14672891674221233, 136.8892130654492405 35.14830399731092569, 136.88029758847102357 35.15418997263552825)))
-------------4301 move-------------
MultiPolygon (((136.86727030099766012 35.17233977695253344, 136.90399722909174329 35.15401142001154966, 136.87390539771936915 35.13610015060631042, 136.86727030099766012 35.17233977695253344),(136.8743951384710158 35.16064758930219369, 136.87737275270677628 35.14995855007554582, 136.88626100711590539 35.15153410842203385, 136.8743951384710158 35.16064758930219369)))
-------------4301 → 新座標-------------
MultiPolygon (((-27271.2323161015738151 -91779.79549598139419686, -23931.23210011895207572 -93822.46488364967808593, -26678.68443605181892053 -95801.75775258826615755, -27271.2323161015738151 -91779.79549598139419686),(-26626.05832836601621239 -93078.77546303598501254, -26358.24300081738692825 -94265.33339575032005087, -25547.91871392306347843 -94092.87040038264240138, -26626.05832836601621239 -93078.77546303598501254)))

-------------4301 move-------------
MultiPolygon (((136.86727030099766012 35.17233977695253344, 136.90399722909174329 35.15401142001154966, 136.87390539771936915 35.13610015060631042, 136.86727030099766012 35.17233977695253344),(136.8743951384710158...

リングが２つあってもリングの1点目と最終点が同じズレ方しちゃうことが判明。
ということはポリゴンの時は地物のvertexを全部移動するような処理では駄目だということですね。

リングの1点目と最終点はどちらかは移動処理をしないように修正が必要。

修正

parts()でパートのイテレタ取得して、リング毎のジオメトリに対してvartexを取得しつつ
vertices_beginなんかで判定して飛ばすのかなぁとか妄想しています。
が、具体的にやり方が判らない。

時間がかかりそうなので一旦ここまでで投稿しておく。

つづく