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

這里推薦兩個比較好的OGR學習資源網(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對應的地址

下載gdal對應python版本的whl文件，這里因為我是python 3.7 64位的，所以下載相應版本就行了：

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

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

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

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

2.1 定義Driver

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

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

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

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

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

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

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

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

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

from osgeo import ogr driver = ogr.GetDriverByName('ESRI Shapefile') #這個函數(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ù)會根據(jù)文件后綴名自動選擇driver進行數(shù)據(jù)讀取。

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

2.4 讀取要素個數(shù)

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

可以通過dir函數(shù)獲取Layer可以使用的方法。這里使用GetLayer方法獲取shp數(shù)據(jù)的圖層，再對其使用GetFeatureCount方法可以獲取shp數(shù)據(jù)中的元素個數(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包。這個包是書《Geoprocessing with python》作者開發(fā)的。可以使用print_attributes來打印出數(shù)據(jù)屬性表，可以直接輸入文件名也可以輸入layer。

不建議使用這個函數(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)用圖層，效果等價

2.5.2 查看圖形

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

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

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

2.6 讀取四至范圍

使用GetExtent方法獲得四至信息，結(jié)果是一個元組，順序為左右下上。若shp數(shù)據(jù)本身含有投影坐標，則輸出的也是投影坐標系的值。

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循環(huán)按ID讀取的低效率；
• 要使用個try except機制，不然再最后一個要素讀完之后，GetNextFeature方法仍然會讀下一個空要素，這時輸出面積會報錯。
• ResetReading函數(shù)是用來復位的，不然下次使用GetNextFeature程序接著上次讀的位置繼續(xù)讀。

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

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

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函數(shù)可以查看geom可以使用的方法；
• GetX和GetY可以直接打印一個個點的坐標；
• 使用geom.Area()可以讀feature的面積，默認單位為㎡。

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

2.10 釋放內(nèi)存

• 要素.Destory是先關(guān)閉單個要素，后面的Destory是關(guān)閉整個DataSource；
• 關(guān)閉數(shù)據(jù)源，相當于文件系統(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)容還不錯，歡迎將生活随笔推薦給好友。