1

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?

2

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"])  
  • I imported the modules and copied the code, but i am unsure if the square bracket in the last line is necessary. – Lauzeruz Sep 3 at 10:01
  • it was a typo, I forgot to close it. I am adding two lists , it could also be written as ["gdalbuildvrt","-te", str(dfGeoXUL), str(dfGeoYLR), str(dfGeoXLR), str(dfGeoYUL), "-tr", xres, yres, "outputVRT.vrt","image_to_align.tif"] – radouxju Sep 3 at 10:12

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.