.npy文件_Python3.GDAL从shp文件生成mask-程序员宅基地

技术标签: input选择完文件操作  gdal不支持fgdbr格式  .npy文件  arcgis打开tif文件  js 读取mid文件  矢量数据shp七个文件介绍  

没有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

你以为这就完了吗?

eab34c893672b7434262de4e624dd536.gif

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。

11300af7d33f014787480fb0e6864112.png

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

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

cb6dd597bca76736c312b0fef6b3d5fe.png

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

快,点赞夸我!

GitHub:Mo-Dabao

微信公众号:碎积云

CSDN:墨大宝

知乎:墨大宝

气象家园:墨家大宝)

版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/weixin_39719476/article/details/110417789

智能推荐

GDKOI 2021 提高组 Day1 第三题 回文(manachar+ST表)_gdkoi2021 day1t3-程序员宅基地

文章浏览阅读175次。GDKOI 2021 提高组 Day1 第三题 回文题目大意给出长为nnn的串,和qqq组询问,每次询问区间中的最长回文串。n,q≤5∗105n,q\le5*10^5n,q≤5∗105题解可以先用manachar求出以每个位置为中心的回文串,询问时二分答案,然后在区间中判断是否存在长度为midmidmid的回文串,用ST表维护区间最值。注意二分判断时并非在整个区间[l,r][l,r][l,r]中找最大值,而需分别将左端点右移和右端点左移大约midmidmid(因奇偶而不同)的位置,以保证找_gdkoi2021 day1t3

课程设计-点餐系统_点餐系统课程设计-程序员宅基地

文章浏览阅读1.2w次,点赞14次,收藏151次。https://open.weixin.qq.com/connect/oauth2/authorize?appid=APPID&redirect_uri=REDIRECT_URI&response_type=code&scope=SCOPE&state=STATE#wechat_redirect_点餐系统课程设计

韦东山全新系列——led驱动程序框架分析_韦东山全系列-led驱动程序框架-程序员宅基地

文章浏览阅读576次。这次,韦东山老师计划同时开讲若干款板子,那么就不能一款板子,写一遍代码,因此,抽象出了一个框架,或者说模板。首先看Makefile,具体如下,和之前的相比,可以看出,这次多次使用了变量,以解决不同板子之间的差异。并且在其中,增加了测试程序的编译# 1. 使用不同的开发板内核时, 一定要修改KERN_DIR# 2. KERN_DIR中的内核要事先配置、编译, 为了能编译内核, 要先设..._韦东山全系列-led驱动程序框架

8-Arm PEG-SCM,8-Arm-PEG-Succinimidyl NHS ester-程序员宅基地

文章浏览阅读86次。英文名称:8-Arm PEG-SCM中文名称:八臂-聚乙二醇-琥珀酰亚胺乙酸酯分子量:1k,2k,3.4k,5k,10k,20k(可按需定制)质量控制:95%+存储条件:-20°C,避光,避湿用 途:仅供科研实验使用,不用于诊治外 观:粘稠液体或者固体粉末,取决于分子量注意事项:取用一定要干燥,避免频繁的溶解和冻干溶解性:溶于大部分有机溶剂,如:DCM、DMF、DMSO、THF等等。在水中有很好的溶解性8臂PEG-SCM是一种多臂PEG衍生物,在连接到一个六甘油核心的八个

Ubuntu16.04 从源码安装python3.6.x并升级pip_ubuntu 源码安装 python pip-程序员宅基地

文章浏览阅读135次。首先不要卸载系统自带的python版本,否则会有一系列麻烦!!!1.下载源码https://www.python.org/downloads/release/python-369/2.从源码安装python执行命令:tar xfz Python-3.6.9.tgzcd Python-3.6.9./configure --with-sslmakesudo make install此时pip不能用了3.更新pip执行命令:sudo python3.6 get-pip._ubuntu 源码安装 python pip

Linux学习笔记_linux历史发展-程序员宅基地

文章浏览阅读1.4k次,点赞8次,收藏40次。Linux学习笔记Linux的历史操作系统,英语Operating System简称为OS。说道操作系统就需要先讲一讲Unix,UNIX操作系统,是一个强大的多用户、多任务操作系统,支持多种处理器架构,按照操作系统的分类,属于分时操作系统,最早由KenThompson、Dennis Ritchie和Douglas McIlroy于1969年在AT&T的贝尔实验室开发。而linux就是..._linux历史发展

随便推点

解读文献(五)------基于阻抗控制_基于力的阻抗控制-程序员宅基地

文章浏览阅读6.9k次,点赞13次,收藏140次。文章目录背景主动柔顺控制阻抗控制基于力的阻抗控制基于位置的阻抗控制力位混合控制背景主动柔顺控制阻抗控制基于力的阻抗控制基于位置的阻抗控制力位混合控制..._基于力的阻抗控制

webservice 的调用 服务端 参数 List<Object> 的实例,-程序员宅基地

文章浏览阅读2k次。2019独角兽企业重金招聘Python工程师标准>>> ..._webservice参数list

Android 中this、getContext()、getApplicationContext()、getApplication()、getBaseContext() 之间的区别...-程序员宅基地

文章浏览阅读157次。 : 知之为知之,不知为不知是知也! 使用this, 说明当前类是context的子类,一般是activity application等; this:代表当前,在Activity当中就是代表当前的Activity,换句话说就是Activity.this在Activity当中可以缩写为this. Activity.this的context 返回当前..._android中getapplicationcontext和getbasecontext和this和getapplication

ajax ssm 页面跳转_SSM整合Ajax-程序员宅基地

文章浏览阅读512次。1.导入jar包: 2.配置spring-servlet.xml:text/html;charset=UTF-8text/json;charset=UTF-8application/json;charset=UTF-83.控制器中用@ResponseBody实现返回json数据格式:@ResponseBody@RequestMapping(value="queryByUser",method=Re..._ssm页面跳转ajax提交数据

pip3 mysql_python3 pip install mysqlclient报错-程序员宅基地

文章浏览阅读200次。报错显示如下,不清楚是哪里的编码问题,环境是win10+python3.7PS C:\Users\zhu\python\mysite> pip3 install mysqlclientCollecting mysqlclientUsing cached mysqlclient-1.3.12.tar.gzInstalling collected packages: mysqlclientRun..._pip3 install mysqlclient 指定源

MindSpore21天实战营(2):基于BERT实现中文新闻分类实战_bert 中文实战-程序员宅基地

文章浏览阅读1k次。转载地址:MindSpore21天实战营(2):基于BERT实现中文新闻分类实战_MindSpore_昇腾论坛_华为云论坛作者:胡琦“ModelArts + MindSpore”实战BERT中文新闻分类Copy攻城狮人狠话不多,学AI就到huaweicloud.ai,和“MM”一起玩转AI。前言大家好,我是Copy攻城狮胡琦,有幸参与华为业界首个全场景AI实战营。今天是MIndSpore 21天实战营的第二次课,光接触的名词就已经碉堡了--一站式AI开发平台ModelArts、全栈全场景_bert 中文实战

推荐文章

热门文章

相关标签