2. GDAL python教程（1）——用OGR讀寫矢量數據?

本教程的講義和源碼都是取自Utah State University的openGIS課程

相關資料，包括講義、源碼、數據樣例，請從此處下載

本人只是做點翻譯，寫寫學習體會而已，版權屬于原作者。

歡迎轉載，不過別忘了上面這段話。

2.1. 為什么用open source？?

優點

  1. 免費，適合個人和小公司
  2. 強大的開發工具，找bug更容易
  3. 跨平臺，windows和linux都能用
  4. 拉風！

缺點

  1. 沒有內嵌地理處理器
  2. 用的人少

2.2. Open source RS/GIS模塊?

  1. OGR矢量庫：簡單的矢量數據讀寫，是GDAL的一部分
  2. GDAL地理空間數據抽象庫： a) 讀寫柵格數據 b) ArcGIS也是基于GDAL開發的 c) C++庫，但是可以用python調用

2.3. 相關模塊?

  1. Numeric：高速的數組處理，對柵格數據尤其重要
  2. NumPy：下一代的Numeric
  3. 更強大的gis庫 http://www.gispython.org/

2.4. 導入庫：?

import ogr

或者：

from osgeo import ogr

萬能的方法是：

try:
    from osgeo import ogr
except:
    import ogr

要讀取某種類型的數據，必須要先載入數據驅動，也就是初始化一個對象，讓它“知道”某種數據結構。

import ogr
driver = ogr.GetDriverByName(‘ESRI Shapefile’)

數據驅動driver的open()方法返回一個數據源對象

open(<filename>, <update>)

其中update為0是只讀，為1是可寫

例如：

from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile')
filename = 'C:/Users/gongwei/Documents/My eBooks/python_and_sage/GDAL python/test/ospy_data1/sites.shp'
dataSource = driver.Open(filename,0)
if dataSource is None:
    print 'could not open'
    sys.exit(1)
print 'done!'

注意filename一定要寫絕對路徑！

因為一定要用絕對路徑，為了簡化代碼，經常會使用到os.chdir()

2.5. 讀取數據層?

layer = dataSource.GetLayer(0)

一般ESRI的shapefile都是填0的，如果不填的話默認也是0.

再看看這個數據層里面有幾個點呢？

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

2.5.1. 讀出上下左右邊界?

extent = layer.GetExtent()
print 'extent:', extent
print 'ul:', extent[0], extent[3]
print 'lr:', extent[1], extent[2]

讀取某一要素feature（總算切入正題了），這里讀取的是一個點

feat = layer.GetFeature(41)
fid = feat.GetField('id')
print fid
feat = layer.GetFeature(0)
fid = feat.GetField('id') #should be a different id
print fid

另外還有按順序讀取feature，循環遍歷所有的feature

feat = layer.GetNextFeature()  #讀取下一個
while feat:
    feat = layer.GetNextFeature()
later.ResetReading()  #復位

2.5.2. 提取feature的幾何形狀?

geom = feat.GetGeometryRef()
geom.GetX()
geom.GetY()
print geom.

2.6. 釋放內存?

feature.Destroy()

關閉數據源，相當于文件系統操作中的關閉文件

dataSource.Destroy()

2.7. 讀完了再說怎么寫?

2.7.1. 創建新文件?

driver.CreateDataSource(<filename>)

但是這個文件不能已經存在了，否則會出錯

創建新的layer

dataSource.CreateLayer(<name>,CreateLayer(<name>, geom_type=<OGRwkbGeometryType>, [srs])

舉個例子：

ds2 = driver.CreateDataSource('test.shp')
layer2 = ds2.CreateLayer('test', geom_type=ogr.wkbPoint)

要刪除一個shp文件

driver.DeleteDataSource(‘test.shp’)

要添加一個新字段，只能在layer里面加，而且還不能有數據

添加的字段如果是字符串，還要設定寬度

fieldDefn = ogr.FieldDefn('id', ogr.OFTString)
fieldDefn.SetWidth(4)
layer.CreateField(fieldDefn)

添加一個新的feature，首先得完成上一步，把字段field都添加齊了

然后從layer中讀取相應的feature類型，并創建feature

featureDefn = layer.GetLayerDefn()
feature = ogr.Feature(featureDefn)

設定幾何形狀

feature.SetGeometry(point)

設定某字段的數值

feature.SetField('id', 23)

將feature寫入layer

layer.CreateFeature(feature)