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.

enter image description here

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?

• Jon: I disagree that "it is not easy to use [gdal] tools together in Python." It is quite simple. Anyway, you could look at this post for ideas on how to clip your rasters with a raster. (I'd go the gdal route, myself.) gis.stackexchange.com/questions/125202/…
• Mar: @Jon that not work for me
• Jon: Well, if you want help you should probably post what you tried and what didn't work.

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:

enter image description here

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