5

I'm writing a QGIS plugin that will reproject a raster.

I have all the pieces and it's writing a new raster with the new projection, but isn't reprojecting the underlying data.

The QGIS documentation says 'CRS to reproject to' but I'm concerned that no reprojection is being done. I'd prefer not to use GDAL, but I'm now at the point where I've exhausted my options in the PyQGIS API, I think...

My code is below:

layer.setCrs(orig_crs)
tr = QgsCoordinateTransform(orig_crs, dest_crs)

file_writer = QgsRasterFileWriter(out_file)
pipe = QgsRasterPipe()
provider = layer.dataProvider()
pipe.set(provider.clone())
result = file_writer.writeRaster(pipe,
                                 provider.xSize(),
                                 provider.ySize(),
                                 tr.transform(provider.extent()),
                                 dest_crs)

Edit: I'm editing now and will add a bounty.

The question really is this: How can you use the built-in QGIS Python API functionality to transform a raster. Answers accepted for QGIS 2.0 or 3.0.

  • Why not use gdal directly. This is what QGIS does for raster processing. – radouxju Apr 11 '18 at 12:33
  • Hey @radouxju, I did do this and it works fine. I think the question remains, though, as it ought to be possible using QGIS. – Alex Leith Apr 12 '18 at 6:05
4
+50

The QGIS documentation for QgsCoordinateTransform also says that 'transform' method is for QgsPoint or QgsRectangle objects; not for raster objects. However, I was able to get a result with your code but, it was totally weird probably for the former reason. My adapted code is:

layer = iface.activeLayer()
orig_crs = layer.crs()

dest_crs = QgsCoordinateReferenceSystem(32142)

tr = QgsCoordinateTransform(orig_crs, dest_crs)

out_file = '/home/zeito/pyqgis_data/reprojected_raster.tif'

file_writer = QgsRasterFileWriter(out_file)
pipe = QgsRasterPipe()
provider = layer.dataProvider()
pipe.set(provider.clone())

result = file_writer.writeRaster(pipe,
                                 provider.xSize(),
                                 provider.xSize(),
                                 tr.transform(provider.extent()),
                                 dest_crs)

and I got next result:

enter image description here

Original raster is (791 rows x 1680 columns):

enter image description here

and reprojected raster by using gdalwarp command in next code:

import os

layer = iface.activeLayer()

provider = layer.dataProvider()

source_path = " -of GTiff " + provider.dataSourceUri()
target_path = "/home/zeito/pyqgis_data/utah_dem_UthNorth2.tif"

CRS = layer.crs()
source_epsg = CRS.postgisSrid()
new_espg = 32142
target_epsg = " -t_srs EPSG:" + str(new_espg)

cmd = "gdalwarp -overwrite -s_srs EPSG:" + \
                        str(source_epsg) + \
                        target_epsg + \
                        source_path + \
                        " " + \
                        target_path

print cmd

os.system(cmd)

it was (practically original indistinguishable):

enter image description here

However, this raster has 800 rows x 1684 columns. So, if some possibility of reprojection exists with 'writeRaster' method, you shouldn't force it to have the same rows and columns that original raster.

  • Thanks for the response. I might be complicating things by including the transform command, but I was testing with transforming the extents to see if that made a difference (it doesn't). I think that the writeraster method is NOT transforming, just defining the projection as the new projection, which is bad. – Alex Leith Jan 9 '17 at 0:04
  • 1
    looking at the "glitchy" version, I see that provider.xSize() is used for both axes, which might not work if the x and y cell sizes differ. Is that a typo? That might explain the weird results. – Steven Kay Jan 14 '17 at 1:13
  • Hi @StevenKay. The cells of the raster don't get shifted, so I doubt that's the cause (my test data was a simple and small, 200 m change). – Alex Leith Jan 19 '17 at 21:27
0

You could try using QgsRasterProjector (see doc) which has a setCrs() method to "set source and destination CRS".
The QgsRasterProjector is accessible by applying the projector() method on a QgsRasterPipe object.

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.