GDAL读写矢量数据-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/