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.org1. GDAL庫的安裝
首先要去下面的網站
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. 讀取矢量數據
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。整理一下可見下圖:
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是用來控制是否填充要素,默認是true2.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格式的要求一致。
2.8 遍歷要素
- 使用GetNextFeature方法可以省去使用For循環按ID讀取的低效率;
- 要使用個try except機制,不然再最后一個要素讀完之后,GetNextFeature方法仍然會讀下一個空要素,這時輸出面積會報錯。
- ResetReading函數是用來復位的,不然下次使用GetNextFeature程序接著上次讀的位置繼續讀。
也可以通過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的面積,默認單位為㎡。
2.10 釋放內存
- 要素.Destory是先關閉單個要素,后面的Destory是關閉整個DataSource;
- 關閉數據源,相當于文件系統操作中的關閉文件。
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')總結
以上是生活随笔為你收集整理的python 读取地震道头数据_【Python】OGR库(1):读取矢量数据的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 永恒python配合什么主武器好_学点p
- 下一篇: 共阳数码管段码表_简单共阴极数码管电路图