[PYTHON] Créer un plug-in tky2jgd sans praticité avec QGIS Partie 1

introduction

Vous voudrez peut-être convertir les coordonnées des anciens et des nouveaux systèmes d'arpentage. Avec ArcGIS, il peut être traité rapidement et rapidement. Il semble que la [conversion] pratique (https://qiita.com/tohka383/items/e73c7235efc15efe2c1b) puisse être effectuée avec divers FOSS4G. génial.

C'est pourquoi @ tohka383 dit que la conversion native de QGIS à l'aide de fichiers de paramètres est possible. Après avoir appris, j'oserai créer mon propre plug-in tky2jgd.

Cependant, avant de faire ce plug-in, dans les éléments suivants Je ne pense pas que ce soit pratique car il sera inférieur aux autres produits. ·la vitesse ·précision · Polyvalence

Bref, ce que je veux dire c'est Je ne veux pas faire de conversion de coordonnées, je veux créer un plug-in de conversion de coordonnées.

Point cible final

・ Conversion avant (système de géographie japonais → système de géographie mondiale) ・ Peut être converti indépendamment de la latitude / longitude et de l'angle du plan ・ Tous les types de géométrie peuvent être convertis -Une couche de post-conversion qui hérite de la valeur d'attribut peut être créée. ・ Interpole la valeur de correction du fichier de paramètres pour améliorer la précision. * Difficulté: moyenne ・ Plusieurs couches peuvent être sélectionnées et converties à la fois ・ La conversion inverse (système d'enquête mondiale → système d'enquête japonais) est également possible * Difficulté: élevée

Ce thème

-Implémentation de la lecture des données par fichier ・ Conversion du système de géographie japonais (EPSG: 4301) en système de géographie du monde (EPSG: 4612) -Le type de géométrie de la couche cible est uniquement multi-polygone. -Le résultat de la conversion est la géométrie uniquement (pas d'attribut) ・ Ne gérez pas les exceptions

Commence travail

Faisons-le maintenant. Commencez par créer la partie de base du plug-in à l'aide de QGIS Plugin Builder. Je vais omettre comment l'utiliser, donc si vous voulez en savoir plus, veuillez vous référer à ce qui suit. J'ai essayé de créer un plug-in python avec QGIS3 Partie 1 Créer une base

Une fois la pièce de base terminée, nous pouvons commencer le montage.

Lire le fichier par

Tout d'abord, lisez le fichier par et créez un dictionnaire (dict) qui peut obtenir la valeur de correction à partir du code de maillage. Téléchargez le fichier par essentiel à partir de ce qui suit. https://www.gsi.go.jp/sokuchikijun/tky2jgd_download.html Après le téléchargement, décompressez-le et placez-le dans le dossier plugins.

Lorsque vous ouvrez le fichier par avec un éditeur, il ressemble à ceci.

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

Vous pouvez voir que les deux premières lignes sont les informations d'en-tête et la troisième ligne sont les données de longueur fixe.

La fonction de lecture ressemble à ce qui suit.

    #lecture de fichier par
    def loadPar(self):
        self.par = {}
        parfile = 'TKY2JGD.par'
        with open(os.path.join(os.path.abspath(os.path.dirname(__file__)), parfile)) as f:
            #Passer deux lignes dans l'en-tête
            skip = 2
            for i in range(skip):
                next(f)

            #Lire ligne par ligne et stocker dans la liste
            while True:
                line = f.readline()
                if not line:
                    # EOF
                    break
                #La valeur est en secondes, alors divisez-la ici
                self.par[int(line[0:8])] = (float(line[8:18]) / 3600, float(line[18:28]) / 3600)

Le dictionnaire des paramètres de conversion est complété pour la variable d'instance par.

Par exemple, accéder à self.par [46303582] renverra le taple (0.0035549972222222222, -0.0022593166666666667).

Transformation des coordonnées implémentée

Le flux de traitement de conversion est à peu près comme ça. J'emprunte le code de @ mormor pour calculer le code du maillage. Script Python pour convertir la latitude et la longitude en code de maillage

    #Coordonner le processus de conversion
    def execTrans(self, layer):
        #Générer une couche après la conversion
        afterLayer = self.createAfterLayer()

        #Commencer l'édition
        afterLayer.startEditing()
        fields = afterLayer.fields()

        inFeat = QgsFeature()
        feats = layer.getFeatures()
        while feats.nextFeature( inFeat ):
            beforeGeom = QgsGeometry( inFeat.geometry() )
            afterGeom = self.moveCorrection(beforeGeom)
            #Correction d'un bug<Supprimer les transformations redondantes>
            # afterGeom.transform(self.crsTrans)

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

        #Fin de l'édition
        afterLayer.commitChanges()
        QgsProject.instance().addMapLayer(afterLayer)

    def createAfterLayer(self):
        # EPSG:Créer 4612 couches de mémoire (multipolygone)
        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

    """Aucun QgsCoordinateTransform requis pour la conversion de latitude / longitude
    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

Exécuter lorsque vous appuyez sur le bouton

Placez l'objet bouton dans le fichier ui et Cliquez dessus pour effectuer la conversion du calque actif.

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

Vérification

Je pense avoir atteint tous les thèmes, donc je vais vérifier si je l'ai vraiment fait. Lorsque j'ai activé le calque EPSG: 4301 approprié et appuyé sur le bouton, le nombre de calques mal alignés a augmenté.

Comparons les coordonnées des sommets. Extrayez les sommets du polygone pour obtenir les coordonnées de n'importe quel point.

Maintenant, comparons-le avec le résultat de la conversion de Version Web TKY2JGD --Kokuchi Ritsuin.

Étant donné que le processus d'achèvement n'est pas effectué, il y a une erreur, mais c'est presque la même chose. (Il y a aussi un arrondi des chiffres sur l'affichage des attributs de QGIS) Tu l'as fait.

La source jusqu'à présent est ici. (Version de correction de bogue 2020/01/14) https://github.com/ozo360/tky2jgd/releases/tag/ddbc2ba

Calendrier suivant

La prochaine fois, implémentons quelque chose comme ça. ・ Peut être converti indépendamment de la latitude / longitude et de l'angle du plan ・ Tous les types de géométrie peuvent être convertis -Une couche de post-conversion qui hérite de la valeur d'attribut peut être créée. ・ Plusieurs couches peuvent être sélectionnées et converties à la fois

Passez un bon Noël à tous.

Licence pour cet article

 Cet article est fourni sous la Creative Commons Attribution 4.0 International License .