将 GeoTIFF 图像划分为规则网格并计算平均值

Question

如何将 GeoTIFF（单波段）图像（例如，1000 x 1000 像素）子集为常规方形网格（例如，50 x 50 像素；总共 400 个网格对象）；并计算每个对象的平均纬度、平均经度和平均 DN 值？任何想法都将受到高度赞赏。

输入数据只是图像，例如从以下位置下载： http://eoimages.gsfc.nasa.gov/images/imagerecords/57000/57752/land_shallow_topo_2048.tif

这可以使用gdal打开到python中，如下所示：

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

执行此操作的一种快速便捷的方法是生成一个积分图像（AKA 总面积表）。然后，您可以使用四次查找来计算任何矩形块中的平均值。

如果您在问题中提供更多详细信息（例如，您的数据作为 Python 代码的格式是什么？）我可以提供更具体的建议。

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

Answer 2

有一个 Numpy 重塑技巧可以做到这一点。您可以将原始 2D 数组重新整形为 4D 数组，其中将落在单个 50*50 网格中的所有单元格放入两个唯一维度。调用这些轴上的任何函数将为您提供聚合结果。

让我们创建一个示例二维数组：

n = 1000
grid_size = 20

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

然后计算聚合因子和两个维度的网格数：

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

然后您可以将原始 grid 数组重新整形为 4 维，并在第二维和第四维上应用 mean。

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

您可以通过自己切片一些网格并对其应用平均值来测试它：

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

结果和差异可视化：