python 读取地震道头数据_【Python】OGR库(1):读取矢量数据

OGR库是一个非常流行的处理地理空间矢量数据的开源库。它可以读取丰富的数据格式,允许用户进行几何处理、属性表操作、数据分析,是个非常强大的开源GIS库。目前OGR已集成在GDAL库中,可以说是GIS的本源之一了,有大量的软件用到了这个库。本篇文章是关于OGR库的一些基础用法汇总,将会持续更新~

这里推荐两个比较好的OGR学习资源网站:

Welcome to the Python GDAL/OGR Cookbook!​pcjericks.github.iohttps://gdal.org/python/​gdal.org

1. GDAL库的安装

首先要去下面的网站

https://www.lfd.uci.edu/~gohlke/pythonlibs/​www.lfd.uci.edu

找到GDAL对应的地址

aa9de1407221070e3b34de178bb3104d.png

下载gdal对应python版本的whl文件,这里因为我是python 3.7 64位的,所以下载相应版本就行了:

0642c2d40d10a6ee7ea3ff4abe022390.png

下载到本地后,对文件存放的文件夹按住shift+鼠标右键,选择“在此处打开Powershell窗口”

ec67119012abbe67a2ff04f5d8b10ac1.png

然后在弹出来的对话框中键入“pip install 文件名”按回车,等待安装即可

afe99fe344626e807304f83cdd3e6aad.png

在安装好之后,会成功安装osgeo这个包。OGR模块是在osgeo包中的。这个包里的所有模块都是以小写字母命名。

2. 读取矢量数据

2.1 定义Driver

在正式读写之前得先理解一下OGR的对象组织形式,弄清楚OGR类之间的关系。

在打开一个矢量数据前,首先需要定义一个driver,用来告诉OGR你想要读取何种格式的数据。每个driver就是对应一种数据格式,比如‘ESRI Shapefile’就是对应的shp格式。OGR可以读取70多个数据格式,基本常用的都包含了。下面是定义driver的一个示例:

from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile') 

2.2 OGR类结构

在讲如何打开数据之前,最好先了解一下OGR的对象组织逻辑。

使用OGR打开一个矢量数据,如shp文件或GeoJSON文件,会产生一个DataSource对象。该对象包含若干个Layer,每个Layer就是一个要素集。值得注意的是,大多数矢量数据格式一般只有一个Layer,少数有多个(如SpatiaLite格式)。既然Layer是个要素集,可想而知它包含的就是一个个的feature了。Feature的概念稍微学过GIS原理的朋友应该都知道,Arcmap里也经常提到,就是几何对象和属性表的集和。所以,feature包含的就是geometry和attribute。整理一下可见下图:

303122d2cb4bbc58f2da5f87e960d0ea.png
OGR类结构

2.3 读取矢量数据

通过对driver对象调用Open方法,可以打开一个文件,打开的结果就是一个Data Souce对象。通过print可以查看到这个信息shp_test是个osgeo.ogr.DataSource对象。

from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile') #这个函数是不区分大小写的filename = r'F:ZhihuDATA'
shp_test = driver.Open(filename,0) #0是只读,1是可写
if shp_test is None:print ('could not open')sys.exit(1)
print(shp_test)

也可以直接使用ogr.open函数打开文件。该函数会根据文件后缀名自动选择driver进行数据读取。

from osgeo import ogrfilename = r'F:ZhihuDATA'
shp_test = ogr.Open(filename)

2.4 读取要素个数

OGR的Layer概念类似于ArcGIS里的FeatureClass,就是多个同类要素(点、线、多边形)的集和。

可以通过dir函数获取Layer可以使用的方法。这里使用GetLayer方法获取shp数据的图层,再对其使用GetFeatureCount方法可以获取shp数据中的元素个数。

layer = shp_test.GetLayer()
dir(layer)
n = layer.GetFeatureCount()
print ('feature count:', n)

2.5 查看数据

2.5.1 查看属性

