Aligning many rasters using PyQGIS or Python?

Question

I want to align many rasters in some pixel size, extent and projection system using Python or PyQGIS.

First think is to use GDAL:

1. gdalinfo (to find information from reference raster)

2. gdalwarp (to convert projection system and pixel size)

3. gdal_translate (to convert extent)

They work, but it is not easy to use those tools together in Python and they need much time to finish this work. Also, it can be done easily using QGIS and Align Rasters Tool.

Does a tool doing this work using PyQGIS or Python exist? (I want to work more programming automatically)

UPDATE

I find this code from this question:

from osgeo import gdal, gdalconst

inputfile = #Path to input file
input = gdal.Open(inputfile, gdalconst.GA_ReadOnly)
inputProj = input.GetProjection()
inputTrans = input.GetGeoTransform()

referencefile = #Path to reference file
reference = gdal.Open(referencefile, gdalconst.GAReadOnly)
referenceProj = reference.GetProjection()
referenceTrans = reference.GetGeoTransform()
bandreference = reference.GetRasterBand(1)    
x = reference.RasterXSize 
y = reference.RasterYSize

outputfile = #Path to output file
driver= gdal.GetDriverByName('GTiff')
output = driver.Create(outputfile, x, y, 1, bandreference.DataType)
output.SetGeoTransform(referenceTrans)
output.SetProjection(referenceProj)

gdal.ReprojectImage(input, output, inputProj, referenceProj, gdalconst.GRA_Bilinear)

del output

That code works fine except raster clip. Any idea how to update this code to clip input raster to extent of reference raster in code?

Answer 1

GDAL approach and time

I agree with @user:78446 - gdal is the best way forwards and you should be able to use the tools you already mention i.e. gdalwarp, gdaltranslate etc. This will be quite a long process as for each raster, each cell needs to be realigned and essentially re-interpolated to account for the cell realignment - even longer based on raster size.

Incidentally in any following analyses, make sure to account for these raster pixel value changes too - especially if the raster realignments are over very large distances.

---

GDAL and Python

Looking around the forum you'll find lots of examples demonstrating how to string gdal and python together e.g. here.

Also, check out the help here for working with gdal and Python together.

Answer 2

It seems you can use gdal.Warp instead of ReprojectImage, as in

OutTile = gdal.Warp(OutTileName, Raster,
  format=RasterFormat, outputBounds=[minX, minY, maxX, maxY], 
  xRes=PixelRes, yRes=PixelRes, dstSRS=Projection, 
  resampleAlg=gdal.GRA_NearestNeighbour, options=['COMPRESS=DEFLATE']
)

(Code is taken from https://gis.stackexchange.com/a/237372/15183).

Answer 3

Works fine for me. Only fixed to GA_ReadOnly in a line

reference = gdal.Open(referencefile, gdalconst.GAReadOnly)

to

reference = gdal.Open(referencefile, gdalconst.GA_ReadOnly)

Script reproject and clip the input raster (whole Europe) to resolution and extent of reference raster (central part), result:

gdalwarp can be also used for aligning rasters (Changing pixel size in 'asc' file using Qgis?)