Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

share|improve this question
    
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
add comment

1 Answer

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 )
share|improve this answer
    
Thank you very much for your input - I will give it a go! –  Becky Aug 14 '12 at 8:53
    
And it did work! –  Becky Aug 14 '12 at 15:17
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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