raster algebra in python with rasters of different extents

Question

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

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 )