Make a tky2jgd plugin with no practicality in QGIS Part 2

This theme

Part 1 If you click the button on the plugin I made a plug-in that can convert the coordinates of EPSG: 4301 multi-polygon layer to EPSG: 4612.

However, it is useless only for multi-polygon layers. Since it does not support layers in the plane orthogonal coordinate system This time, I would like to strengthen that area and make it a little more practical.

・ All geometry types can be converted ・ It can be converted regardless of latitude / longitude and plane angle. -A post-conversion layer that inherits the attribute value can be created. ~~ ・ You can select multiple layers and convert them at once ~~ ・ Do not handle exceptions

Since it is a UI-related matter to convert multiple layers at once, I decided to separate them. Let's do our best.

Final target point

・ Forward conversion (Japanese geodetic system → World geodetic system) What you can do (finished) </ font> ・ 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

Start work

If you don't have the source created in Part 1, drop it from the repository. https://github.com/ozo360/tky2jgd/releases/tag/ddbc2ba

Compatible with all geometry types

I'm not saying that it corresponds to all. The vector layer geometry types that QGIS can handle are QgsWkbTypes :: geometryType Another type is QgsWkbTypes :: Type.

The target of this plugin is limited to 2D data with points, lines, and polygons. This should cover most of the data. I don't know if it can convert Curve and Triangle data properly, so I will test it after making it. If it doesn't work, change it to exclude it.

To get the conversion source layer wkbType and embed it in the uri of the new layer, do as follows. Exception handling is not performed, but a branch is included to limit the target.

layer = self.iface.activeLayer()
wkbType = layer.wkbType()
#target line polygon
if not QgsWkbTypes.geometryType(wkbType) in [
        QgsWkbTypes.PointGeometry, 
        QgsWkbTypes.LineGeometry, 
        QgsWkbTypes.PolygonGeometry]:
    #Not applicable geometryType
    self.iface.messageBar().pushMessage("tky2jgd <b>{}</b>".format('activeLayer is out of scope geometryType'), Qgis.Warning)

#Target 2D layers
if QgsWkbTypes.hasZ(wkbType) or QgsWkbTypes.hasM(wkbType):
    #Not applicable wkbType
    self.iface.messageBar().pushMessage("tky2jgd <b>{}</b>".format('activeLayer is out of scope wkbType'), Qgis.Warning)

self.wkbType = QgsWkbTypes.displayString(wkbType)

layerUri = self.wkbType +"?crs=postgis:4612"
newlayer = QgsVectorLayer(layerUri, "exchange", "memory")

It took a long time because I didn't know how to get the enum string (displayString), but I managed to create a layer that corresponds to the geometry type of the target layer.

Convert regardless of latitude / longitude or plane orthogonality

Since the target is the Japanese geodetic system, it is necessary to be able to perform the following conversions.

system Conversion source SRID Conversion destination SRID
4301 4612
I 30161 2443
II 30162 2444
III 30163 2445
IV 30164 2446
XVIII 30178 2460
XIX 30179 2461

Also, the parameter file shows the correction value in EPSG: 4301. The plane right angle must be converted to latitude and longitude before moving the correction value. Therefore, the conversion to the conversion destination SRID will be from EPSG: 4301 for any SRID. Generate a coordinate transformation class for each.

layer = self.iface.activeLayer()
srid = layer.crs().postgisSrid()
if srid == 4301:
    self.isLonlat = True
    self.toSrid = 4612
elif (srid > 30160 and srid < 30180):
    self.isLonlat = False
    self.toSrid = 2442 + srid - 30160
else:
    #Not applicable SRID
    self.iface.messageBar().pushMessage("tky2jgd <b>{}</b>".format('activeLayer is not Tokyo Datum'), Qgis.Warning)

#Generate coordinate transformation class
crs4301 = QgsCoordinateReferenceSystem(4301, QgsCoordinateReferenceSystem.EpsgCrsId)
if not self.isLonlat:
    crsFrom = QgsCoordinateReferenceSystem(fromSrid, QgsCoordinateReferenceSystem.EpsgCrsId)
    self.trans4301 = QgsCoordinateTransform(crsFrom, crs4301, QgsProject.instance())

crsTo = QgsCoordinateReferenceSystem(self.toSrid, QgsCoordinateReferenceSystem.EpsgCrsId)
self.transTo = QgsCoordinateTransform(crs4301, crsTo, QgsProject.instance())

As long as the coordinate conversion class is created, if the plane is perpendicular to the plane during the conversion process, the preprocessing for converting to latitude and longitude is just included.

Inherit attribute value

You need to add a column to the new layer to inherit the attribute values. Copy the columns from the original layer to the new layer.

#When newLayer is a new layer
layer = self.iface.activeLayer()
newLayer.dataProvider().addAttributes(layer.fields())

When making a feature, the original attribute value is inherited.

