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.
・ 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
-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
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.
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.
Contenu | Moule | Taille |
---|---|---|
Code de maillage | int | 8 |
Valeur de correction:dB | float | 10 |
Valeur de correction:dL | float | 10 |
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)
.
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
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')
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.
couche | X | Y |
---|---|---|
Avant la conversion | 136.863130 | 35.204678 |
Après la conversion | 136.860179 | 35.207900 |
Maintenant, comparons-le avec le résultat de la conversion de Version Web TKY2JGD --Kokuchi Ritsuin.
couche | X | Y |
---|---|---|
Brancher | 136.860179 | 35.207900 |
Version Web | 136.860179203 | 35.207899483 |
É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
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 Creative Commons](https://i.creativecommons.org/l/by/4.0 /88x31.png) Cet article est fourni sous la Creative Commons Attribution 4.0 International License .
Recommended Posts