[PYTHON] Make a tky2jgd plugin with no practicality in QGIS Part 1

Introduction

You may want to convert the coordinates of the old and new geodetic datums. With ArcGIS, you can process quickly and quickly, It seems that practical conversion can be done with various FOSS4G. great.

That's why @ tohka383 says that QGIS native conversion using parameter files is possible. After being taught, I will dare to make my own tky2jgd plug-in.

However, before making this plugin, in the following elements I don't think it is practical because it will be inferior to other products. ·speed ·accuracy ・ Versatility

In short, what I want to say is I don't want to do coordinate conversion, I want to make a coordinate conversion plug-in.

Final target point

・ Forward conversion (Japanese geodetic system → World geodetic system) ・ It can be converted regardless of latitude / longitude and plane angle. ・ All geometry types can be converted -A post-conversion layer that inherits the attribute value can be created. -Interpolate the correction value of the parameter file to improve the accuracy * Difficulty: Medium ・ Multiple layers can be selected and converted at once ・ Inverse conversion (world geodetic system → Japanese geodetic system) is also possible * Difficulty: High

This theme

-Implemented reading of par file data ・ Convert from Japanese geodetic system (EPSG: 4301) to world geodetic system (EPSG: 4612) -The geometry type of the target layer is multi-polygon only. -Conversion result is geometry only (no attributes) ・ Do not handle exceptions

Start work

Let's make it now. First, create the base part of the plugin using QGIS Plugin Builder. I will omit how to use this, so if you want to know more, please refer to the following. I tried to make a python plugin with QGIS3 Part 1 Create a base

Once the base part is complete, we can start mounting.

Read par file

First, read the par file and create a dictionary (dict) that can get the correction value from the mesh code. Download the essential par file from the following. https://www.gsi.go.jp/sokuchikijun/tky2jgd_download.html After downloading, unzip it and put it in the plugins folder.

When you open the par file with an editor, it looks like this.

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

You can see that the first two lines are the header information and the third line is the fixed length data.

The read function looks like this:

    #par file reading
    def loadPar(self):
        self.par = {}
        parfile = 'TKY2JGD.par'
        with open(os.path.join(os.path.abspath(os.path.dirname(__file__)), parfile)) as f:
            #Skip two lines in the header
            skip = 2
            for i in range(skip):
                next(f)

            #Read line by line and store in list
            while True:
                line = f.readline()
                if not line:
                    # EOF
                    break
                #The value is seconds, so divide it here
                self.par[int(line[0:8])] = (float(line[8:18]) / 3600, float(line[18:28]) / 3600)

The conversion parameter dictionary is completed for the instance variable par.

For example, accessing self.par [46303582] will return the tuple (0.0035549972222222222, -0.0022593166666666667).

Implemented coordinate transformation

The flow of conversion processing is roughly like this. I borrow @ mormor's code to calculate the mesh code. Python script to convert latitude / longitude to mesh code

    #Coordinate conversion process
    def execTrans(self, layer):
        #Generate layer after conversion
        afterLayer = self.createAfterLayer()

        #Start editing
        afterLayer.startEditing()
        fields = afterLayer.fields()

        inFeat = QgsFeature()
        feats = layer.getFeatures()
        while feats.nextFeature( inFeat ):
            beforeGeom = QgsGeometry( inFeat.geometry() )
            afterGeom = self.moveCorrection(beforeGeom)
            #Bug fixes<Remove redundant transformations>
            # afterGeom.transform(self.crsTrans)

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

        #End of editing
        afterLayer.commitChanges()
        QgsProject.instance().addMapLayer(afterLayer)

    def createAfterLayer(self):
        # EPSG:Create 4612 memory layers (multipolygon)
        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

    """No QgsCoordinateTransform required for latitude / longitude conversion
    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

Execute when you press the button

Place the button object in the ui file and Click on it to perform the conversion on the active layer.

    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')

Verification

I think I've achieved all the themes, so I'll verify if I really did. When I activated the appropriate EPSG: 4301 layer and pressed the button, the number of layers that were misaligned increased.

Let's compare the vertex coordinates. Extract the vertices of the polygon to get the coordinates of any point.

Now, let's compare it with the conversion result of Web version TKY2JGD --Geographical Survey Institute.

Since the complement processing is not performed, there is an error, but it is almost the same. (There is also rounding on the QGIS attribute display) You did it.

The source so far is here. (2020/01/14 Bug fix version) https://github.com/ozo360/tky2jgd/releases/tag/ddbc2ba

Next schedule

Next time, let's implement something like this. ・ It can be converted regardless of latitude / longitude and plane angle. ・ All geometry types can be converted -A post-conversion layer that inherits the attribute value can be created. ・ Multiple layers can be selected and converted at once

Have a nice Christmas, everyone.

License for this article

 This article is provided under Creative Commons Attribution 4.0 International License .