TransWikia.com

Using Python to convert rasters (.asc) into rendered image GeoTIFFs (.tif) with new projection

Geographic Information Systems Asked by Josh Terry on July 2, 2021

I am trying to use the Python Console to essentially perform the same function as right clicking the layer, then going Export > Save As, and setting the following:

"Save Raster As..." settings to apply

The raster layer is an ascii with a CRS of EPSG:28356 that I have applied styles to. I’d like to reproject it to lon/lat (EPSG:4326), and save the styling (hence needing to convert to rendered image rather than raw data). The reason I am trying to do this is because I have approximately 100 rasters (all in exactly the same location) to convert and I’d like to develop this code further to batch process all layers in my QGIS project.

I looked at a few solutions from previous questions (Solution1) (Solution2) and pieced together their code in order to test a single layer (before making a batch):

layer = iface.activeLayer()
renderer = layer.renderer()
provider=layer.dataProvider()
pipe = QgsRasterPipe()
pipe.set(provider.clone())
pipe.set(renderer.clone())
width, height = layer.width(), layer.height()
extent = layer.extent()
crsDest = QgsCoordinateReferenceSystem(4326) #epsg of target crs (WGS84)
#(extra code here)
file_writer = QgsRasterFileWriter('D:/Output/outputName.tif')
file_writer.writeRaster(pipe, width, height, extent, crsDest)

I believe the problem I am having is that the extents of the original raster are in metres, while the new CRS expects degrees. I am not sure how to update the extents to suit the newly created TIFF. If I replace "crsDest" in the final line with "layer.extent()", a TIFF is successfully made, however it is in the original CRS (EPSG:28356).

In the line "#(extra code here)", I have tried two methods. For Method 1 I adopted code from another solution (Solution3). To get extents, I manually exported/reprojected the original raster with "right click > Export > Save As" (shown in the picture above) and copied the newly created extents to produce:

layer = iface.activeLayer()
renderer = layer.renderer()
provider=layer.dataProvider()
pipe = QgsRasterPipe()
pipe.set(provider.clone())
pipe.set(renderer.clone())
width, height = layer.width(), layer.height()
crsDest = QgsCoordinateReferenceSystem(4326)

xmin = 151.1754921309999986
ymin = -33.7468109630000015
xmax = 151.2245917170000098
ymax = -33.7199599579999969
bb = QgsRectangle(xmin, ymin, xmax, ymax)

file_writer = QgsRasterFileWriter('D:/Output/method1.tif')
file_writer.writeRaster(pipe, width, height, bb, crsDest) 

This code produces a TIFF with the correct extent and CRS is produced, but is blank – it’s black when viewing in an image viewing application and will load QGIS as a layer but with nothing in it. Using the info tool on the area where data should be gives this:

Info tool results for method 1 output raster

For Method 2, I used this solution Solution4 to produce:

import os
layer = iface.activeLayer()
renderer = layer.renderer()
provider=layer.dataProvider()
pipe = QgsRasterPipe()
pipe.set(provider.clone())
pipe.set(renderer.clone())
width, height = layer.width(), layer.height()
crsDest = QgsCoordinateReferenceSystem(4326)

crsSrc = QgsCoordinateReferenceSystem(layer.crs())
xform = QgsCoordinateTransform(crsSrc, crsDest, QgsProject.instance())
projected_extent = xform.transformBoundingBox(provider.extent())

file_writer = QgsRasterFileWriter('D:/Output/method2.tif')
file_writer.writeRaster(pipe, width, height, projected_extent, crsDest)

The result is the same as Method 1. A comparison between the produced TIFFs is below:

Comparison of QGIS "Information from provider" for manually exported TIFF vs TIFFs created with Python code

Windows Explorer view of the three TIFFs produced

How can I make either of Method 1 or 2 work, or is there a better way to approach this problem?

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP