Python+GDAL/OGR矢量数据读写

       常见的矢量数据格式有Shapefile、GeoJSON、CSV，及文件数据库gdb和空间数据库PostGIS，不论是何种格式的数据或如何存储，一旦打开数据源、获取矢量图层后（详情参考OGR操作矢量数据的类结构图），对数据的操作都一样。下面对矢量数据的读写进行详细的介绍。

    一、打开不同的矢量文件

    1、定义打开数据源的函数，并遍历所有的图层，输出他们的名字和图层'''
输出数据源中的图层
参数：fn 数据源的路径
is_write 打开数据源的模式，0 表示只读模式，1 表示读写模式
'''
def print_layers(fn, is_write):
ds = ogr.Open(fn, is_write)
if ds is None:
#raise OSError('Could not open %s', fn)
raise OSError('Could not open {}'.format(fn))
for i in range(ds.GetLayerCount()):
lyr = ds.GetLayer(i)
print('{0}:{1}'.format(i, lyr.GetName()))    2、打开PostGIS数据源

    利用上面打开数据源并输出图层的方法，打开本地的Post GIS数据库 postgis_24_sample，输出图层如下：#打开postgis数据库
print_layers('PG:user=postgres password=postgis dbname=postgis_24_sample', 0)

0:tiger.county
1:tiger.state
2:tiger.place
3:tiger.cousub
4:tiger.edges
5:tiger.addrfeat
6:tiger.faces
7:tiger.zcta5
8:tiger.tract
9:tiger.tabblock
10:tiger.bg

    3、打开文件夹数据源（shapefiles和CSV)

    打开shapefiles文件夹，输出的结果如下所示：#打开shapefile文件夹
shp_fn = r'E:\HZZ_GXDP_Spatial_Data'
print_layers(shp_fn, 0)

0:CommunityRiver
1:CountyRiver
2:StreetRiver
3:Video
4:WRY

    二、读取矢量数据

        读取矢量数据的步骤包括：
        1、打开数据源；
        2、获取图层；
        3、获取图层中的要素；
        4、获取要素的属性和几何信息。
        下面是读取矢量数据的代码片段：import sys
from osgeo import ogr
import ospybook as pb

fn = r'D:\soft\geoprocessing-with-python\china_basic_map'
ds = ogr.Open(fn, 0)
if ds is None:
sys.exit('Could not open {0}.'.format(fn))

'''索引获取数据源中的图层方式'''
lyr = ds.GetLayer(0)
print(lyr.GetName())

i = 0
for fea in lyr:
pt = fea.geometry()
x = pt.GetX()
y = pt.GetY()
'''字段名称 获取属性值'''
code = fea.GetField('GBCODE')
'''对象方式获取属性值'''
length = fea.LENGTH
'''对象方式获取属性值2'''
lpoly_ = fea['LPOLY_']
'''索引方式获取属性值'''
fnode_ = fea.GetField(0)
'''获取特定类型的属性值'''
str_len = fea.GetFieldAsString('LENGTH')
print(code, length, fnode_, lpoly_, str_len, x, y)
i += 1
if i == 20:
break

'''名称获取数据源中的图层'''
lyr2 = ds.GetLayer('国家')
print(lyr2.GetName())    三、获取数据元数据
        矢量数据的元数据，包括数据范围、几何类型、空间参考及图层对象的概要信息等，下面的代码片段展示如何获取这些元数据信息：import sys
from osgeo import ogr

fn = r'D:\soft\geoprocessing-with-python\china_basic_map'
ds = ogr.Open(fn, 0)
if ds is None:
sys.exit('Could not open {0}.'.format(fn))

lyr = ds.GetLayer(0)
'''获取图层范围'''
extent = lyr.GetExtent()
print(extent)

'''(73.44696044921875, 135.08583068847656, 3.408477306365967, 53.557926177978516)'''
'''(min_x, max_x, min_y, max_y)'''

'''获取图层集合类型，返回值为数字'''
'''1:点， 2：线， 3：面'''
geom_type = lyr.GetGeomType()
print(geom_type) #2
print(geom_type == ogr.wkbPoint) #False
print(geom_type == ogr.wkbLineString) #True
print(geom_type == ogr.wkbPolygon) #False

