6

I need to transform QgsVectorLayers in python from one crs into another crs. By googling I found the following solution (https://qgis.readthedocs.org/en/latest/docs/pyqgis_developer_cookbook/05_crs.html), 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 3 '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 4 '14 at 16:21
  • 2
    If you have QGIS you also have OGR2OGR, it's in the .\bin folder. – Michael Stimson Aug 4 '14 at 21:43
  • 2
    @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 3 '15 at 8:24
  • 2
    Go right ahead @PavelV. , just remember to credit the comments (out of courtesy) or others may accuse you of plagiarism. – Michael Stimson Aug 4 '15 at 3:59
4

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 4 '17 at 8:09
2
+250

You can use the reproject layer algorithm of QGIS. See documentation.

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

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.

0

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.