GDAL-矢量数据读取与写入

1.数据读取

fn = r'E:\gdal\osgeopy-data\osgeopy-data-global\osgeopy-data\global\ne_50m_populated_places.shp'
ds = ogr.Open(fn, 0)
if ds is None:
    sys.exit('Could not open {0}.'.format(fn))
lyr = ds.GetLayer(0)
lyr.ResetReading() 

获取图层中feature的数量
例如图层中有100个点.num_features=100

num_features = lyr.GetFeatureCount()
print('图层中特征的数量{}'.format(num_features))

根据feature索引提取最后一个feature

last_feature = lyr.GetFeature(num_features - 1)
print(last_feature.NAME)

2.获取图层属性与几何数据

遍历feature,并访问feature的属性与几何
获取图层字段名称与类型

for field in lyr.schema:
    print(field.name, field.GetTypeName())
lyr.ResetReading() 

遍历图层中每个feature的字段属性

for feat in lyr:
    print(feat.GetField('NAME'), feat.GetField('POP2050'))

遍历图层中每个feature的几何属性

for feat in lyr:
    pt = feat.geometry()
    print(feat.GetField('NAME'), pt.GetX(), pt.GetY())

3.访问图层的元数据

获取空间范围.

lyr.ResetReading()
extent = lyr.GetExtent()
print(extent)
print('Upper left corner: {}, {}'.format(extent[0], extent[3]))
print('Lower right corner: {}, {}'.format(extent[1], extent[2]))

图层属性：获取几何类型

print(lyr.GetGeomType())
print(lyr.GetGeomType() == ogr.wkbPoint)
print(lyr.GetGeomType() == ogr.wkbPolygon)

获取feature的几何类型名称

feat = lyr.GetFeature(0)
print(feat.geometry().GetGeometryName())

获取图层的空间参考

print(lyr.GetSpatialRef())

获取图层中的字段与类型名称

for field in lyr.schema:
    print(field.name, field.GetTypeName())

4. 矢量文件创建

4.1根据数据驱动类型，创建数据驱动文件

json_driver = ogr.GetDriverByName('GeoJSON')
print(json_driver.name)
esri_driver = ogr.GetDriverByName('ESRI Shapefile')
print(esri_driver.name)

4.2根据数据驱动文件创建对用数据类型的地理数据文件

json_fn = os.path.join(data_dir, 'example.geojson')
json_ds = json_driver.CreateDataSource(json_fn)
if json_ds is None:# 无法创建
    #  数据创建的断言，判定错误类型
    ogr.UseExceptions()
    fn = os.path.join(data_dir,'africa.geojson')
    driver = ogr.GetDriverByName('GeoJSON')
    print('Doing some preliminary analysis...')

    try:
        # This will fail if the file already exists
        ds = driver.CreateDataSource(fn)
        lyr = ds.CreateLayer('layer')
        # Do more stuff, like fields and save data
    except RuntimeError as e:
        # This runs if the data source already exists and an error was raised
        print(e)

    # print('Doing some more analysis...')
    # 处理代码，若文件已存在
    if os.path.exists(json_fn):
        json_driver.DeleteDataSource(json_fn)
    else:
    # 若文件不存在
        print('文件路径是否可以创建')
        sys.exit('Could not create {0}.'.format(json_fn))

4.3将已有矢量文件另存为新的矢量文件

import os
os.environ['PROJ_LIB'] = r'D:\Programfile\conda3\envs\RS\lib\site-packages\osgeo\data\proj'

上述两行代码用于处理bug;

RuntimeError: PROJ: proj_identify: D:\Programfile\conda3\envs\RS\Library\share\proj\proj.db lacks DATABASE.LAYOUT.VERSION.MAJOR / DATABASE.LAYOUT.VERSION.MINOR metadata.
It comes from another PROJ installation.
读取文件

in_fn = r'E:\gdal\osgeopy-data\osgeopy-data-washington\osgeopy-data\Washington\large_cities.shp'
in_ds = ogr.Open(in_fn, 0)
if in_ds is None:
    sys.exit('Could not open {0}.'.format(in_fn))
in_lyr = in_ds.GetLayer(0)

创建输出文件

# 获取原始文件的数据驱动类型
~~~ python
driver = in_ds.GetDriver()
out_fn = os.path.join(data_dir, 'precision_test.shp')
if os.path.exists(out_fn):
    driver.DeleteDataSource(out_fn)
out_ds = driver.CreateDataSource(out_fn)
if out_ds is None:
    sys.exit('Could not create {0}.'.format(out_fn))

在数据源中创建先创建图层,并设置图层属性

out_lyr = out_ds.CreateLayer('precision_test',
                             in_lyr.GetSpatialRef(),
                             ogr.wkbPoint)

然后为图层创建属性字段
先定义字段类型

name_fld = ogr.FieldDefn('Name', ogr.OFTString)
name_fld.SetWidth(6)

将字段添加到图层中

out_lyr.CreateField(name_fld)
#创建带有默认精度的属性字段.
coord_fld = ogr.FieldDefn('X_default', ogr.OFTReal)
out_lyr.CreateField(coord_fld)
coord_fld.SetName('Y_default')
out_lyr.CreateField(coord_fld)

创建feature，并将图层元数据赋予feature

out_feat = ogr.Feature(out_lyr.GetLayerDefn())

遍历输入特征

for in_feat in in_lyr:
    # 获取输入特征的几何与属性
    pt = in_feat.geometry()
    name = in_feat.GetField('NAME')
    # 为输出特征添加属性
    out_feat.SetGeometry(in_feat.geometry())
    out_feat.SetField('Name', name)
    out_feat.SetField('X_default', pt.GetX())
    out_feat.SetField('Y_default', pt.GetY())
    # out_feat.SetField('X_short', pt.GetX())
    # out_feat.SetField('Y_short', pt.GetY())
    # 为图层添加特征（图层的属性字段会与feature的属性字段相匹配）
    out_lyr.CreateFeature(out_feat)
del out_ds

总结

矢量数据写入的操作相对复杂，在已有参考矢量元数据与属性数据已知的情况下，需要对新建是来那个输出数据实现包括是数据驱动、数据输出、创建图层（投影、几何类型）、未涂层添加字段属性、创建feature，将feature添加到图层中。
需要明确source、layer、feature、字段与几何之间的关系；其中feature的字段名与类型要与layer中的字段名与类型相匹配，layer中的字段类型与feature中的字段名类型是公用集与字子集之间的关系。
注：其中feature设独立的变量（arcgis中的属性表），包含几何与字段属性，创建feature后需要将其添加到图层中，并删除数据源，在完成矢量数据的写入操作。

reference
[python地理数据处理]