请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 42600|回复: 23

[源代码] Python3.GDAL从shp文件生成mask

[复制链接]

新浪微博达人勋

发表于 2020-6-29 20:21:13 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
没有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.     从shp_name的多边形生成形如description的mask

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

  8.     Returns:
  9.         None
  10.     """
  11.     (upper_left_x, upper_left_y, pixel_width, pixel_height, rows, cols) = description
  12.     mem = gdal.GetDriverByName('MEM')
  13.     mid_ds = mem.Create('', cols, rows, 1, gdal.GDT_Byte)
  14.     mid_ds.SetGeoTransform([upper_left_x, pixel_width, 0, upper_left_y, 0, pixel_height])
  15.     mid_ds.SetMetadataItem('AREA_OR_POINT', 'Point')
  16.     mid_ds.GetRasterBand(1).WriteArray(np.ones((rows, cols), dtype=np.bool))
  17.     srs = osr.SpatialReference()
  18.     srs.SetWellKnownGeogCS('WGS84')
  19.     mid_ds.SetProjection(srs.ExportToWkt())
  20.     mask_ds = gdal.Warp('', mid_ds, format='MEM', cutlineDSName=shp_name)
  21.     out_format = mask_name[-3:]
  22.     if out_format == 'tif':
  23.         ds2GTiff(mask_ds, mask_name)
  24.     elif out_format == 'npy':
  25.         ds2npy(mask_ds, mask_name)
  26.     else:
  27.         ds2GTiff(mask_ds, mask_name+'.tif')
  28.         ds2npy(mask_ds, mask_name+'.npy')
  29.     del mid_ds, mask_ds

  30. def ds2GTiff(ds, tif_name):
  31.     gtiff = gdal.GetDriverByName('GTiff')
  32.     result = gtiff.CreateCopy(tif_name, ds)
  33.     result.FlushCache()
  34.     del result

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

  37. if __name__ == '__main__':
  38.     shp_name = r'E:\data\map\cn1984\cn1984.shp'
  39.     mask_name = r'mask.tif'
  40.     description = (70, 60, 0.05, -0.05, 1201, 1401)
  41.     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相关的那段),这样就减少了硬盘的读写,节约空间提高效率。

快,点赞夸我!

点评

楼主用心了  发表于 2020-6-29 22:29

评分

参与人数 1金钱 +10 收起 理由
暗藏 + 10 厉害啦~

查看全部评分

本帖被以下淘专辑推荐:

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2020-6-30 15:36:45 | 显示全部楼层
学习一下,谢谢分享!
不知道我这么理解对不对:您这种方法是针对成图后的tiff进行剪裁,相当于先正常画很大范围的图,然后再用shp针对图像剪裁,和set_clip_path那种对某一类(如contourf)的collection的剪裁还是有原理上的不同的。
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2020-6-30 09:37:31 | 显示全部楼层
楼主太厉害啦!这个正好需要,感谢~~{:5_213:}
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

0
早起挑战累计收入
发表于 2020-6-30 13:16:25 | 显示全部楼层
建议多来论坛分享,公众号的受众太有限,而且不利于后期搜索查找。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-6-30 13:24:12 | 显示全部楼层
mofangbao 发表于 2020-6-30 13:16
建议多来论坛分享,公众号的受众太有限,而且不利于后期搜索查找。

是的,公众号只能当传播途径,很难有讨论的氛围
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-6-30 17:23:41 | 显示全部楼层
平流层的萝卜 发表于 2020-6-30 15:36
学习一下,谢谢分享!
不知道我这么理解对不对:您这种方法是针对成图后的tiff进行剪裁,相当于先正常画很 ...

前辈您好

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

合影
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-7-15 20:53:10 | 显示全部楼层
{:5_213:}{:5_213:}大佬厉害
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-7-31 13:41:23 | 显示全部楼层
楼主,我想将最终形成的mask与我的格点数据分辨率匹配,从全球资料中提取中国区域的数据,所以将mask的分辨率和行列数修改成“description = (70.5, 60.5, 0.1, -0.1, 41, 66)”,但最后的结果都是0是怎么回事呀?希望各位大佬能帮我解答一下
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-7-31 16:47:22 | 显示全部楼层
不用回头 发表于 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
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

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

就是这个原因,我修改了一下分辨率之后就好了,谢谢楼主
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表