Python3.GDAL从shp文件生成mask

墨家大宝

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

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

复制代码

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

你以为这就完了吗？

                               
登录/注册后可看大图

I'll do you one better!

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

举个例子：

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

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

1. def 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格式，否则两种格式都生成
   
10. 
    
11.     Returns:
    
12.         None
    
13.     """
    
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. 
    
34. def ds2GTiff(ds, tif_name):
    
35.     gtiff = gdal.GetDriverByName('GTiff')
    
36.     result = gtiff.CreateCopy(tif_name, ds)
    
37.     result.FlushCache()
    
38.     del result
    
39. 
    
40. def ds2npy(ds, npy_name):
    
41.     np.save(npy_name, ds.ReadAsArray().astype(np.bool), allow_pickle=False)
    
42. 
    
43. if __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相关的那段），这样就减少了硬盘的读写，节约空间提高效率。

快，点赞夸我！

平流层的萝卜

学习一下，谢谢分享！
不知道我这么理解对不对：您这种方法是针对成图后的tiff进行剪裁，相当于先正常画很大范围的图，然后再用shp针对图像剪裁，和set_clip_path那种对某一类（如contourf）的collection的剪裁还是有原理上的不同的。

暗藏

楼主太厉害啦！这个正好需要，感谢~~{:5_213:}

mofangbao

建议多来论坛分享，公众号的受众太有限，而且不利于后期搜索查找。

墨家大宝

mofangbao 发表于 2020-6-30 13:16
建议多来论坛分享，公众号的受众太有限，而且不利于后期搜索查找。

是的，公众号只能当传播途径，很难有讨论的氛围

墨家大宝

平流层的萝卜 发表于 2020-6-30 15:36
学习一下，谢谢分享！
不知道我这么理解对不对：您这种方法是针对成图后的tiff进行剪裁，相当于先正常画很 ...

前辈您好

用我所述的方法，绘图过程中就不需要shp文件参与裁剪了，相当于用mask掩膜代替shp进行了白化。而mask就是和要绘图数据同形的二维矩阵，是在绘图之前由shp文件通过本文所述方法生成的。

合影

倔强青铜

{:5_213:}{:5_213:}大佬厉害

不用回头

楼主，我想将最终形成的mask与我的格点数据分辨率匹配，从全球资料中提取中国区域的数据，所以将mask的分辨率和行列数修改成“description = (70.5, 60.5, 0.1, -0.1, 41, 66)”，但最后的结果都是0是怎么回事呀？希望各位大佬能帮我解答一下

墨家大宝

不用回头 发表于 2020-7-31 13:41
楼主，我想将最终形成的mask与我的格点数据分辨率匹配，从全球资料中提取中国区域的数据，所以将mask的分辨 ...

你设置的参数表示生成的mask只有41行66列，左上角为70.5°E60.5°N，空间分辨率为0.1°，就是说径向跨度只有6.6°。纬向跨度只有4.1°，这个区域好像和中国领土没有交集吧，自然都是0了。需要修改一下参数，并确保地图的shp是未投影过的地理坐标系，地图参考https://blog.csdn.net/modabao/article/details/103027316

不用回头

墨家大宝 发表于 2020-7-31 16:47
你设置的参数表示生成的mask只有41行66列，左上角为70.5°E60.5°N，空间分辨率为0.1°，就是说径向跨度 ...

就是这个原因，我修改了一下分辨率之后就好了，谢谢楼主