Making CRS transformation of QgsVectorLayer with PyQGIS

Question

I need to transform QgsVectorLayers in PyQGIS from one CRS into another. By googling I found the following solution PyQGIS Developer Cookbook, which works for points only:

crsSrc = QgsCoordinateReferenceSystem(4326)    # WGS 84
crsDest = QgsCoordinateReferenceSystem(32633)  # WGS 84 / UTM zone 33N
xform = QgsCoordinateTransform(crsSrc, crsDest)
pt1 = xform.transform(QgsPoint(18, 5))

Is there a way to make a CRS transformation of an entire layer?

Answer 1

You can directly use a loop in pyqgis to apply the same transformation xform to all features in a QgsVectorLayer and create a transformed layer in the new CRS, just as you transformed the single point in the question.

The loop to create e.g. a transformed scratchLayer2 from scratchLayer1:

feats = []
for f in scratchLayer1.getFeatures():
    g = f.geometry()
    g.transform(xform)
    f.setGeometry(g)
    feats.append(f)

scratchLayer2.dataProvider().addFeatures(feats)

The code uses QgsGeometry::transform , so it applies to points, linestring as well as other feature types.

Below is a test (QGIS 3) on transforming a memory polyline layer.

crsSrc = QgsCoordinateReferenceSystem(4326)    # WGS 84
crsDest = QgsCoordinateReferenceSystem(32633)  # WGS 84 / UTM zone 33N
xform = QgsCoordinateTransform(crsSrc, crsDest)

#create test layers    
uri1 = "linestring?crs=epsg:4326&field=id:integer"
scratchLayer1 = QgsVectorLayer(uri1, "Scratch point layer1",  "memory")

feat = QgsFeature(scratchLayer1.pendingFields())
feat.setGeometry(QgsGeometry.fromPolyline([QgsPointXY(18, 5),QgsPointXY(18, 6)]))
(res, outFeats) = scratchLayer1.dataProvider().addFeatures([feat])

uri2 = "linestring?crs=epsg:32633&field=id:integer"
scratchLayer2 = QgsVectorLayer(uri2, "Scratch point layer2",  "memory")

#CRS transformation
feats = []
for f in scratchLayer1.getFeatures():
    g = f.geometry()
    g.transform(xform)
    f.setGeometry(g)
    feats.append(f)


scratchLayer2.dataProvider().addFeatures(feats)

#show results:    
for f in scratchLayer2.getFeatures():
    print(f.geometry().exportToWkt())

Output of the test:

LineString (832713.79872411652468145 553423.98688296671025455, 832157.79170280916150659 664114.16206580714788288)

Answer 2

You can use the "Reproject layer" algorithm.

QGIS 3

See the documentation for more details.

params = {
    'INPUT': inputLayer, 
    'TARGET_CRS': 'EPSG:32633',
    'OUTPUT': "path/to/layerInUtm33N.shp"
}

processing.run("native:reprojectlayer", params)

QGIS 2

See the documentation for more details.

processing.runalg('qgis:reprojectlayer', inputLayer, 'EPSG:32633', "path/to/layerInUtm33N.shp")

Answer 3

I do this in a plugin that I wrote, source code is here: https://github.com/icsm-au/icsm_qgis_transformer

The important part is this:

    dest_crs = QgsCoordinateReferenceSystem()
    dest_crs.createFromProj4(proj_string)

    temp_dir = tempfile.mkdtemp()
    temp_outfilename = os.path.join(temp_dir, 'temp_file.shp')
    error = QgsVectorFileWriter.writeAsVectorFormat(layer, temp_outfilename, 'utf-8', dest_crs, 'ESRI Shapefile')

That is doing a transform while writing out to file. If you want the file to be temporary, you could transform into a virtual layer, I think.

Edit: super quick edit, but reading comments makes me think that what you want is to do a transform like above, resulting in a virtual (in memory) file. That's how QGIS works, you can't just have it magically converting, unless you use the built in on-the-fly reprojection.

Answer 4

@MichaelMiles-Stimson suggested OGR2OGR, which is packed with QGIS. One option is to use it through command line, with the following command:

ogr2ogr -t_srs epsg:32633 new.shp old.shp

You might prefer to call this command directly from python. Elevine's answer elsewhere shows how:

1) download http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/ogr2ogr.py and store it somewhere in your python import path (/usr/lib/pythonX.X/dist-packages or /usr/local/lib/pythonX.X/dist-packages on Linux, not sure about Windows)

2) write following code:

import ogr2ogr

ogr2ogr.main(["","-f", "ESRI Shapefile", "-s_srs", "epsg:4326", "-t_srs", "epsg:32633", "new.shp", "old.shp"])

In case of bugs you might need to add full paths to the file names.

EDIT: alternative way, less bug-prone:

from os import system

#define oldPath and newPath
cmd = 'ogr2ogr -f ESRI Shapefile -s_srs epsg:4326 -t_srs epsg:32633 '+newPath+' '+oldPath
os.system(cmd)

Perhaps you can use system(cmd) instead, not sure which syntax is "better".

EDIT2: now I think the better way is to write the command and the call it, but function call() is prefered for it. See a detailed post on SO on this topic.