可在cmd中使用pip install https://github.com/cgarrard/osgeopy-code/raw/master/ospybook-latest.zip安装ospybook包。这个包是书《Geoprocessing with python》作者开发的。可以使用print_attributes来打印出数据属性表,可以直接输入文件名也可以输入layer。

不建议使用这个函数来打印大型数据的全部属性表,可使用数字来控制打印几行并指定打印的字段名。

import ospybook as pbpb.print_attributes(filename,3,['Shape_Area','ID','long','lat']) pb.print_attributes(lyr,3,['Shape_Area','ID','long','lat']) #既可以调用文件名,也可以调用图层,效果等价

2.5.2 查看图形

ospybook 提供了用于绘图的类。它是基于matplotlib的,所以必须安装matplotlib。

交互式图像在创建之后无需使用draw函数调用就可呈现。声明方法是建立一个VectorPlotter类的新实例。

import os 
os.chdir(r'F:ZhihuDATA') #设置工作空间
from ospybook.vectorplotter import VectorPlottervp = VectorPlotter(True) #声明图像为交互式图像
vp.plot(filename,fill=False) #fill是用来控制是否填充要素,默认是true

2.6 读取四至范围

使用GetExtent方法获得四至信息,结果是一个元组,顺序为左右下上。若shp数据本身含有投影坐标,则输出的也是投影坐标系的值。

extent = layer.GetExtent()
print ('extent:', extent)
print ('x range:', extent[0], extent[1])
print ('y range:', extent[2], extent[3])

2.7 读取单个要素

  • 使用GetFeature方法按照FID读取要素,这里读取的第二个要素,即FID=1的那个要素。
  • 通过GetField方法可读取要素指定列信息,值得注意的是这里需要输入的列名不分大小写,同shp格式的要求一致。
feat = layer.GetFeature(1)
fid = feat.GetField('id') 
area = feat.GetField('shape_area') print (fid)
print (area)
dir(feat)

2.8 遍历要素

  • 使用GetNextFeature方法可以省去使用For循环按ID读取的低效率;
  • 要使用个try except机制,不然再最后一个要素读完之后,GetNextFeature方法仍然会读下一个空要素,这时输出面积会报错。
  • ResetReading函数是用来复位的,不然下次使用GetNextFeature程序接着上次读的位置继续读。
feat = layer.GetNextFeature()  #读取下一个
while feat:feat = layer.GetNextFeature()try:area = feat.GetField('shape_area') print(area)except:print('Done!')
layer.ResetReading()  #复位

也可以通过for循环直接遍历每个要素。但值得注意的一点是当前要素问题。比如我们用For循环遍历了一遍Layer中的所有feature并输出了它们一个字段的值。若想再通过for循环遍历一遍输出另一个字段的值,你会发现得不到任何结果。这是因为在第一次for循环后,指针停在了最后一个feature上,必须使用Layer.ResetReading()函数来进行重置。

for i,feat in enumerate(layer):if i >= 5:breakx=feat.geometry().GetX()y=feat.geometry().GetY()fid = feat.IDarea = feat.GetField('shape_area') print (fid,area,x,y)
layer.ResetReading()  #复位

2.9 提取要素几何信息

  • 使用GetGeometryRef方法读取要素几何信息,通过dir函数可以查看geom可以使用的方法;
  • GetX和GetY可以直接打印一个个点的坐标;
  • 使用geom.Area()可以读feature的面积,默认单位为㎡。
feat = layer.GetFeature(1)
geom = feat.GetGeometryRef()print(geom)
print(geom.Area())

2.10 释放内存

  • 要素.Destory是先关闭单个要素,后面的Destory是关闭整个DataSource;
  • 关闭数据源,相当于文件系统操作中的关闭文件。
feat.Destroy()
shp_test.Destroy()

2.11 删除文件

使用DeleteDataSource可以删除shp文件及其附属文件(如dbf,poj等文件)。