inFeat = QgsFeature()
feats = layer.getFeatures()
while feats.nextFeature( inFeat ):
    attributes = inFeat.attributes()
    geometry = QgsGeometry( inFeat.geometry() )
    #Actually, the coordinates are converted here
    feature = QgsFeature(fields)
    feature.setAttributes(attributes)
    feature.setGeometry(geometry)
    
    newLayer.addFeature(feature)

Embed in plugins

Now that each process has been completed, we will embed it in the plug-in.

###################################
#Click the process start button
###################################
def buttonClicked(self):
    self.layer = self.iface.activeLayer()
    if not self.isScope():
        return

    self.wkbType = QgsWkbTypes.displayString(self.layer.wkbType())

    #Generate coordinate conversion parameter dictionary
    self.loadPar()
    #Generate coordinate transformation class
    self.setCrsTrans()

    #Conversion process
    self.execTrans()
    QMessageBox.information(self.dlg, 'tky2jgd', 'finished')
    
###################################
#Read conversion parameter file
###################################
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)

###################################
#Determine if it is a conversion target layer
#(By the way, keep the srid and latitude / longitude flags before and after conversion)
###################################
def isScope(self):
    try:
        if not isinstance(self.layer, QgsVectorLayer):
            raise ValueError("activeLayer is not QgsVectorLayer")

        self.srid = self.layer.crs().postgisSrid()
        if self.srid == 4301:
            self.isLonlat = True
            self.toSrid = 4612
        elif (self.srid > 30160 and self.srid < 30180):
            self.isLonlat = False
            self.toSrid = 2442 + self.srid - 30160
        else:
            #Not applicable SRID
            raise ValueError("activeLayer is not Tokyo Datum")

        wkbType = self.layer.wkbType()
        #target line polygon
        if not QgsWkbTypes.geometryType(wkbType) in [
                QgsWkbTypes.PointGeometry, 
                QgsWkbTypes.LineGeometry, 
                QgsWkbTypes.PolygonGeometry]:
            #Not applicable geometryType
            raise ValueError("activeLayer is out of scope geometryType")

        #Target 2D layers
        if QgsWkbTypes.hasZ(wkbType) or QgsWkbTypes.hasM(wkbType):
            #Not applicable wkbType
            raise ValueError("activeLayer is out of scope wkbType")

    except ValueError as e:
        #Not applicable
        self.iface.messageBar().pushMessage("tky2jgd <b>{}</b>".format(e), Qgis.Warning)
        return False
    else:
        #Target
        return True

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

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

    inFeat = QgsFeature()
    feats = self.layer.getFeatures()
    while feats.nextFeature( inFeat ):
        attributes = inFeat.attributes()
        beforeGeom = QgsGeometry( inFeat.geometry() )
        if not self.isLonlat:
            beforeGeom.transform(self.trans4301)
        afterGeom = self.moveCorrection(beforeGeom)
        if not self.isLonlat:
            afterGeom.transform(self.transTo)

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

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

###################################
#Create a converted memory layer
###################################
def createAfterLayer(self):
    layerUri = self.wkbType +"?crs=postgis:" + str(self.toSrid)
    afterLayer = QgsVectorLayer(layerUri, "exchange", "memory")
    afterLayer.dataProvider().addAttributes(self.layer.fields())
    return afterLayer

###################################
#Move the top of a feature
###################################
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

