Determine spatial average of several tiff files

Question

I am relatively new to GIS.

I have a set of tiff files containing a raster each. I would like to average them, though their extents are slightly different. They each have their geospatial data recorded in their GeoTransform attribute. I am working in Python.

I have attempted to use GDAL to build a vrt and then average the resulting array, but I am given an error I cannot solve.

My attempt is based off this question here: raster algebra in python with rasters of different extents

#Conglomerate tiff files into a vrt
os.system('gdalbuildvrt obs_avg.vrt %s'%('Obs_tiffs/*.tif'))

print '-----Attempting to open'
g = gdal.Open('obs_avg.vrt')
data = g.ReadAsArray()
mdata = np.ma.array(data, mask=(data==0))
mean = mdata.mean(axis=0)

#Checking output
print type(mean)
print data.shape
plt.imshow(mean)
plt.show()

However, I am given the following error:

0...10...20...30...40...50...60...70...80...90...100 - done. gdalbuildvrt does not support positive NS resolution. Skipping outdir/Obs_tiffs/inputfile.tif

Which causes the vrt to not be created.

What is the cause of this error?

Answer 1

The gdalbuildvrt program isn't the right tool for this task. It's for making virtual mosaics and multiband images.

I suggest this approach:

1. Define the domain in which you want results (geotransform, rows and columns). We used to call this a "canvas" back in the day.
2. Warp (gdalwarp, nearest neighbor) your source TIFFs to that domain, producing TIFFs with identical extent and size.
3. Read and average those. You could use osgeo.gdal. I recommend using Rasterio, which is a newer Python interface to the GDAL library, to read and write. Do all your numerical analysis using Numpy arrays.