# 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.

So I want to average the values from each pass....baring 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 ... and 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)? ...or is there a better way than this?

I am new to using GDAL/numpy ....my impression is that calculations using multiple rasters is best done in numpy arrays? again, is this correct.

Thank you in advance

Becky

-
GRASS seems a better tool for this. When you import a raster into GRASS it enforces a 'region extent' on each raster, filling remaining values as no-data. You can then use r.mapcalc to do the averaging. – spatialthoughts Aug 10 '12 at 15:23

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" )