python矢量裁剪栅格代码_Python中使用面状矢量裁剪栅格影像,并依据Value值更改矢量属性...
本文整体思路:在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依赖库下载对应的依赖库并安装就可以使用了。
本程序完整的代码如下:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
# -*- coding: utf-8 -*-
from geopandasimport *;
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 iin 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 kin range(len(out_list)):
for jin 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数据进行投影变换,将其投影参数变为栅格数据的投影参数,由以下代码实现:
1
shpdata=shpdata.to_crs(rasterdata.crs)
其中crs表示数据的投影参数,to_crs为投影参数转换函数。
2)裁剪函数rasterio.mask.mask的参数问题,需要传入的矢量数据为GeoJSON数据,因此读入的每一个面状feature都需要用__geo_interface__函数进行格式转换,这一步可以通过调试来查看具体的数据格式是否正确;此函数有两个返回值,第一个记录裁剪栅格的数据值,第二个记录其坐标转换信息,一般用到第一个返回值比较多;裁剪得到的栅格形状其实是一个矩形,与feature的外接矩形区域一致,只是位于feature外部的像素值默认被设置为Nodata,当然这可以通过传入的参数进行设置。
3)获取到裁剪的栅格后,通过.data来获取其Value值,但此时还不能直接用于统计分析,需要将数据转化为List,函数如下:
1
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。如下代码所示:
1
out_file='D:/PythonConda/Data/OutShpData'
最后将生成的数据在Arcmap中打开,设置显示的Labels后可以看到效果如下:
至此全部完成!
转载自:https://blog.csdn.net/qq_42564846/article/details/84768242
python矢量裁剪栅格代码_Python中使用面状矢量裁剪栅格影像,并依据Value值更改矢量属性...相关推荐
- python画立体图形代码_python中如何画三维的图形?
之前教了大家不少简单的画图方法,有的小伙伴已经想尝试更漂亮的高难度图形.小编找到了一个简单容易上手的分享给大家. 如果需要绘制三维的图形,matlab可以,python同样可以,代码风格也是极其的与m ...
- python控制台清屏代码_python中代码清屏的方法是什么
python中代码清屏的方法是什么 发布时间:2020-09-03 10:52:03 来源:亿速云 阅读:97 作者:小新 这篇文章给大家分享的是有关python中代码清屏的方法是什么的内容.小编觉得 ...
- python播放视频代码_Python中操作各种多媒体,视频、音频到图片的代码详解
我们经常会遇到一些对于多媒体文件修改的操作,像是对视频文件的操作:视频剪辑.字幕编辑.分离音频.视频音频混流等.又比如对音频文件的操作:音频剪辑,音频格式转换.再比如我们最常用的图片文件,格式转换.各 ...
- python添加横线代码_Python中的五种下划线
1 1.单前导下划线 _var2 3 单个下划线是一个Python命名约定,表示这个名称是供内部使用的. 它通常不由Python解释器强制执行,仅仅作为一种对程序员的提示4 5 程序员使用名称前的单下 ...
- python pandas excel 修改列_Python中Pandas读取修改excel操作攻略(代码示例)
本篇文章给大家带来的内容是关于Python中Pandas读取修改excel操作攻略(代码示例),有一定的参考价值,有需要的朋友可以参考一下,希望对你有所帮助. 环境:python 3.6.8 以某米赛 ...
- python变量需要声明吗_python中可以声明变量类型吗
变量(variable)是Python语言中一个非常重要的概念.变量的主要作用就是为Python程序中的某个值起一个名字.类似于"张三"."李四"." ...
- python元类的使用_python中元类用法实例
本文实例讲述了python中元类用法,分享给大家供大家参考.具体方法分析如下: 1.元类(metaclass)是用来创建类的类 2.type(object):返回一个对象的类型,与object.__c ...
- python可以实现哪些功能_Python中实现机器学习功能的四种方法介绍
本篇文章给大家带来的内容是关于Python中实现机器学习功能的四种方法介绍,有一定的参考价值,有需要的朋友可以参考一下,希望对你有所帮助. 在本文中,我们将介绍从数据集中选择要素的不同方法; 并使用S ...
- python报错怎么看_python中的错误如何查看
python常见的错误有 1.NameError变量名错误 2.IndentationError代码缩进错误 3.AttributeError对象属性错误 4.TypeError类型错误 5.IOEr ...
最新文章
- 在eclipse中搭建maven工程(第二种方法)
- 李彦宏喊你来坐出租车,无人驾驶的那种;百度还要继续搞芯片,联手华为的那种...
- 内核层 inlinehook 隐藏进程
- ROS multi-master——multimaster_fkie配置
- 【LeetCode】剑指 Offer 47. 礼物的最大价值
- STM8单片机复位源判断
- linux C 语言的 system
- 单片机课程设计:基于STM32智能交通灯的设计
- visio2007序列号
- mongo分组查询,统计 Group.grouping(),Group.first(),Accumulator(“$sum“, new Integer(1) 含义解释
- 手机app测试点有哪些?手机app测试用例怎么写?错过绝对后悔!
- 分词并去停用词自定义函数:seg_word(sentence)
- php is numeric用法,PHP使用 is_numeric的实例解析
- SSM+天山产业园访客与疫情防控系统 毕业设计-附源码191123
- 3DMAX快速入门 界面介绍【上】
- 24.多线程(等待唤醒机制,volatile,CAS 算法,线程池,定时器,设计模式)
- 购买iPhone手机时的专业术语名词解析
- 孙猴子装饰七十二变篇--装饰器模式C++实现
- 学英语最好的计算机美剧推荐,推荐四部经典的美剧!适合学习英语口语的美剧有哪些?...
- STM32F767--->时钟