Aligning rasters with Python (QGIS)

Question

Does someone know if there is a SAGA or GDAL tool that is capable of aligning two raster layers?

For the creation of a plugin I want to automate this process so the raster alignment option of QGIS won't do. A few online sources advice to use 'SAGA Resampling' but until now I did not succeed to create a set of aligned rasters. I thought this was due to the unconventional Amersfoort RD (EPSG:28992) projection, but unaligned results are also created when using the more conventional 'EPSG:3857 - WGS 84 / Pseudo-Mercator' projection. So my question is what processing tool would you recommend to align rasters? Or is there a way to write some code which uses 'Align rasters' from QGIS itself?

Answer 1

one solution is to use gdalbuildvrt to make virtual raster aligned with one of your raster. With the -te option you can specify an extent that matches the extent of your choice, and with the -tr option you can make sure to have exactly the same resolution. The default resampling is nearest neighbour, but you can modify it with the -r option.

import subprocess, os,  sys, glob
try:
    from osgeo import gdal
    from osgeo import osr
except:
    import gdal
    import osr

command = ["gdalbuildvrt","-te"]


hDataset = gdal.Open( "image_as_template.tif", gdal.GA_ReadOnly )
adfGeoTransform = hDataset.GetGeoTransform(can_return_null = True)
if adfGeoTransform is not None:
    dfGeoXUL = adfGeoTransform[0] 
    dfGeoYUL = adfGeoTransform[3] 
    dfGeoXLR = adfGeoTransform[0] + adfGeoTransform[1] * hDataset.RasterXSize + adfGeoTransform[2] * hDataset.RasterYSize
    dfGeoYLR = adfGeoTransform[3] + adfGeoTransform[4] * hDataset.RasterXSize + adfGeoTransform[5] * hDataset.RasterYSize
    xres = str(abs(adfGeoTransform[1]))
    yres = str(abs(adfGeoTransform[5]))
    subprocess.call(command +[ str(dfGeoXUL), str(dfGeoYLR), str(dfGeoXLR), str(dfGeoYUL), "-tr", xres, yres, "outputVRT.vrt","image_to_align.tif"])