.npy文件_Python3.GDAL从shp文件生成mask

没有ArcGIS的矢量转栅格工具的时候用shp多边形从栅格数据中抠出一块来？

from osgeo import gdal
result = gdal.Warp('masked.tif', 'input.tif', cutlineDSName='input.shp')
result.FlushCache()
del result

BOOM！完成！input.tif 被 input.shp 抠出来的部分就保存为了 masked.tif

你以为这就完了吗？

I'll do you one better!

有时候有没有觉得想要一个和你的栅格数据同形的掩膜mask？

举个例子：

比如要分析大量栅格数据中固定感兴趣区域的某个统计量，原来每个栅格都要扣一遍生成一个扣过的文件再分析，如果有了感兴趣区域的mask只要读取每个栅格再mask一下数据就出来了。通常mask比用分辨率很高的shp扣要快，而且这个mask可以反复用反复节约时间、硬盘空间。

我写成了函数可以直接调用：

def shp2mask(shp_name, description, mask_name='mask'):
    """从shp文件生成mask

    从shp_name的多边形生成形如description的mask

    Args:
        shp_name: shapfile文件名
        description: (upper_left_x, upper_left_y, pixel_width, pixel_height, rows, cols)
        mask_name: 生成的mask文件名，支持.tif和.npy格式，否则两种格式都生成

    Returns:
        None
    """
    (upper_left_x, upper_left_y, pixel_width, pixel_height, rows, cols) = description
    mem = gdal.GetDriverByName('MEM')
    mid_ds = mem.Create('', cols, rows, 1, gdal.GDT_Byte)
    mid_ds.SetGeoTransform([upper_left_x, pixel_width, 0, upper_left_y, 0, pixel_height])
    mid_ds.SetMetadataItem('AREA_OR_POINT', 'Point')
    mid_ds.GetRasterBand(1).WriteArray(np.ones((rows, cols), dtype=np.bool))
    srs = osr.SpatialReference()
    srs.SetWellKnownGeogCS('WGS84')
    mid_ds.SetProjection(srs.ExportToWkt())
    mask_ds = gdal.Warp('', mid_ds, format='MEM', cutlineDSName=shp_name)
    out_format = mask_name[-3:]
    if out_format == 'tif':
        ds2GTiff(mask_ds, mask_name)
    elif out_format == 'npy':
        ds2npy(mask_ds, mask_name)
    else:
        ds2GTiff(mask_ds, mask_name+'.tif')
        ds2npy(mask_ds, mask_name+'.npy')
    del mid_ds, mask_ds

def ds2GTiff(ds, tif_name):
    gtiff = gdal.GetDriverByName('GTiff')
    result = gtiff.CreateCopy(tif_name, ds)
    result.FlushCache()
    del result

def ds2npy(ds, npy_name):
    np.save(npy_name, ds.ReadAsArray().astype(np.bool), allow_pickle=False)

if __name__ == '__main__':
    shp_name = r'E:datamapcn1984cn1984.shp'
    mask_name = r'mask.tif'
    description = (70, 60, 0.05, -0.05, 1201, 1401)
    shp2mask(shp_name, description, mask_name)

调用方法可以看最后4行，传入你的shp_name，用含6个数字的列表描述想要的mask左上角的经纬度、水平垂直分辨率、行列数，指定输出mask的文件名就行，可以输出tif和npy格式的mask。

自己实现这个操作的动机是画图。前一段时间几个群里总有人问如何“白化”，当时我还不知道啥叫白化，后来自己也实现了作图过程中用set_clip_path方法白化，甚至也独立尝试出了同时白化多个区域的方法。实现后才发现气象家园有个平流层的蔬菜早就做出来了，而且我和他拼接Path的想法与方法都一样。

在文档尚不十分丰富的情况下通过自己的理解充分发挥出这些代码工具的真实威力，真的是人生一大快事，而发现前路有人更让我热血沸腾，带酒赶路，在此分享出白化的另一种方法（构造np.ma.masked_array就行了）。

回过来说说这篇教程，shp文件裁剪栅格，非常基础的操作，不知是不是我找得姿势不对，网上尤其是中文网络中并没有用Python实现的介绍，只有几篇借助PIL库和用Python调用系统命令操作GDAL的相关介绍，费了老鼻子劲儿才找到问题的核心gdal.Warp。另外值得吹一吹的是，我是读网上一篇C++的GDAL教程才找到在内存中建立gdal.Dataset的方法（就是MEM相关的那段），这样就减少了硬盘的读写，节约空间提高效率。

快，点赞夸我！

GitHub：Mo-Dabao

微信公众号：碎积云

CSDN：墨大宝

知乎：墨大宝

气象家园：墨家大宝)