python - 如何在 Python 中基于 polyon 对栅格文件进行区域统计？

Question

我有一个光栅文件和一个多边形 shapefile。我想获得多边形覆盖区域的栅格文件的平均值。我喜欢在 python 独立脚本中执行此操作。所以 QGIS 和 Starspan 不起作用。Arcpy 也不可用。我喜欢使用GDAL。什么 Python 包？你能推荐这样做的方法吗？

Answer 1

有一个 gdal.RasterizeLayer 函数可以让您栅格化矢量图层。它有一些缺点，您需要有一个栅格化的输出数据集。此外，如果您有重叠的几何图形，您希望首先在单独的矢量图层上隔离每个几何图形，这意味着您必须遍历所有几何图形。

使用 gdal，您可以使用 MEM 驱动程序创建内存中的文件，这使它更容易一些，但仍然存在大量数据集创建开销。

对于每个几何图形，步骤或多或少类似于：

driver = gdal.GetDriverByName('MEM')
outds = driver.Create('', pixelxsize, pixelysize, 1, GDT_Byte)
outds.SetProjection(target_proj)
outds.SetGeoTransform(target_gt)

gdal.RasterizeLayer(outds, [1], vectorlayer, burn_values=[1])

现在 outds 包含几何图形的掩码，例如将其与 np.masked_where 一起使用，您可以隔离几何图形中的像素。

它不如它好，但一旦你有一个多边形的掩码数组，它很容易通过使用 numpy/scipy 获得统计数据。

编辑：有关更详细的示例，请参阅此脚本：http: //svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py