QGIS Python Commands – Reprojecting Layer in QGIS with Python Commands

pyqgispythonqgis

I would like to reproject a layer from one coordinate system to another in qgis using python. The first answer to this post provides something similar to what I am trying to achieve, except that in my case:

  • I am using qgis 2.18
  • I would like to use actual layers instead of “memory` layers and save the outputs to disk.

Here is what I wrote so far:

wkt = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]'
in_crs = QgsCoordinateReferenceSystem(wkt)
wkt = 'PROJCS["WGS 84 / UTM zone 34N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",21],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32634"]]'
out_crs = QgsCoordinateReferenceSystem(wkt)

points = iface.activeLayer()
in_features = points.getFeatures()
xform = QgsCoordinateTransform(in_crs, out_crs)

out_layer = QgsVectorLayer('Point', 'points_out', 'memory')
out_features = []

#reproject
for feature in in_features:
    geometry = feature.geometry()
    geometry.transform(xform)
    feature.setGeometry(geometry)
    out_features.append(feature)

out_layer.dataProvider().addFeatures(out_features)

for f in points.getFeatures():
    print f.geometry().exportToWkt()

for f in out_layer.getFeatures():
    print f.geometry().exportToWkt()

#write out_layer
error = QgsVectorFileWriter.writeAsVectorFormat(out_layer, r'D:\11202750-002 - risicoanalyse en adaptatie experiment\Basis\point_new_crs.shp', 'CP1250', None, 'ESRI Shapefile')
iface.addVectorLayer(r'D:\11202750-002 - risicoanalyse en adaptatie experiment\Basis\point_new_crs.shp', 'point_new_crs', 'ogr')
print 'write file error = ' + str(error)

Which creates the following output:

Point (20.09595631708375407 39.81595508370384806)
Point (20.09814206032278605 39.81498734462326894)
Point (20.10044513203965622 39.81704382842367806)
Point (422623.74030901060905308 4407721.2567776320502162)
Point (422809.73621926747728139 4407611.96021181158721447)
Point (423009.15451480081537738 4407838.22327789757400751)
write file error = 0

Looking only at the console output it appears that the layer has indeed been reprojected. However, the layer that gets loaded by the last line is still in the original in_crs coordinate system:

enter image description here

And as it loads a warning is shown by the GUI:
enter image description here

How can I get the written file to be in the right CRS?

EDIT

I have also noticed that the "reprojected" point layer is placed at the right coordinates but in the wrong projection. That is the points are not in the same location. It seems to simply draw the points without taking into consideration that they are in a different coordinate system.

I have also tried replacing None in the last line by an actual QgsCoordinateReferenceSystem object as indicated by the API docs

error = QgsVectorFileWriter.writeAsVectorFormat(out_layer, r'D:\11202750-002 - risicoanalyse en adaptatie experiment\Basis\point_new_crs.shp', 'CP1250', out_crs, 'ESRI Shapefile')

but this throws a write_file_error = 6 and the same GUI warning as shown above. I'm still not sure why I can't set the coordinate system of my output file

Best Answer

The simplest way to do it is to use the processing tools :

import processing

processing.runalg("qgis:reprojectlayer",
                   "path_to_your_layer",
                   "the_coordinate_system_of the reprojection",
                   "path_where_the_reprojected_layer_will_be_saved")
Related Question