[PYTHON] Erstellen eines tky2jgd-Plug-Ins ohne praktische Anwendbarkeit mit QGIS Teil 1

Einführung

Möglicherweise möchten Sie die Koordinaten des alten und des neuen Vermessungssystems konvertieren. Mit ArcGIS kann es schnell und schnell verarbeitet werden. Es scheint, dass praktische Konvertierung mit verschiedenen FOSS4G durchgeführt werden kann. groß.

Aus diesem Grund sagt @ tohka383, dass eine native QGIS-Konvertierung mithilfe von Parameterdateien möglich ist. Nachdem ich unterrichtet wurde, werde ich es wagen, mein eigenes tky2jgd-Plug-In zu erstellen.

Bevor Sie dieses Plug-In erstellen, lesen Sie jedoch die folgenden Elemente Ich denke nicht, dass es praktisch ist, weil es anderen Produkten unterlegen sein wird. ·Geschwindigkeit ·Richtigkeit ・ Vielseitigkeit

Kurz gesagt, was ich sagen möchte, ist Ich möchte keine Koordinatenkonvertierung durchführen, sondern ein Plug-In für die Koordinatenkonvertierung erstellen.

Letzter Zielpunkt

・ Vorwärtskonvertierung (japanisches Geografiesystem → Weltgeografiesystem) ・ Kann unabhängig von Breite / Länge und Ebenenwinkel konvertiert werden ・ Alle Geometrietypen können konvertiert werden Es kann eine Nachkonvertierungsschicht erstellt werden, die den Attributwert erbt. ・ Interpolieren Sie den Korrekturwert der Parameterdatei, um die Genauigkeit zu verbessern. * Schwierigkeitsgrad: Mittel ・ Es können mehrere Ebenen gleichzeitig ausgewählt und konvertiert werden ・ Reverse Conversion (Weltvermessungssystem → Japanisches Vermessungssystem) ist ebenfalls möglich. * Schwierigkeitsgrad: Hoch

Dieses Thema

-Implementiertes Lesen von Par-Datei-Daten ・ Umstellung vom japanischen Geografiesystem (EPSG: 4301) auf das Weltgeografiesystem (EPSG: 4612) -Der Geometrietyp der Zielebene ist nur Multi-Polygon. -Konvertierungsergebnis ist nur Geometrie (keine Attribute) ・ Behandeln Sie keine Ausnahmen

Fang an zu arbeiten

Lass es uns jetzt machen. Erstellen Sie zunächst den Basisteil des Plug mit QGIS Plugin Builder. Ich werde auf die Verwendung verzichten. Wenn Sie also mehr wissen möchten, lesen Sie bitte die folgenden Informationen. Ich habe versucht, ein Python-Plug-In mit QGIS3 Teil 1 zu erstellen. Erstellen Sie eine Basis

Sobald das Basisteil fertig ist, können wir mit der Montage beginnen.

Par-Datei lesen

Lesen Sie zuerst die Par-Datei und erstellen Sie ein Wörterbuch (dict), das den Korrekturwert aus dem Netzcode abrufen kann. Laden Sie die wesentliche Par-Datei aus dem Folgenden herunter. https://www.gsi.go.jp/sokuchikijun/tky2jgd_download.html Entpacken Sie es nach dem Herunterladen und legen Sie es im Plugins-Ordner ab.

Wenn Sie die Par-Datei mit einem Editor öffnen, sieht es so aus.

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

Sie können sehen, dass die ersten beiden Zeilen die Kopfzeileninformationen und die dritte Zeile die Daten fester Länge sind.

Die Lesefunktion sieht wie folgt aus.

    #Par-Datei lesen
    def loadPar(self):
        self.par = {}
        parfile = 'TKY2JGD.par'
        with open(os.path.join(os.path.abspath(os.path.dirname(__file__)), parfile)) as f:
            #Überspringen Sie zwei Zeilen in der Kopfzeile
            skip = 2
            for i in range(skip):
                next(f)

            #Zeile für Zeile lesen und in Liste speichern
            while True:
                line = f.readline()
                if not line:
                    # EOF
                    break
                #Der Wert ist Sekunden, also teilen Sie ihn hier
                self.par[int(line[0:8])] = (float(line[8:18]) / 3600, float(line[18:28]) / 3600)

Das Konvertierungsparameterwörterbuch ist für die Instanzvariable par abgeschlossen.

Wenn Sie beispielsweise auf "self.par [46303582]" zugreifen, wird der Taple "(0.0035549972222222222, -0.0022593166666666667)" zurückgegeben.

Koordinatentransformation implementiert

Der Ablauf der Konvertierungsverarbeitung ist ungefähr so. Ich leihe mir den Code von @ mormor aus, um den Netzcode zu berechnen. Python-Skript zum Konvertieren von Längen- und Breitengraden in Mesh-Code

    #Konvertierungsprozess koordinieren
    def execTrans(self, layer):
        #Ebene nach der Konvertierung generieren
        afterLayer = self.createAfterLayer()

        #Beginnen Sie die Bearbeitung
        afterLayer.startEditing()
        fields = afterLayer.fields()

        inFeat = QgsFeature()
        feats = layer.getFeatures()
        while feats.nextFeature( inFeat ):
            beforeGeom = QgsGeometry( inFeat.geometry() )
            afterGeom = self.moveCorrection(beforeGeom)
            #Bug-Fix<Redundante Transformationen entfernen>
            # afterGeom.transform(self.crsTrans)

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

        #Ende der Bearbeitung
        afterLayer.commitChanges()
        QgsProject.instance().addMapLayer(afterLayer)

    def createAfterLayer(self):
        # EPSG:Erstellen Sie 4612 Speicherebenen (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

    """Für die Breiten- / Längengradkonvertierung ist keine QgsCoordinateTransform erforderlich
    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

Ausführen, wenn Sie die Taste drücken

Platzieren Sie das Button-Objekt in der UI-Datei und Klicken Sie darauf, um die Konvertierung für die aktive Ebene durchzuführen.

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

Überprüfung

Ich denke, ich habe alle Themen erreicht, also werde ich überprüfen, ob ich es wirklich getan habe. Als ich die entsprechende EPSG: 4301-Ebene aktivierte und die Taste drückte, erhöhte sich die Anzahl der falsch ausgerichteten Ebenen.

Vergleichen wir die Scheitelpunktkoordinaten. Extrahieren Sie die Eckpunkte des Polygons, um die Koordinaten eines beliebigen Punkts zu erhalten.

Vergleichen wir es nun mit dem Konvertierungsergebnis von Webversion TKY2JGD --Kokuchi Ritsuin.

Da der Abschlussvorgang nicht ausgeführt wird, liegt ein Fehler vor, der jedoch nahezu identisch ist. (Es gibt auch eine Ziffernrundung in der Attributanzeige von QGIS.) Du hast es geschafft.

Die Quelle ist bisher hier. (2020/01/14 Bugfix Version) https://github.com/ozo360/tky2jgd/releases/tag/ddbc2ba

Nächster Zeitplan

Lassen Sie uns das nächste Mal so etwas implementieren. ・ Kann unabhängig von Breite / Länge und Ebenenwinkel konvertiert werden ・ Alle Geometrietypen können konvertiert werden Es kann eine Nachkonvertierungsschicht erstellt werden, die den Attributwert erbt. ・ Es können mehrere Ebenen gleichzeitig ausgewählt und konvertiert werden

Ich wünsche allen ein schönes Weihnachtsfest.

Lizenz für diesen Artikel

 Dieser Artikel wird unter Creative Commons Attribution 4.0 International License bereitgestellt.