# 2. GDAL python教程（1）——用OGR读写矢量数据¶

## 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’)
```

```open(<filename>, <update>)
```

```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!'
```

## 2.5. 读取数据层¶

```layer = dataSource.GetLayer(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]
```

```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
```

```feat = layer.GetNextFeature()  #读取下一个
while feat:
feat = layer.GetNextFeature()
```

### 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>)
```

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

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

driver.DeleteDataSource(‘test.shp’)

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

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

```feature.SetGeometry(point)
```

```feature.SetField('id', 23)
```

```layer.CreateFeature(feature)
```