7

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.

Improve this question
2
  • Why not use gdal directly. This is what QGIS does for raster processing.
    – radouxju
    Commented Apr 11, 2018 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
    Commented Apr 12, 2018 at 6:05

2 Answers 2

Reset to default
6
+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.

Improve this answer
5
  • 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
    Commented Jan 9, 2017 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
    Commented Jan 14, 2017 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
    Commented Jan 19, 2017 at 21:27
  • When using writeRaster in 3.10, the result looks fine but is not projected properly for me (using the UI Save As... results in a good projection however). Also note since 3.8 writeRaster is asking for 'transformContext'
    – cefect
    Commented Dec 14, 2020 at 19:11
  • Is there anyway of making the os.system(cmd) quiet? i.e. avoiding the cmd window to open and just running the function in the background?
    – Juan Ossa
    Commented Sep 27, 2021 at 7:09
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.

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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