arcgis裁剪多个shp文件_Python3.GDAL从shp文件生成mask

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

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

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

你以为这就完了吗？

I'll do you one better!

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

举个例子：

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

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

 1def shp2mask(shp_name, description, mask_name='mask'):
 2    """从shp文件生成mask 3 4    从shp_name的多边形生成形如description的mask 5 6    Args: 7        shp_name: shapfile文件名 8        description: (upper_left_x, upper_left_y, pixel_width, pixel_height, rows, cols) 9        mask_name: 生成的mask文件名，支持.tif和.npy格式，否则两种格式都生成1011    Returns:12        None13    """
14    (upper_left_x, upper_left_y, pixel_width, pixel_height, rows, cols) = description
15    mem = gdal.GetDriverByName('MEM')
16    mid_ds = mem.Create('', cols, rows, 1, gdal.GDT_Byte)
17    mid_ds.SetGeoTransform([upper_left_x, pixel_width, 0, upper_left_y, 0, pixel_height])
18    mid_ds.SetMetadataItem('AREA_OR_POINT', 'Point')
19    mid_ds.GetRasterBand(1).WriteArray(np.ones((rows, cols), dtype=np.bool))
20    srs = osr.SpatialReference()
21    srs.SetWellKnownGeogCS('WGS84')
22    mid_ds.SetProjection(srs.ExportToWkt())
23    mask_ds = gdal.Warp('', mid_ds, format='MEM', cutlineDSName=shp_name)
24    out_format = mask_name[-3:]
25    if out_format == 'tif':
26        ds2GTiff(mask_ds, mask_name)
27    elif out_format == 'npy':
28        ds2npy(mask_ds, mask_name)
29    else:
30        ds2GTiff(mask_ds, mask_name+'.tif')
31        ds2npy(mask_ds, mask_name+'.npy')
32    del mid_ds, mask_ds
33
34def ds2GTiff(ds, tif_name):
35    gtiff = gdal.GetDriverByName('GTiff')
36    result = gtiff.CreateCopy(tif_name, ds)
37    result.FlushCache()
38    del result
39
40def ds2npy(ds, npy_name):
41    np.save(npy_name, ds.ReadAsArray().astype(np.bool), allow_pickle=False)
42
43if __name__ == '__main__':
44    shp_name = r'E:\data\map\cn1984\cn1984.shp'
45    mask_name = r'mask.tif'
46    description = (70, 60, 0.05, -0.05, 1201, 1401)
47    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 相关的那段)，这样就减少了硬盘的读写，节约空间提高效率。

快，点赞夸我！

CSDN：	墨大宝	blog.csdn.net/modabao
GitHub：	Mo-Dabao	blog.csdn.net/modabao
知乎：	墨大宝	www.zhihu.com/people/mo_dabao
气象家园：	墨家大宝	bbs.06climate.com/?86416