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.

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. 1. Created a line layer of Japan Geodetic System 16th EPSG: 30176 with scratch layer.
2. Read OSM and scribble on Iriomote Island.
3. 3. Converted with the tky2jgd plugin.

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

4. Make points and check using Web version TKY2JGD --Geographical Survey Institute.

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. 1. Create a multi-point layer this time with the Japanese geodetic system 9th EPSG: 30169 with the scratch layer.
2. Write a multipoint near the Imperial Palace.
3. 3. Converted with the tky2jgd plugin.

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

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.

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

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