# Raster algebra in Python with rasters of different extents

I am trying to find out how to use GDAL & numpy modules to average a set of rasters which are Sigma0 values from satellite passes at different times within a year.

Each raster has been mosaicked from a number of smaller images using GDAL_merge. Because the orbits differ each time they pass over the area of interest the merged datasets are different extents and not regular squares.

I want to average the values from each pass bearing in mind that when theoretically stacked on top of one another, there is sometimes overlaps of one/two/three images.

I imagine the best way of me getting around this is making all rasters the same extent, using a 'no data' value for areas of the raster where there is no data and then ignoring these values during the calc of the average.

If this indeed the best way, how do I go about making them all the same extent when they are not regular (rectangle)?

I am new to using GDAL/numpy. My impression is that calculations using multiple rasters is best done in numpy arrays.

I think you can do that pretty easily with GDAL and numpy. Mind you, I think that you will need to do more complex analysis afterwards (speckle reduction etc), but in principle, you could stack all your observations into a single multiband file using eg

``````gdalbuildvrt -separate -te xmin ymin xmax ymax -input_file_list my_filenames.txt output_file.vrt
``````

`output_file.vrt` is a dataset that gdal should understand. `xmin ymin xmax ymax` are the maximum extent (and I'm assuing they all your observations share the same spatial resolution). In python, you load up the data, set a mask for your no data value (I'm assuming 0 here, but do check), and then average:

``````from osgeo import gdal
g = gdal.Open ( "output_file.vrt" )