###################################
#Generate coordinate transformation class
###################################
def setCrsTrans(self):
    crs4301 = QgsCoordinateReferenceSystem(4301, 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(crs4301, crsTo, QgsProject.instance())

I didn't start with class design properly, so it was kind of messy. This area will be cut out into classes at the end and organized, and we will give priority to moving first. It's so dull that the isScope function makes me want to die.

Verification

It is troublesome to make sample data soberly.

Let's verify that the sample data created by ourselves can be converted even at right angles to the plane.

    1. Created a line layer of Japan Geodetic System 16th EPSG: 30176 with scratch layer.
  1. Read OSM and scribble on Iriomote Island.
    1. Converted with the tky2jgd plugin.

The SRID of the converted layer is properly 2458. © OpenStreetMap contributors image.png

  1. Make points and check using Web version TKY2JGD --Geographical Survey Institute.
layer X Y
Before conversion -21165 -177520
After conversion -21313 -177066
Web version TKY2JGD -21315.4052 -177083.9327

The error in the Y direction is great, but if the calculation is wrong, this error should not be enough, so maybe it's okay?

Let's continue to try it in Tokyo. At the same time, the geometry type is verified.

    1. Create a multi-point layer this time with the Japanese geodetic system 9th EPSG: 30169 with the scratch layer.
  1. Write a multipoint near the Imperial Palace.
    1. Converted with the tky2jgd plugin.

Multipoint layers can also be converted, and the SRID of this layer is 2451. © OpenStreetMap contributors image.png

layer X Y
Before conversion -6984 -34730
After conversion -7276 -34371
Web version TKY2JGD -7277.1503 -34374.4408

I don't have a bad feeling, but I close my eyes and move on.

Now verify if the attributes are inherited Add character, integer, and floating point columns to the memory layer before adding features. Let's give an attribute value. When the attributes of the converted layer are displayed, they are retained properly. image.png

With this, all the functions that were the subject of this time have been implemented.

Click here for the source so far

  • However, there is a big bug https://github.com/ozo360/tky2jgd/releases/tag/sono2

The identity of the unpleasant premonition

Try saving the converted layer in a Shapefile. Of course, the specified CRS is that of the converted layer. When I load the saved Shapefile into QGIS again, it shifts. I have to get on the fly near the original layer, but I get on the converted location. When I check the coordinates, they are really out of alignment. A value that is similar to the value obtained by converting the old coordinates to the new coordinates and then moving them to the new coordinates.

You may have made the wrong way to add it to your project, or you may have made a fundamental mistake before that. Normally, I think that you should post after fixing the problem, but it's okay to have such trial and error. That's why I'm hibernating, so I can't update it for a while, but next time I'll investigate this bug.

Postscript 2020/01/14

After moving in the parameter file by converting latitude and longitude, it was a mistake to further convert from the old coordinates to the new coordinates with QgsCoordinateTransform, so I corrected it including "Part 1". This made the latitude / longitude conversion normal. It was found that the parameter file was moved twice as much when converting the plane orthogonality. The cause of this has not yet been found.

License for this article

![Creative Commons License](https://i.creativecommons.org/l/by/4.0 /88x31.png) This article is provided under Creative Commons Attribution 4.0 International License .

Recommended Posts

Make a tky2jgd plugin with no practicality in QGIS Part 2
Make a tky2jgd plugin with no practicality in QGIS Part 1
Make a Yes No Popup with Kivy
How to make a QGIS plugin (package generation)
Make a simple Slackbot with interactive button in python
How to make a shooting game with toio (Part 1)
Make a Blueqat backend ~ Part 1
Make a Blueqat backend ~ Part 2
Make a bookmarklet in Python
Make a fortune with Python
Make a fire with kdeplot
Let's make a WEB application for phone book with flask Part 1
Make a cat detector with Google Colabratory (Part 2) [Python] ~ Use OpenCV ~
Let's make a WEB application for phone book with flask Part 2
Make a Spinbox that can be displayed in Binary with Tkinter
Make a thermometer with Raspberry Pi and make it viewable with a browser Part 4
Let's make a WEB application for phone book with flask Part 3
Let's make a WEB application for phone book with flask Part 4
Make a Spinbox that can be displayed in HEX with Tkinter
Play with a turtle with turtle graphics (Part 1)
Let's make a GUI with python.
Make a sound with Jupyter notebook
Let's make a breakout with wxPython
Make a recommender system with python
Write a vim plugin in Python
[Blender] How to make a Blender plugin
Make a filter with a django template
Let's make a graph with python! !!
Let's make a supercomputer with xCAT
Make a model iterator with PySide
Make a nice graph with plotly
Make a curtain generator in Blender
I made a plug-in "EZPrinter" that easily outputs map PDF with QGIS.
Spiral book in Python! Python with a spiral book! (Chapter 14 ~)
Make a video player with PySimpleGUI + OpenCV
Draw a heart in Ruby with PyCall
Make a rare gacha simulator with Flask
Make a Notebook Pipeline with Kedro + Papermill
Draw a heart in Python Part 2 (SymPy)
Make a partially zoomed figure with matplotlib
Make a drawing quiz with kivy + PyTorch
Let's make a voice slowly with Python
Make a cascade classifier with google colaboratory
Let's make a simple language with PLY 1
Make a logic circuit with a perceptron (multilayer perceptron)
Make a wash-drying timer with a Raspberry Pi
Write a simple Vim Plugin in Python 3
Make a GIF animation with folder monitoring
Let's make a web framework with Python! (1)
Let's make a tic-tac-toe AI with Pylearn 2
Let's make a combination calculation in Python
Make a desktop app with Python with Electron
Let's make a Twitter Bot with Python!
Let's make a web framework with Python! (2)
Draw a graph with PyQtGraph Part 1-Drawing
Let's make a Backend plugin for Errbot
Make a thermometer with Raspberry Pi and make it visible on the browser Part 3
Let's make a nervous breakdown application with Vue.js and Django-Rest-Framework [Part 3] ~ Implementation of nervous breakdown ~
Let's make a nervous breakdown app with Vue.js and Django-Rest-Framework [Part 2] ~ Vue setup ~
Let's make a nervous breakdown app with Vue.js and Django-Rest-Framework [Part 1] ~ Django setup ~
Let's make a nervous breakdown application with Vue.js and Django-Rest-Framework [Part 6] ~ User Authentication 2 ~