OGR庫是一個(gè)非常流行的處理地理空間矢量數(shù)據(jù)的開源庫。它可以讀取豐富的數(shù)據(jù)格式，允許用戶進(jìn)行幾何處理、屬性表操作、數(shù)據(jù)分析，是個(gè)非常強(qiáng)大的開源GIS庫。目前OGR已集成在GDAL庫中，可以說是GIS的本源之一了，有大量的軟件用到了這個(gè)庫。本篇文章是關(guān)于OGR庫的一些基礎(chǔ)用法匯總，將會(huì)持續(xù)更新~

這里推薦兩個(gè)比較好的OGR學(xué)習(xí)資源網(wǎng)站：

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

1. GDAL庫的安裝

首先要去下面的網(wǎng)站

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

找到GDAL對(duì)應(yīng)的地址

下載gdal對(duì)應(yīng)python版本的whl文件，這里因?yàn)槲沂莗ython 3.7 64位的，所以下載相應(yīng)版本就行了：

下載到本地后，對(duì)文件存放的文件夾按住shift+鼠標(biāo)右鍵，選擇“在此處打開Powershell窗口”

然后在彈出來的對(duì)話框中鍵入“pip install 文件名”按回車，等待安裝即可

在安裝好之后，會(huì)成功安裝osgeo這個(gè)包。OGR模塊是在osgeo包中的。這個(gè)包里的所有模塊都是以小寫字母命名。

2. 讀取矢量數(shù)據(jù)

2.1 定義Driver

在正式讀寫之前得先理解一下OGR的對(duì)象組織形式，弄清楚OGR類之間的關(guān)系。

在打開一個(gè)矢量數(shù)據(jù)前，首先需要定義一個(gè)driver，用來告訴OGR你想要讀取何種格式的數(shù)據(jù)。每個(gè)driver就是對(duì)應(yīng)一種數(shù)據(jù)格式，比如‘ESRI Shapefile’就是對(duì)應(yīng)的shp格式。OGR可以讀取70多個(gè)數(shù)據(jù)格式，基本常用的都包含了。下面是定義driver的一個(gè)示例：

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

2.2 OGR類結(jié)構(gòu)

在講如何打開數(shù)據(jù)之前，最好先了解一下OGR的對(duì)象組織邏輯。

使用OGR打開一個(gè)矢量數(shù)據(jù)，如shp文件或GeoJSON文件，會(huì)產(chǎn)生一個(gè)DataSource對(duì)象。該對(duì)象包含若干個(gè)Layer，每個(gè)Layer就是一個(gè)要素集。值得注意的是，大多數(shù)矢量數(shù)據(jù)格式一般只有一個(gè)Layer，少數(shù)有多個(gè)（如SpatiaLite格式）。既然Layer是個(gè)要素集，可想而知它包含的就是一個(gè)個(gè)的feature了。Feature的概念稍微學(xué)過GIS原理的朋友應(yīng)該都知道，Arcmap里也經(jīng)常提到，就是幾何對(duì)象和屬性表的集和。所以，feature包含的就是geometry和attribute。整理一下可見下圖：

OGR類結(jié)構(gòu)

2.3 讀取矢量數(shù)據(jù)

通過對(duì)driver對(duì)象調(diào)用Open方法，可以打開一個(gè)文件，打開的結(jié)果就是一個(gè)Data Souce對(duì)象。通過print可以查看到這個(gè)信息shp_test是個(gè)osgeo.ogr.DataSource對(duì)象。

from osgeo import ogr driver = ogr.GetDriverByName('ESRI Shapefile') #這個(gè)函數(shù)是不區(qū)分大小寫的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函數(shù)打開文件。該函數(shù)會(huì)根據(jù)文件后綴名自動(dòng)選擇driver進(jìn)行數(shù)據(jù)讀取。

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

2.4 讀取要素個(gè)數(shù)

OGR的Layer概念類似于ArcGIS里的FeatureClass，就是多個(gè)同類要素(點(diǎn)、線、多邊形)的集和。

可以通過dir函數(shù)獲取Layer可以使用的方法。這里使用GetLayer方法獲取shp數(shù)據(jù)的圖層，再對(duì)其使用GetFeatureCount方法可以獲取shp數(shù)據(jù)中的元素個(gè)數(shù)。

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