'''通过要素获取几何类型'''
fea = lyr.GetFeature(0)
print(fea.geometry().GetGeometryType()) #2
print(fea.geometry().GetGeometryName()) #LINESTRING

'''获取图层的空间参考'''
print(lyr.GetSpatialRef())
# GEOGCS["GCS_WGS_1984",
# DATUM["WGS_1984",
# SPHEROID["WGS_84",6378137.0,298.257223563]],
# PRIMEM["Greenwich",0.0],
# UNIT["Degree",0.0174532925199433],
# AUTHORITY["EPSG","4326"]]

print(lyr.schema)

'''获取图层属性名称及数据类型'''
for field in lyr.schema:
print(field.name, field.GetTypeName())
# FNODE_ Integer64
# TNODE_ Integer64
# LPOLY_ Integer64
# RPOLY_ Integer64
# LENGTH Real
# BOU2_4M_ Integer64
# BOU2_4M_ID Integer64
# GBCODE Integer

    四、矢量数据写入

        1、矢量数据写入的思路如下：
        （1）以读写模式打开（创建）数据源；
        （2）获取待添加要素的图层（或新创建图层以添加要素）；
        （3）创建空要素，并为几何对象和属性赋值；
        （4）将创建的要素插入图层中。
        下面是利用一个图层中的要素创建一个新图层的详细代码：# -*- coding:utf-8 -*-
import sys
from osgeo import ogr

'''创建一个图层（根据一个图层中的要素创建）'''

fn = r'D:\soft\geoprocessing-with-python\china_basic_map'
'''以读写模式打开数据源'''
ds = ogr.Open(fn, 1)
if ds is None:
sys.exit('Could not open {0}.'.format(fn))

in_lyr = ds.GetLayer('省会城市')

# lyr = ds.GetLayer('国家')
#
# from ospybook.vectorplotter import VectorPlotter
# vp = VectorPlotter(False)
# vp.plot(lyr, fill=False)
# vp.draw()

'''若存在同名图层，则先删除再创键'''
if ds.GetLayer('capital_city'):
ds.DeleteLayer('capital_city')
out_lyr = ds.CreateLayer('capital_city',
in_lyr.GetSpatialRef(),
ogr.wkbPoint)

'''根据元图层字段对象创建新图层的字段'''
out_lyr.CreateFields(in_lyr.schema)

'''根据图层定义创建一个空的要素'''
out_defn = out_lyr.GetLayerDefn()
out_fea = ogr.Feature(out_defn)
for in_fea in in_lyr:
geom = in_fea.geometry()
out_fea.SetGeometry(geom)
for i in range(in_fea.GetFieldCount()):
value = in_fea.GetField(i)
if i != 5:
out_fea.SetField(i, value)
else:
print(value)
out_lyr.CreateFeature(out_fea)        2、创建新的数据源
        创建新数据源的关键在使用正确的驱动程序，每种驱动程序只处理操作一种类型的矢量数据。有两种方式获取驱动程序：（1）从一个已经打开的数据集中获取，这将允许创建一个新的数据源，它和已存在的数据源矢量数据格式一样；（2）使用OGR中的GetDriverByName函数，传递给它驱动程序的简称，驱动程序的名称在OGR网站有介绍。下面是代码片段：''''''''''创建新的数据源'''''''''
import sys
from osgeo import ogr

fn = r'D:\soft\geoprocessing-with-python\china_basic_map'

'''获取数据驱动的第一种方式：根据已有数据源获取'''
ds = ogr.Open(fn, 0)
if ds is None:
sys.exit('Could not open {0}.'.format(fn))

driver = ds.GetDriver()
print(driver.GetName())

'''获取数据源的第二种方式：GetDriverByName'''
json_driver = ogr.GetDriverByName('GeoJSON')
print(json_driver.GetName())

'''创建GeoJson数据源，数据源路径应到具体文件名'''
json_fn = r'D:\soft\geoprocessing-with-python\china_basic_map\json_fn.json'
json_ds = json_driver.CreateDataSource(json_fn)
if json_ds is None:
sys.exit('Could not create {0}'.format(json_fn))
print(json_ds)        3、新建属性字段
        要将一个字段添加到图层中，需要一个包含字段名称、数据类型、字段和精度等重要信息的FieldDefn对象，下面是新建属性字段的代码片段：
