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?

Are you wanting to make a new layer (shapefile etc.) with the new coordinate system? — Michael Stimson, Aug 03 '14 at 22:06
no, i have a layer read from a shapefile, whereby the crs of the shapefile differs from the crs of the project. Now I want to transform the geodate (crs shapefile to crs project) just after reading so that I do not have to concern it later on (eg. geom.exportToWkt()) — gallay, Aug 04 '14 at 16:21
Why not just project the whole layer using OGR2OGR, that way it's already in the right coordinate system, no need for fancy code. I've not seen anything in any API to project a whole layer, it's done a single geometry at a time - the same way OGR2OGR does it. — Michael Stimson, Aug 04 '14 at 21:37
ok, good idea, however, at the office, i am not allowed to install any software :( thus, ogr2ogr is not an option — gallay, Aug 04 '14 at 21:42
If you have QGIS you also have OGR2OGR, it's in the .\bin folder. — Michael Stimson, Aug 04 '14 at 21:43
@MichaelMiles-Stimson: this answer seems answered by your comments, but it doesn't have a proper answer. Would you like to expand the comments into an answer, or should I do it myself? — Pavel V., Aug 03 '15 at 08:24
Go right ahead @PavelV. , just remember to credit the comments (out of courtesy) or others may accuse you of plagiarism. — Michael Stimson, Aug 04 '15 at 03:59

score 6 · Answer 1 · edited Dec 23 '23 at 11:06

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

score 6 · Answer 2 · answered Sep 03 '17 at 18:50

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)

Thank you! That is the straightforward solution that I was hoping to avoid. Vincent_v_E's answer uses a builtin processing module which seems to be as clean and perfect as it can be. — bugmenot123, Sep 04 '17 at 08:09

score 2 · Answer 3 · answered Aug 30 '17 at 05:10

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.

score 0 · Answer 4 · edited May 23 '17 at 12:39

@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.

Making CRS transformation of QgsVectorLayer with PyQGIS

4 Answers4

QGIS 3

QGIS 2

Linked

Related