如何在Python中基于多边形对栅格文件进行区域统计？
我有一个栅格文件和一个多边形的形状文件。我想要计算这个多边形覆盖区域的栅格文件的平均值。我想用一个独立的Python脚本来完成这个，所以QGIS和Starspan都不适用。而且Arcpy也不能用。我想用GDAL。你能推荐什么Python包或者方法来实现这个吗？
2 个回答
0
虽然这个问题已经过了很久，但现在有一个专门的Python包可以解决这个问题，叫做 xagg，你可以通过以下命令安装它：conda install xagg。
还有一个选择是 regionmask，它可以提供多边形的掩膜（可以在使用 xarray 的 .where().mean() 调用中使用），但它并不会根据网格单元和多边形之间的重叠来计算面积平均，所以当栅格数据分辨率较低时，可能会不太准确。
2
有一个叫做 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