2

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

2

I think I know what I was doing wrong. It looks like the coordinates carried by a layer and the coordinate reference system carried by a layer are two different things. In the script above I was only transforming the coordinates to the new CRS but not the actual CRS itself. This can only be done upon instantiation of the QgsVectorLayer object. Therefore the script now looks like this:

in_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(in_wkt)
out_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(out_wkt)

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

out_layer = QgsVectorLayer('Point?crs='+out_wkt, '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 does not generate any of the warnings shown above, and which causes the points to plot at exactly the same location, although they have different CRS's and coordinates.

To summarize, I was facing a problem in my question because I was plotting at the same coordinate system but with different coordinates.... (since qgis was automatically assigning the default CRS to the incoming layer)

1

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")
  • 1
    I am trying to gain a more fundamental understanding of the native QGIS API, and I have intentionally avoided using processing.runalg methods.... But I might resort to this solution in the future. Thanks anyway! – user32882 Sep 14 '18 at 14:35

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.