import os
filename = 'F:/Zhihu/DATA/testCopy.shp'
if os.path.exists(filename):driver.DeleteDataSource(filename)print('File was deleted!')
else:print('File not exist')

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/567095.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

进入其他网络共享计算机,局域网内如何进入其他电脑,两个电脑利用无线建立局域网-...

虽然大家平时都在使用电脑,但是大家中的相当一部分朋友们从来都没想过应该如何通过局域网进入到别人的电脑这个问题。怎么样?听起来是不是很神奇呢?其实这种方法从电脑能狗互联的时候就已经诞生了,只是大家一般不经常使用这种工具,所以就不…

appcrash事件怎么解决_解决问题的最佳办法,是让问题不再是问题

我们常常陷入迷局,绕来绕去却怎么也找不到出路。因为,当事者迷,旁观者清。身在局中,我们的思维很容易就被固定在既定的逻辑里。有的是以往的经验总结,有些是从他人处习得。对于我们已经习得的东西,在遇到事…

台式计算机的硬件组成部分,台式电脑主机的硬件组成部分简介

计算机硬件系统中用于放置主板及其他主要部件的容器(Mainframe)。通常包括 CPU、内存、硬盘、光驱、电源、以及其他输入输出控制器和接口,如 USB 控制器、显卡、网卡、声卡等等。位于主机箱内的通常称为内设,而位于主机箱之外的通常称为外设(如显示器、键…

python 伪造源ip_Swaks伪造邮件

0x00 swaksswaks - Swiss Army Knife SMTP, the all-purpose smtp transaction tester.swaks堪称SMTP协议的瑞士军刀,使用它我们可以灵活的操作SMTP协议报文,这篇文章主要是记录一下我是如何伪造一封邮件绕过gmail的检测。通常最简单的发送命令&#xff…

clientmacaddr进不去系统win10_教你一分钟搞定戴尔电脑WIN10改WIN7

最近有很多人问小编,戴尔的新款电脑WIN10(win8)改WIN7电脑不认U盘,不知道怎么设置,今天小编就给大家分享一个快速进入的方法。首先把装有系统的U盘插入电脑,开机一直按F12进入Bios菜单第二步:进入菜单界面后&#xff0…

win7笔记本外接显示器html,window7笔记本外接显示器只显示一个屏幕怎么设置

许多用户都会偏向于入手win7笔记本电脑,这样电脑携带起来也是非常方便,不过由于屏幕较小的缘故,有用户就会选择外接一台显示器,不过在给win7笔记本外接显示器之后就需要对于其进行设置只显示一个屏幕,接下来小编就来教…

苹果手机怎么查看足迹_用了5年苹果手机!才知道查看一个字母就能辨别手机真假...

苹果手机的价格一般都比较贵,所以大家都怕买到假货!今天笔者就教大家如何快速分辨苹果手机真假,只需查看一个字母就能知道手中的苹果手机是什么型号。方法一:桌面图标iPhone手机桌面上的时钟图标比较特别,它的时针会随…

java实现矩阵谱峰搜索算法

矩阵谱峰搜索算法,也称为矩阵谱峰查找算法,是一种用于搜索二维矩阵中谱峰的方法。谱峰是指在矩阵中的一个元素,它比其上下左右四个相邻元素都大或相等。 该算法的基本思想是从矩阵的中间列开始,找到该列中的最大元素,…

电脑中计算机右键管理无法打开,win8系统计算机右键菜单中的管理打不开怎么办...

‍‍计算机管理一组Windows管理工具,这些工具被组合到一个控制台中,方便我们操作。最近有些雨林木风win8旗舰版用户遇到了计算机管理打不开的情况,在右键点击计算机打开菜单后,点击管理打不开,遇到这种问题该怎么办呢&…

中海达gps软件wince_应用|无人机航测15分钟能做啥?中海达PPK告诉你答案

标星置顶,一秒找到中海达讯点击上方“中海达讯”→点击右上角“…”→点选“设为星标 ★”在航测作业中快速现场成图生成快拼成果报告快速通过内方位元素精度评估完成以上步骤你最快要多长时间?15分钟这是中海达PPK套装给出的答案点击视频查看中海达PPK套…

机械制造工艺基础_机械制造工艺基础知识,錾削与锯削加工工艺

一、錾削用锤子打击錾子对金属工件进行切削加工1.錾削工具(1)錾子錾子的种类及用途(2)锤子2.錾削时的几何角度(1)楔角(βo)錾削硬度较高材料: βo60~70錾削软材料:βo30~50錾削中等硬度材料:βo…

excel怎么更改坐标轴刻度_如何用excel制作帕累托图

帕累托图,也叫排列图/帕拉图/主次图,是一种将出现的质量问题和质量改进项目按照重要程度依次排列而采用的图表。当我们的帕累托图完成时,便可辅助我们直观的找到造成问题的主要原因,进而针对问题实施对策,最终达到改善…

c++矩阵连乘的动态规划算法并输出_你在Java中用过动态规划吗?

1. 介绍动态规划典型的被用于优化递归算法,因为它们倾向于以指数的方式进行扩展。动态规划主要思想是将复杂问题(带有许多递归调用)分解为更小的子问题,然后将它们保存到内存中,这样我们就不必在每次使用它们时重新计算它们。要理解动态规划的…

spark中dataframe解析_SparkSql 中 JOIN的实现

Join作为SQL中一个重要语法特性,几乎所有稍微复杂一点的数据分析场景都离不开Join,如今Spark SQL(Dataset/DataFrame)已经成为Spark应用程序开发的主流,作为开发者,我们有必要了解Join在Spark中是如何组织运行的。SparkSQL总体流程…

含枚举类型的函数声明_02Golang基础类型

基础类型命名Go语言中的函数名、变量名、常量名、类型名、语句标号和包名等所有的命名,都遵循一个简单的命名规则:一个名字必须以一个字母(Unicode字母)或下划线开头,后面可以跟任意数量的字母、数字或下划线。大写字母…

将xscj指定为当前数据库_通过网络连接数据库模式Hive的搭建过程详解

最近在搭建通过网络直接连接数据库模式的Hive时总是在启动的时候报各种错误,所以今天,我们来总结一下这种模式的Hive的搭建过程。【数据库安装】安装mysqlyum install mysql-server -y配置:启动mysql服务:service mysqld start设置…

游戏脚本代码大全_按键精灵】一个很好学的脚本

这【按键精灵】一个很好学的脚本命令名称:GetPixelColor 得到指定点颜色命令功能:得到指定位置的点的颜色命令参数:参数1 整数型,屏幕X坐标参数2 整数型,屏幕Y坐标返 回 值:【按键精灵】一个很好学的脚本字…

标记三维点_便携式3D扫描仪全自动三坐标测量机三维扫描设计扫描测量摄影

项目简介客户产品该客户的产品是铝铸件,铸件的很多位置没有太高的精度要求,但是铸件加工出来的孔位需要严格对上装配的位置,精加工的面要求却又特别高,并且孔位之间相距较远。客户的困难 由于产品要求较高,一般卡尺量具…

python画图显示不了中文_Python使用matplotlib绘图无法显示中文问题的解决方法

本文实例讲述了Python使用matplotlib绘图无法显示中文问题的解决方法。分享给大家供大家参考,具体如下: 在python中,默认情况下是无法显示中文的,如下代码: import matplotlib.pyplot as plt # 定义文本框和箭头格式 d…

软件测试语句覆盖,软件测试中的语句覆盖,分支覆盖,条件覆盖以及路径覆盖...

我举一个简单的例子来解释一下语句覆盖,分支覆盖,条件覆盖以及路径覆盖的相关知识,如果有不对的地方,恳请各位同行指正:举例说明:if Atrue and Btrue then Action1if Ctrue or Dtrue then Action2这是一个很…