Python中使用面状矢量裁剪栅格影像,并依据Value值更改矢量属性
内容导读
互联网集市收集整理的这篇技术教程文章主要介绍了Python中使用面状矢量裁剪栅格影像,并依据Value值更改矢量属性,小编现在分享给大家,供广大互联网技能从业者学习和参考。文章包含3019字,纯文字阅读大概需要5分钟。
内容图文
![Python中使用面状矢量裁剪栅格影像,并依据Value值更改矢量属性](/upload/InfoBanner/zyjiaocheng/1331/62ad025e8d4748298413d0de9fc68250.jpg)
本文整体思路:在Python中使用Geopandas库,依次读取shp文件的每一个面状要素,获取其空间边界信息并裁剪对应的栅格影像,计算所裁剪影像Value值的众数,将其设置为对应面状要素的NewTYPE值,所有要素属性值都改好之后保存为新的shp文件。
使用Python处理空间数据确实用的不多,所以一个星期以来一直深受这个程序的折磨,官方文档、博客、谷歌、百度、论文,能用的方法都给用了,但是进度还是很慢,特别是当看到这篇博客的时候。。。好气啊。。
不过幸亏头比较铁,虽败不馁,慢慢一步一步调试找问题,一个一个解决,终于啃下了这根硬骨头。(对写程序来说,调试真的是第一法宝啊!)
好吧进入正题。。。
Pandas是一款高性能的Python数据分析库,而Geopandas是由Shapely、Fiona、PyProj、matplotlib以及其他必须的库共同构建的Pandas地理空间扩展,它既可以处理空间数据、也可以处理属性数据。看到有些博客说在读取shp文件的时候用Geopandas库,而在编辑、导出的时候用pyshp比较好,其实不然,Geopandas也提供了功能完备的导出接口,而且使用特别方便,只不过需要注意一个小的细节问题,否则就会报错。Rasterio是基于GDAL库二次封装的主要用于空间栅格数据处理的Python库,本程序需要对栅格影像进行裁剪,因此也需要引用这个库。两个库的官方文档如下,参考的时候要注意版本问题,不同的版本有些接口可能已经改变。Geopandas参考文档;Rasterio参考文档。
我是Windows 10系统,在Python中安装Geopandas库比较麻烦,不能用pip命令直接安装,而需要先下载Anaconda再用conda命令安装,这部分网上有很多参考资料,就不多赘述了。但是安装完成之后发现它的一些依赖库不能使用,需要pip命令将其卸载之后,再通过此地址:python依赖库下载对应的依赖库并安装就可以使用了。
本程序完整的代码如下:
# -*- coding: utf-8 -*- from geopandas import *; import rasterio as rio; import rasterio.mask; # 读入矢量和栅格文件 shpdatafile=‘D:/PythonConda/Data/ShpData.shp‘ rasterfile=‘D:/PythonConda/Data/Raster.tif‘ out_file=‘D:/PythonConda/Data/OutShpData‘ shpdata=GeoDataFrame.from_file(shpdatafile) rasterdata=rio.open(rasterfile) out_shpdata = shpdata.copy() #投影变换,使矢量数据与栅格数据投影参数一致 shpdata=shpdata.to_crs(rasterdata.crs) for i in range(0, len(shpdata)): # 获取矢量数据的features geo = shpdata.geometry[i] feature = [geo.__geo_interface__] #通过feature裁剪栅格影像 out_image, out_transform = rio.mask.mask(rasterdata, feature, all_touched=True, crop=True, nodata=rasterdata.nodata) #获取影像Value值,并转化为list out_list = out_image.data.tolist() #除去list中的Nodata值 out_list = out_list[0] out_data = [] for k in range(len(out_list)): for j in range(len(out_list[k])): if out_list[k][j] >=0: out_data.append(out_list[k][j]) #求数据中的众数 if len(out_data): counts = np.bincount(out_data) new_type = np.argmax(counts) else: new_type = None #依据众数值更改feature的NewTYPE属性值 out_shpdata.NewTYPE[i] = new_type #将属性更改过的GeodataFrame导出为shp文件 out_shpdata.to_file(out_file)
需要注意的问题:
1)两个文件需要在同一坐标系下,需要将用于裁剪的shp数据进行投影变换,将其投影参数变为栅格数据的投影参数,由以下代码实现:
shpdata=shpdata.to_crs(rasterdata.crs)
其中crs表示数据的投影参数,to_crs为投影参数转换函数。
2)裁剪函数rasterio.mask.mask的参数问题,需要传入的矢量数据为GeoJSON数据,因此读入的每一个面状feature都需要用__geo_interface__函数进行格式转换,这一步可以通过调试来查看具体的数据格式是否正确;此函数有两个返回值,第一个记录裁剪栅格的数据值,第二个记录其坐标转换信息,一般用到第一个返回值比较多;裁剪得到的栅格形状其实是一个矩形,与feature的外接矩形区域一致,只是位于feature外部的像素值默认被设置为Nodata,当然这可以通过传入的参数进行设置。
3)获取到裁剪的栅格后,通过.data来获取其Value值,但此时还不能直接用于统计分析,需要将数据转化为List,函数如下:
out_list = out_image.data.tolist()
此时还需要调试来进一步确定out_list的数据内容,发现out_list[0]才是我们真正能用到的数据值。
4)获取众数的时候需要清除Nodata值的影响,因此用for循环把out_list中的非Nodata数据再组成一个新的List,用numpy的自带函数求其众数。因为所有的编辑并不能对shpdata源数据进行改变,所以需要构建一个shpdata的copy即out_shpdata,将求得的众数赋给out_shpdata的NewTYPE,编辑完成之后再将out_shpdata导出为完整的shp文件。
5)GeoDataFrame.to_file函数可以将out_shpdata直接导出为shp文件,仅需要传入一个路径参数即可,但需要注意由于shp文件包含.shp、.shx、.dbf和.prj,因此路径只能是一个文件夹,而不能具体到.shp。如下代码所示:
out_file=‘D:/PythonConda/Data/OutShpData‘
最后将生成的数据在Arcmap中打开,设置显示的Labels后可以看到效果如下:
至此全部完成!
原文:https://www.cnblogs.com/MatthewHome/p/10057770.html
内容总结
以上是互联网集市为您收集整理的Python中使用面状矢量裁剪栅格影像,并依据Value值更改矢量属性全部内容,希望文章能够帮你解决Python中使用面状矢量裁剪栅格影像,并依据Value值更改矢量属性所遇到的程序开发问题。 如果觉得互联网集市技术教程内容还不错,欢迎将互联网集市网站推荐给程序员好友。
内容备注
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 gblab@vip.qq.com 举报,一经查实,本站将立刻删除。
内容手机端
扫描二维码推送至手机访问。