GDAL讀寫向量資料-Python__Python

1、GDAL（Geospatial Data Abstraction Library），源柵格空間資料轉換庫。它利用抽象資料模型來表達所支援的各種檔案格式。它還有一系列命令列工具來進行資料轉換和處理
OGR提供對向量資料格式的讀寫支援，它所支援的檔案格式包括：ESRI Shapefiles， S-57， SDTS， PostGIS，Oracle Spatial， Mapinfo mid/mif ， Mapinfo TAB。
OGR包括如下幾部分：
Geometry：類Geometry (包括OGRGeometry等類)封裝了OpenGIS的向量資料模型，並提供了一些幾何操作，WKB(Well Knows Binary)和WKT(Well Known Text)格式之間的相互轉換，以及空間參考系統(投影)。
Spatial Reference：類OGRSpatialReference封裝了投影和基準面的定義。
Feature：類OGRFeature封裝了一個完整Feature的定義，一個完整的Feature包括一個Geometry和Geometry的一系列屬性。
Feature Definition：類OGRFeatureDefn裡面封裝了feature的屬性，類型、名稱及其預設的空間參考系統等。一個OGRFeatureDefn對象通常與一個層(layer)對應。
Layer：類OGRLayer是一個抽象基類，表示資料來源類OGRDataSource裡面的一層要素(Feature)。
Data Source：類OGRDataSource是一個抽象基類，表示含有OGRLayer對象的一個檔案或一個資料庫。

2、eclipse+python+gdal
(1)、python及PyDev安裝
http://blog.csdn.net/cdl2008sky/article/details/8295176

(2)、搭建GDAL開發環境。 

(2.1)、到Python的GDAL首頁 http://pypi.python.org/pypi/GDAL/1.6.1 下載GDAL的python包——GDAL-1.6.1.win32-py2.6.exe。
  安裝GDAL-1.6.1.win32-py2.6.exe檔案，會自動安裝到python的安裝目錄（$:\Python26\Lib\site-packages）。

(2.2)、到GDAL官網 http://www.gdal.org/的http://download.osgeo.org/gdal/win32/1.6/下載gdalwin32exe160.zip，解壓
縮到任意目錄下，建議採用英文目錄，比如我們解壓到D盤根目錄下，如D:\gdalwin32-1.6。

(2.3)、按照README_EXE.TXT檔案的說明，在“我的電腦”->"進階"->"環境變數"。將“D:\gdalwin32-1.6\bin”添加到"Path"變數
中。並新增一個環境變數GDAL_DATA，變數值設定為 D:\gdalwin32-1.6\data。

(2.4)、 到 http://pypi.python.org/pypi/numpy下載python的numpy模組numpy-1.7.1.win32-py2.6.exe ， 安裝這個包（如果你的python環境未包括這個包，則一定要安裝，因為gdal的一些功能依賴這個包）。

Test:
>>>from osgeo import gdal
>>>from osgeo.gdalconst import *
>>>dataset=gdal.Open(‘d:/test.jpg’,GA_ReadOnly)
>>>dataset.GetDriver().ShortName
>>>'JPEG’

如果沒有錯誤提示，就表明庫已經安裝好了。 

eg:讀寫shap檔案

讀寫步驟
1. import module
2. register driver
3. open datasource
4. get layer
5. get feature
6. ...
7. destroy feature
8. destroy layer
9. destroy datasource

#!/usr/bin/python# -*- coding: gbk -*-'''Created on 2013-8-27@author: chenll'''import os,sysfrom osgeo import gdalfrom osgeo import ogrfrom osgeo import osrimport numpy#讀取shap檔案def readShap():    #為了支援中文路徑，請添加下面這句代碼     gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8","NO")      #為了使屬性工作表欄位支援中文，請添加下面這句    gdal.SetConfigOption("SHAPE_ENCODING","")      #註冊所有的驅動     ogr.RegisterAll()      #資料格式的驅動    driver = ogr.GetDriverByName('ESRI Shapefile')    ds = driver.Open('E:\\arcgis\\data\\Export_Output.shp');    if ds is None:        print 'Could not open 1km地鐵站實驗2.shp'        sys.exit(1)    #擷取第0個圖層    layer0 = ds.GetLayerByIndex(0);    #投影    spatialRef = layer0.GetSpatialRef();    # 輸出圖層中的要素個數      print('要素個數=%d', layer0.GetFeatureCount(0))      print('屬性工作表結構資訊')      defn = layer0.GetLayerDefn()      iFieldCount = defn.GetFieldCount()      for index in range(iFieldCount):          oField =defn.GetFieldDefn(index)          print( '%s: %s(%d.%d)' % (oField.GetNameRef(),oField.GetFieldTypeName(oField.GetType()),oField.GetWidth(),oField.GetPrecision()))         feature = layer0.GetNextFeature()      # 下面開始遍曆圖層中的要素      while feature is not None:          # 擷取要素中的屬性工作表內容          for index in range(iFieldCount):              oField =defn.GetFieldDefn(index)              line =  " %s (%s) = " % (oField.GetNameRef(),oField.GetFieldTypeName(oField.GetType()))              if feature.IsFieldSet( index ):                  line = line+ "%s" % (feature.GetFieldAsString(index))              else:                  line = line+ "(null)"              print(line)          # 擷取要素中的幾何體          geometry =feature.GetGeometryRef()          print geometry        # 為了示範，只輸出一個要素資訊          break      feature.Destroy()    ds.Destroy()    #建立shap檔案def createShap():    #為了支援中文路徑，請添加下面這句代碼     gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8","NO")      #為了使屬性工作表欄位支援中文，請添加下面這句    gdal.SetConfigOption("SHAPE_ENCODING","")      #註冊所有的驅動     ogr.RegisterAll()      #資料格式的驅動    driver = ogr.GetDriverByName('ESRI Shapefile')    ds=driver.CreateDataSource("E:\\arcgis\\point")    shapLayer=ds.CreateLayer("poi",geom_type=ogr.wkbPoint);    #添加欄位    fieldDefn = ogr.FieldDefn('id', ogr.OFTString)     fieldDefn.SetWidth(4)    shapLayer.CreateField(fieldDefn);    #建立feature    defn = shapLayer.GetLayerDefn()    feature = ogr.Feature(defn) ;    #添加屬性    feature.SetField("id","liu")    #添加座標    point = ogr.Geometry(ogr.wkbPoint)    point.AddPoint(float(113.56647912),float(22.16128203))    feature.SetGeometry(point);    shapLayer.CreateFeature(feature)    feature.Destroy()    #指定投影    sr = osr.SpatialReference();    sr.ImportFromEPSG(32612);    prjFile = open("E:\\arcgis\\point\\poi.prj",'w');    sr.MorphToESRI();    prjFile.write(sr.ExportToWkt());    prjFile.close();    ds.Destroy()    def main():    readShap();    createShap();    if __name__ == "__main__":    main();

Python API:http://gdal.org/python/

http://pcjericks.github.io/py-gdalogr-cookbook/