Define the extent of QgsRasterCalculator

Question

I have several raster files (.tif) in a folder that have different resolution and extent but are all in the same study area.

Before merging them, I want to bring them to the same resolution (5 m) and extent. For this, I used the Raster Calculator inside of QGIS and manually defined Xmin, Xmax, Ymin and Ymax and the number of columns and rows. This solution worked fine and I could merge the resulting files.

I want to automate this processing with the Python Plugin, but so far I am not getting the correct results, because I struggling with defining the extent in the RasterCalculator. My code looks like this:

import os

pathin =  r'C:\Documents\Raster_files'
pathout = r'C:\Documents\Processed'
filelist = os.listdir(pathin)

llc = QgsPointXY(319660,4624088)# lower left corner
urc = QgsPointXY(396425,4742843)# upper right corner
aoi = QgsRectangle(llc,urc)

for f in filelist:
    lyr1 = QgsRasterLayer(f)
    output = os.path.join(pathout,'5m_'+f)
    entries = []
    ras = QgsRasterCalculatorEntry()
    ras.ref = 'ras@1'
    ras.raster = lyr1
    ras.bandNumber = 1
    entries.append(ras)
    calc = QgsRasterCalculator('ras@1', output, 'GTiff', aoi, 15353, 23751, entries)
    calc.processCalculation()

So far, he is creating new files like expected but they are in a completetly wrong location and have values from -inf to +inf (probably because he does not use the UTM crs of the raster files). Does anybode have suggestions?

After this part is working, I intend to include gdal merging of the resulting files with very high compression.

Answer 1

Since you are not actually performing any raster calculations, I would recommend using a more appropriate tool such as gdal_translate for this. You can call the command line utility as a library function in Python (Geo tips & tricks- GDAL and OGR utilities as library functions).

Here is an example of how to use it on your raster directory to clip your rasters to your study area and adjust the resolution to 5m with a bilinear re-sampling method (You can use a different method or if you don't want to resample, simply omit that argument). *just double check that the upper-left and lower-right corner coordinates are correct.

from osgeo import gdal
import os

in_folder = 'C:\\Documents\\Raster_files'
out_folder = 'C:\\Documents\\Processed'

ulx = 319660
uly = 4742843
lrx = 396425
lry = 4624088

testdir = os.scandir(in_folder)

for file in testdir:
    if file.name.endswith('.tif'):
        name = f'{os.path.splitext(file.name)[0]}_5m.tif'
        input = file.path
        output = os.path.join(out_folder, name)
        ds = gdal.Open(input)
        ds = gdal.Translate(output, ds, projWin = [ulx, uly, lrx, lry], xRes = 5, yRes = -5, resampleAlg = 'bilinear')
        ds = None

You can see more about translateOptions here.