2.5 查看數(shù)據(jù)

2.5.1 查看屬性

可在cmd中使用pip install https://github.com/cgarrard/osgeopy-code/raw/master/ospybook-latest.zip安裝ospybook包。這個(gè)包是書《Geoprocessing with python》作者開發(fā)的。可以使用print_attributes來打印出數(shù)據(jù)屬性表，可以直接輸入文件名也可以輸入layer。

不建議使用這個(gè)函數(shù)來打印大型數(shù)據(jù)的全部屬性表，可使用數(shù)字來控制打印幾行并指定打印的字段名。

import ospybook as pbpb.print_attributes(filename,3,['Shape_Area','ID','long','lat']) pb.print_attributes(lyr,3,['Shape_Area','ID','long','lat']) #既可以調(diào)用文件名，也可以調(diào)用圖層，效果等價(jià)

2.5.2 查看圖形

ospybook 提供了用于繪圖的類。它是基于matplotlib的，所以必須安裝matplotlib。

交互式圖像在創(chuàng)建之后無需使用draw函數(shù)調(diào)用就可呈現(xiàn)。聲明方法是建立一個(gè)VectorPlotter類的新實(shí)例。

import os os.chdir(r'F:ZhihuDATA') #設(shè)置工作空間 from ospybook.vectorplotter import VectorPlottervp = VectorPlotter(True) #聲明圖像為交互式圖像 vp.plot(filename,fill=False) #fill是用來控制是否填充要素，默認(rèn)是true

2.6 讀取四至范圍

使用GetExtent方法獲得四至信息，結(jié)果是一個(gè)元組，順序?yàn)樽笥蚁律?。若shp數(shù)據(jù)本身含有投影坐標(biāo)，則輸出的也是投影坐標(biāo)系的值。

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

2.7 讀取單個(gè)要素

• 使用GetFeature方法按照FID讀取要素，這里讀取的第二個(gè)要素，即FID=1的那個(gè)要素。
• 通過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循環(huán)按ID讀取的低效率；
• 要使用個(gè)try except機(jī)制，不然再最后一個(gè)要素讀完之后，GetNextFeature方法仍然會(huì)讀下一個(gè)空要素，這時(shí)輸出面積會(huì)報(bào)錯(cuò)。
• ResetReading函數(shù)是用來復(fù)位的，不然下次使用GetNextFeature程序接著上次讀的位置繼續(xù)讀。

feat = layer.GetNextFeature() #讀取下一個(gè) while feat:feat = layer.GetNextFeature()try:area = feat.GetField('shape_area') print(area)except:print('Done!') layer.ResetReading() #復(fù)位

也可以通過for循環(huán)直接遍歷每個(gè)要素。但值得注意的一點(diǎn)是當(dāng)前要素問題。比如我們用For循環(huán)遍歷了一遍Layer中的所有feature并輸出了它們一個(gè)字段的值。若想再通過for循環(huán)遍歷一遍輸出另一個(gè)字段的值，你會(huì)發(fā)現(xiàn)得不到任何結(jié)果。這是因?yàn)樵诘谝淮蝔or循環(huán)后，指針停在了最后一個(gè)feature上，必須使用Layer.ResetReading()函數(shù)來進(jìn)行重置。

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() #復(fù)位

2.9 提取要素幾何信息

• 使用GetGeometryRef方法讀取要素幾何信息，通過dir函數(shù)可以查看geom可以使用的方法；
• GetX和GetY可以直接打印一個(gè)個(gè)點(diǎn)的坐標(biāo)；
• 使用geom.Area()可以讀feature的面積，默認(rèn)單位為㎡。

feat = layer.GetFeature(1) geom = feat.GetGeometryRef()print(geom) print(geom.Area())

2.10 釋放內(nèi)存

• 要素.Destory是先關(guān)閉單個(gè)要素，后面的Destory是關(guān)閉整個(gè)DataSource；
• 關(guān)閉數(shù)據(jù)源，相當(dāng)于文件系統(tǒng)操作中的關(guān)閉文件。

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')

總結(jié)

以上是生活随笔為你收集整理的python 读取地震道头数据_【Python】OGR库（1）：读取矢量数据的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。