'''UseExceptions'''
ogr.UseExceptions()
json_fn = r'D:\soft\geoprocessing-with-python\china_basic_map\json_fn4.json'
json_driver = ogr.GetDriverByName('GeoJSON')
print('start')

try:
json_ds = json_driver.CreateDataSource(json_fn)
'''新建属性数据字段'''
lyr = json_ds.CreateLayer('layer')
coord_fld = ogr.FieldDefn('X', ogr.OFTReal)
coord_fld.SetWidth(8)
coord_fld.SetPrecision(3)
lyr.CreateField(coord_fld)
coord_fld.SetName('
d044
Y')
lyr.CreateField(coord_fld)

#新建要素
fea = ogr.Feature(lyr.GetLayerDefn())
fea.SetField('X', 12.34)
fea.SetFID(1)
lyr.CreateFeature(fea)
#结果同步至硬盘（保存到文件中）
json_ds.SyncToDisk()
except RuntimeError as e:
print(e)

print('end')

    五、更新现有数据

        在进行矢量数据处理时，有时需要更新现有数据，而不是创建一个全新的数据集。是否能够更新以及支持的编辑操作，取决于数据格式。常见的更新数据操作包括 改变图层定义（属性字段）、要素添加、更新和删除等。下面分别进行说明：
        1、改变图层定义

        包括修改属性字段的定义、新增属性字段及删除属性字段，具体代码片段如下：
import sys
from osgeo import ogr
'''''''更新现有数据'''''''
fn = r'D:\soft\geoprocessing-with-python\china_basic_map'
ds = ogr.Open(fn, 1)
if ds is None:
sys.exit('Could not open {0}.'.format(fn))

'''更改图层定义(新增、删除、修改字段)'''
lyr = ds.GetLayer(0)
lyrdefn = lyr.GetLayerDefn()
i = lyrdefn.GetFieldIndex('GBCODE')
'''获取字段的数据类型'''
fld_type = lyrdefn.GetFieldDefn(i).GetType()
print(fld_type)
'''定义新的属性字段'''
fld_defn = ogr.FieldDefn('GBCODECODE', fld_type)
fld_defn2 = ogr.FieldDefn('NEWFIELD4', ogr.OFTInteger)
'''使用ALTER_NAME_FLAG更改属性字段名称，必须保证字段的数据类型一致'''
lyr.AlterFieldDefn(i, fld_defn, ogr.ALTER_NAME_FLAG)
'''新建属性字段'''
lyr.CreateField(fld_defn2)
'''删除字段(参数为字段索引值)'''
lyr.DeleteField(lyr.FindFieldIndex('NEWFIELD', 0))
lyr.DeleteField(lyrdefn.GetFieldIndex('NEWFIELD2'))        2、要素添加、更新和删除'''增加、删除和更新要素'''
fn = r'D:\soft\geoprocessing-with-python\china_basic_map'
ds = ogr.Open(fn, 1)
if ds is None:
sys.exit('Could not open {0}.'.format(fn))

#更新要素（新增属性字段）
lyr = ds.GetLayer(0)
lyr.CreateField(ogr.FieldDefn('NEWFIELD', ogr.OFTInteger))
n = 1
for fea in lyr:
fea.SetField('NEWFIELD', n*n)
lyr.SetFeature(fea)
n += 1

#删除要素
for fea in lyr:
print(fea.GetField('NEWFIELD'))
if fea.GetField('NEWFIELD') == 3179089:
lyr.DeleteFeature(fea.GetFID())

print(lyr.GetFeatureCount())        对数据进行更新操作后，一般需要对数据库进行压缩、重组或重新计算数据的空间范围等，如下：ds.ExecuteSQL('REPACK ' + lyr.GetName()) #压缩数据库
ds.ExecuteSQL('VACUUM') #重组功能

ds.ExecuteSQL('RECOMPUTE EXTENT ON ' + lyr.GetName()) #重新计算图层的空间范围
print(lyr.GetExtent())