Divide a GeoTIFF image into regular grids and compute mean values

Question

How can I subset a GeoTIFF (single band) image (say, 1000 by 1000 pixels) into regular square grids (say, 50 by 50 pixels; 400 grid objects in total); and calculate mean latitude, mean longitude, and mean DN value for each object? Any ideas would be highly appreciated.

The input data is just the imagery,for example as downloaded from: http://eoimages.gsfc.nasa.gov/images/imagerecords/57000/57752/land_shallow_topo_2048.tif

This can be opened into python using gdal as folows:

import gdal
geotiff = gdal.Open ('land_shallow_topo_2048.tif')
colum_numbers,row_numbers,band_numbers=geotiff.RasterXSize,geotiff.RasterYSize,geotiff.RasterCount
print (colum_numbers,row_numbers,band_numbers)
2048 1024 3

Answer 1

There is a Numpy reshaping trick to do this. You can reshape your original 2D array to a 4D array where all cells which would fall in a single 50*50 grid are put into two unique dimensions. Calling any function on these axis will give you the aggregated result.

Lets make a sample 2D array:

n = 1000
grid_size = 20

grid = np.arange(n*n).reshape(n,n)

Then calculate the aggregation factor and the number of grids in both dimensions:

factor = n / grid_size
yblocks, xblocks = np.array(grid.shape) / factor

You can then reshape the original grid array to 4 dimensions and apply the mean on the second and fourth dimension.

grid_small = grid.reshape(yblocks, factor, xblocks, factor).mean(axis=3).mean(axis=1)

You can test it by slicing a few of the grids yourself and applying the mean on them with:

assert(grid[:factor,:factor].mean() == grid_small[0,0])
assert(grid[-factor:,-factor:].mean() == grid_small[-1,-1])

The results and difference visualized:

Answer 2

A fast and convenient way to do this is to generate an integral image (AKA summed area table). You can then compute the mean in any rectangular patch with four look-ups.

If you give more details in your question (for example, what is the format of your data as Python code?) I can give more specific advice.

http://en.wikipedia.org/wiki/Summed_area_table