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


  • 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. Aug 10, 2012 at 15:23

1 Answer 1


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" )
data = g.ReadAsArray()
mdata = np.ma.array ( data, mask=( data == 0 ) )
mean_s0 = mdata.mean ( axis=0 )
  • Thank you very much for your input - I will give it a go!
    – Becky
    Aug 14, 2012 at 8:53
  • And it did work!
    – Becky
    Aug 14, 2012 at 15:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.