Python地理数据处理 四：矢量数据读写（二）

1. 获取数据的元数据

  地理空间元数据，即关于数据的数据，它在地理信息中用于描述地理数据集的内容、质量、表示方式、空间参照系、管理方式以及数据集的其他特征等，帮助和促进人们有效地定位、评价、比较、获取和使用地理相关数据，它是实现地理空间数据集共享的核心内容之一。
  可以使用 GetFeatureCount() 函数获取一个图层的要素数量。但只适用于图层，不适用于数据源，因为数据源中的每个图层可以有不同数量的要素、集合类型、空间范围及属性等。
图层的空间范围：上、下、左、右4个方向的最大和最小边界坐标构成的矩形，可以使用GetExtent（）函数来获取图层的边界坐标，返回 （min_x, max_x, min_y, max_y） 的数字元组。
  例如：

>>> import os
>>> from osgeo import ogr
>>> ds = ogr.Open(r'E:\Google chrome\Download\GIS with python\osgeopy data\osgeopy-data-washington\large_cities.geojson')
>>> lyr = ds.GetLayer(0)
>>> extent = lyr.GetExtent()
>>> print(extent)
(-122.66148376464844, -117.4260482788086, 45.638729095458984, 48.759552001953125)

>>> lyr = ds.GetLayer(0)
>>> extent = lyr.GetExtent()
>>> print(extent)
(-175.22056447761656, 179.21664709402887, -89.99999981438727, 78.21668438639699)
>>> print('Upper left corner: {}, {}'.format(extent[0], extent[3]))
Upper left corner: -175.22056447761656, 78.21668438639699
>>> print('Lower right corner: {}, {}'.format(extent[1], extent[2]))
Lower right corner: 179.21664709402887, -89.99999981438727

  获取图层的几何类型，用 GetGeoType（） 函数返回一个整数值，不是可理解的字符串。
  OGR模块具有多个常量， 可以通过 GetGeoType（） 函数将获取的数值与表中的数值进行对比，来检查它属于哪种数据类型。

几何要素类型	OGR常量
Point	wkbPoint
Mulitpoint	wkbMultiPoint
Line	wkbLineString
Multiline	wkbMultiLineString
Polygon	wkbPolygon
Multipolygon	wkbMultiPolygon
Unknown geometry type	wkbUnknown
No geometry	wkbNone

  检查文件类型：

# 交互式
>>> lyr = ds.GetLayer(0)
>>> print(lyr.GetGeomType())
1
>>> print(lyr.GetGeomType() == ogr.wkbPoint)   # 与常量进行比较，点数据
True
>>> print(lyr.GetGeomType() == ogr.wkbPolygon)  # 与常量进行比较，面数据
False

  若图层有多种几何类型，则函数返回 wkbUnknown。

import sys
from osgeo import ogr

fn = r'E:\Google chrome\Download\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)

feat = lyr.GetFeature(0)  # 获取图层第一个要素
print(feat.geometry().GetGeometryName())  # 获取几何类型，并返回名字

POINT

1.1 获取空间参考系统

  lyr.GetSpatialRef()

>>> print(lyr.GetSpatialRef())
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST],
    AUTHORITY["EPSG","4326"]]

1.2 获取图层自身属性字段信息

  lyr.schema（）

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

SCALERANK Integer
NATSCALE Integer
LABELRANK Integer
FEATURECLA String
NAME String
NAMEPAR String
<snip>

1.3 OGR库（补充）

2. 矢量数据写入

import sys
from osgeo import ogr

# 打开写入的数据源
ds = ogr.Open(r'E:\Google chrome\Download\global', 1) # 传入参数1，创建新的文件（.shp文件）
if ds is None:
    sys.exit('Could not open folder.')

# 获取shaefile文件
in_lyr = ds.GetLayer('ne_50m_populated_places') # 将.shp命名，可以不带后缀

# 若图层存在，则删除该图层
if ds.GetLayer('capital_cities'):
    ds.DeleteLayer('capital_cities')

# 创建一个点图层
out_lyr = ds.CreateLayer('capital_cities',
                         in_lyr.GetSpatialRef(),
                         ogr.wkbPoint) # 图层名在文件夹中是唯一的
out_lyr.CreateFields(in_lyr.schema)

# 创建一个空要素
out_defn = out_lyr.GetLayerDefn()
out_feat = ogr.Feature(out_defn)
for in_feat in in_lyr:
    if in_feat.GetField('FEATURECLA') == 'Admin-0 capital':

        # 复制几何要素和属性
        geom = in_feat.geometry()
        out_feat.SetGeometry(geom)
        for i in range(in_feat.GetFieldCount()):
            value = in_feat.GetField(i)
            out_feat.SetField(i, value)

        # 插入该要素
        out_lyr.CreateFeature(out_feat)

# 关闭文件
del ds

  CreateLayer（name，[srs], [geom_type], [option]）:

1. name：新建图层名。
2. srs：新建图层空间参考系统，默认为空，表示未设置任何空间参考系统。数据如果没有空间参考信息，将难以确定要素在地球上的位置。 如KML只支持WGS84基准面。
3. geom_type：几何常量类型，表示图层中要素的几何类型，默认为wkbUnknown。
4. options：只适用于特定的矢量数据类型。

  基于已有.shp文件，创建一个新的.shp文件，并使用相同的空间参考系统。

import sys
from osgeo import ogr

ds = ogr.Open(r'E:\Google chrome\Download\global',1)
if ds is None:
    sys.exit('Could not open folder.')
    
in_lyr = ds.GetLayer('ne_50m_populated_places')   #  获取shapefile文件

if ds.GetLayer('captial_cities'):
    ds.DeleteLayer('captial_cities')

# 创建一个点图层
out_lyr = ds.CreateLayer('captial_cities',in_lyr.GetSpatialRef(),ogr.wkbPoint)
out_lyr.CreateFields(in_lyr.schema)  #  返回属性字段定义列表

# 创建一个空要素
out_defn = out_lyr.GetLayerDefn()
out_feat = ogr.Feature(out_defn)

# 复制几何要素和属性
for in_feat in in_lyr:
    if in_feat.GetField('FEATURECLA') == 'Admin-0 capital':  # 判断是否为首都城市
        geom = in_feat.geometry()
        out_feat.SetGeometry(geom)
        for i in range(in_feat.GetFieldCount()):
            value = in_feat.GetField(i)
            out_feat.SetField(i,value)
        out_lyr.CreateFeature(out_feat)
        
del ds
# ds.SyncToDisk

2.1 创建新数据源

  创建新数据源的关键是正确使用驱动程序，且一种驱动程序只能处理一种类型的矢量数据。如使用GeoJSON驱动程序创建一个shp文件任然是创建一个GeoJSON文件。
  获取驱动程序
  方法一： 从已打开的数据集中获取

from osgeo import ogr
ds = ogr.Open(r'E:\Google chrome\Download\california\california_50m_wind_albers.shp')
driver = ds.GetDriver()

print(driver.name)

ESRI Shapefile

  方法二：使用 OGR中GetDriverByName()函数，传递给它驱动程序的简称。可以使用：ogrinfo --formats 获得驱动程序。

import sys
from osgeo import ogr

esri_fn = r'E:\Google chrome\Download\california'

esri_driver = ogr.GetDriverByName('ESRI Shapefile')
esri_ds = esri_driver.CreateDataSource(esri_fn)
if esri_ds is None:
    sys.exit('Could not create {0}.'.format(esri_fn))

ds = ogr.Open(r'E:\Google chrome\Download\global',1)

if ds is None:
    sys.exit('Could not open folder.')
    
in_lyr = ds.GetLayer('ne_50m_populated_places')

if esri_ds.GetLayer('captial_cities'):
    esri_ds.DeleteLayer('captial_cities')

out_lyr = esri_ds.CreateLayer('captial_cities',in_lyr.GetSpatialRef(),ogr.wkbPoint)
out_lyr.CreateFields(in_lyr.schema)

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

for in_feat in in_lyr:
    if in_feat.GetField('FEATURECLA') == 'Admin-0 capital':
        geom = in_feat.geometry()
        out_feat.SetGeometry(geom)
        for i in range(in_feat.GetFieldCount()):
            value = in_feat.GetField(i)
            out_feat.SetField(i,value)
        out_lyr.CreateFeature(out_feat)

del esri_ds

  数据源和图层创建选项是完全不同的，使用数据源选项创建一个SpatialLite数据源，而不使用SQLite：

import sys
from osgeo import ogr

json_fn = r'E:\Google chrome\Download\global'

json_driver = ogr.GetDriverByName('GeoJSON')
json_ds = json_driver.CreateDataSource(json_fn)
if json_ds is None:
    sys.exit('Could not create {0}.'.format(json_fn))

>>> from osgeo import ogr
>>> driver = ogr.GetDriverByName('SQLite')
>>> ds = driver.CreateDataSource(r'E:\Google chrome\earth.sqlite',['SPATIALITE=yes'])

  注意：创建鑫都数据源时，不能覆盖现有的数据源。如果可能会覆盖，需要在创建之前删除就数据。所以要使用驱动程序删除现有数据源，而不能使用python内置函数。因为驱动程序可以确保所有必要的文件都被删除掉，如shp、dbf、shx等。
  方法：

if os.path.exists(json_fn):
    json_driver.DeleteDataSource(json_fn)
json_ds = json_driver.CreateDataSource(json_fn)
if json_ds is None:
    sys.exit('Could not create {0}.'.format(json_fn))

  OGR例外，ogr.UseExceptions() 将在程序错误的地方自动报错。

import sys
import os
from osgeo import ogr
data_dir = r'E:\Google chrome\Download\california'

# 开启例外
ogr.UseExceptions()

fn = os.path.join(data_dir, 'output', 'natural_earth_50m.sqlite')
driver = ogr.GetDriverByName('SQLite')
print('Doing some preliminary analysis...')

try:
    # 如果文件已经存在，这里将会报错
    ds = driver.CreateDataSource(fn)
    lyr = ds.CreateLayer('layer')

# 试着保存一部分数据
except RuntimeError as e:
    # 输出错误信息，并继续
    print(e)

print('Doing some more analysis...')

2.2 新建属性字段

  通过FieldDefn构造函数提供名字和数据类型来创建一个字段。

字段数据类型	OGR常量
Integer	OFTInteger
List of integers	OFTIntegerList
Floating point number	OFTReal
String	OFTString
List of String	OFTStringList
Date	OFTDate
Date and time	OFTDateTime

  创建一个基本字段定义后，使用它添加一个字段到图层前，并且可以添加一些其他的限制，比如浮点精度或者字段宽度。
  在shapefile文件中设置字段精度，必须设置字段宽度。

import sys
import os
from osgeo import ogr
data_dir = r'E:\Google chrome\Download\GIS with python\osgeopy data\osgeopy-data-washington\osgeopy-data'

in_fn = os.path.join(data_dir, '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)

# 创建输出shp文件
driver = in_ds.GetDriver()
out_fn = os.path.join(data_dir, 'output', '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))

# 创建一个shp图层
out_lyr = out_ds.CreateLayer('precision_test',
                             in_lyr.GetSpatialRef(),
                             ogr.wkbPoint)

# 将name字段设置为宽度为6，但它将被展开。
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)

# 使用更小的精度创建两个属性字段
# 创建并添加第一个属性字段
coord_fld = ogr.FieldDefn('X_short', ogr.OFTReal)
coord_fld.SetWidth(8)      # 设置字段宽度
coord_fld.SetPrecision(3)  # 设置字段精度
# 重用FieldDefn来创建第二个字段
# 这里创建两个相同的字段定义对象
# 使创建更加容易
out_lyr.CreateField(coord_fld)
coord_fld.SetName('Y_short')
out_lyr.CreateField(coord_fld)

# 使用更小的精度创建两个属性字段。
coord_fld = ogr.FieldDefn('X_short', ogr.OFTReal)
coord_fld.SetWidth(8)
coord_fld.SetPrecision(3)
out_lyr.CreateField(coord_fld)
coord_fld.SetName('Y_short')
out_lyr.CreateField(coord_fld)

# 复制数据
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())
    out_lyr.CreateFeature(out_feat)

3. 更新数据

3.1 改变图层定义

  根据数据类型的不同，可以添加、删除、更改字段来编辑图层定义。对于一个满意的字段定义，可以使用AlterFiledDefn函数用新字段替换现有字段：
AlterFiledDefn（iField，filed_def, nFlags）

1. AlterFiledDefn是改变的属性字段对应的索引
2. filed_def是新属性字段的定义对象
3. nFlags是一个整数，是表中一个或多个常量的总和

需要更改的字段属性	OGR常量
Field name only	ALTER_NAME_FLAG
Field type only	ALTER_TYPE_FLAG
Field width and/or precision only	ALTER_WIDTH_PRECISION_FLAG
All of the above	ALTER_ALL_FLAG

  步骤 ：创建一个包含新的属性的字段定义，找到现有字段索引，确定要使用的OGR常量。如将“Name”改为“CITY_NAME”:

  首先，查看原有的数据集。

import sys
import os
import ospybook as pb
from osgeo import ogr
data_dir = r'E:\Google chrome\Download\GIS with python\osgeopy data\osgeopy-data-washington\osgeopy-data'

original_fn = os.path.join(data_dir, 'Washington', 'large_cities.shp')
new_fn = os.path.join(data_dir, 'Washington', 'large_cities2.shp')

# 复制.shp文件
pb.copy_datasource(original_fn, new_fn)

# 打开复制的shp文件准备写入
ds = ogr.Open(new_fn, 1)
if ds is None:
    sys.exit('Could not open {0}.'.format(new_fn))
lyr = ds.GetLayer(0)

# 在做改变之前，先查看数据集
print('Original attributes')
pb.print_attributes(lyr, geom=False)

  然后，改变字段属性

# 获取图层中Name的所有相关属性信息
# 创建新的文件
# 定义并使用改变新的文件属性
i = lyr.GetLayerDefn().GetFieldIndex('Name')
fld_defn = ogr.FieldDefn('City_Name', ogr.OFTString)
lyr.AlterFieldDefn(i, fld_defn, ogr.ALTER_NAME_FLAG)

# 改变多个字段的多个属性信息
# 改变POINT_X为X_coord，并设置精度为4
# 创建新的字段定义时，使用的是原始字段宽度。如果设置宽度不是足够大，结果将错误
# 所以需要使用与原始数据一样的字段宽度，这样会比较安全
# 所以字段必须重写
lyr_defn = lyr.GetLayerDefn()
i = lyr_defn.GetFieldIndex('X')
width = lyr_defn.GetFieldDefn(i).GetWidth()
fld_defn = ogr.FieldDefn('X_coord', ogr.OFTReal)
fld_defn.SetWidth(width)
fld_defn.SetPrecision(4)
flag = ogr.ALTER_NAME_FLAG + ogr.ALTER_WIDTH_PRECISION_FLAG
lyr.AlterFieldDefn(i, fld_defn, flag)

3.2 添加、更新和删除要素

  添加要素：首先基于字段创建一个空要素，填充它，然后，把它插入到该图层中去。

# 为图层中的每一个要素添加一个唯一的ID
lyr.ResetReading()
lyr.CreateField(ogr.FieldDefn('ID', ogr.OFTInteger))
n = 1
for feat in lyr:
    feat.SetField('ID', n)
    lyr.SetFeature(feat)
    n += 1
print('\nID has been added and precision has taken effect')
pb.print_attributes(lyr, geom=False)

  删除要素：

lyr.ResetReading()
for feat in lyr:
    if feat.GetField('City_Name') == 'Seattle':
        lyr.DeleteFeature(feat.GetFID())
print('\nSeattle deleted')
pb.print_attributes(lyr, geom=False)

  回收空间：

ds.ExecuteSQL('REPACK ' + lyr.GetName())
ds.ExecuteSQL('RECOMPUTE EXTENT ON ' + lyr.GetName())
print('\nDatabase packed')
pb.print_attributes(lyr, geom=False)

4. 例子

  1.使用默认精度创建json文件

import sys
import os
from osgeo import ogr
data_dir = r'E:\Google chrome\Download'

shp_fn = os.path.join(data_dir, 'global', 'ne_50m_admin_0_countries.shp')
shp_ds = ogr.Open(shp_fn, 0)
if shp_ds is None:
    sys.exit('Could not open {0}'.format(shp_fn))
shp_lyr = shp_ds.GetLayer(0)
json_driver = ogr.GetDriverByName('GeoJSON')


# 使用默认精度创建一个json文件
# 创建数据源
json_fn = os.path.join(data_dir, 'GIS with python', 'africa-default.geojson')
if os.path.exists(json_fn):
    json_driver.DeleteDataSource(json_fn)
json_ds = json_driver.CreateDataSource(json_fn)
if json_ds is None:
    sys.exit('Could not create {0}.'.format(json_fn))

# 创建无选择项的图层
json_lyr = json_ds.CreateLayer('africa',
                               shp_lyr.GetSpatialRef(),
                               ogr.wkbMultiPolygon)

# 写入数据
shp_lyr.ResetReading()
json_feat = ogr.Feature(json_lyr.GetLayerDefn())
for shp_feat in shp_lyr:
    if shp_feat.GetField('CONTINENT') == 'Africa':
        json_feat.SetGeometry(shp_feat.geometry())
        json_lyr.CreateFeature(json_feat)
del json_ds

  结果：

  2.使用精度为6，创建json文件

import sys
import os
from osgeo import ogr
data_dir = r'E:\Google chrome\Download'

shp_fn = os.path.join(data_dir, 'global', 'ne_50m_admin_0_countries.shp')
shp_ds = ogr.Open(shp_fn, 0)
if shp_ds is None:
    sys.exit('Could not open {0}'.format(shp_fn))
shp_lyr = shp_ds.GetLayer(0)
json_driver = ogr.GetDriverByName('GeoJSON')


# 使用精度为6，创建一个json文件
# 用选择项COORDINATE_PRECISION创建一个json文件，并且设置精度为为6
# 创建数据源
json_fn = os.path.join(data_dir, 'GIS with python', 'africa-6digit.geojson')
if os.path.exists(json_fn):
    json_driver.DeleteDataSource(json_fn)
json_ds = json_driver.CreateDataSource(json_fn)
if json_ds is None:
    sys.exit('Could not create {0}.'.format(json_fn))

lyr_options = ['COORDINATE_PRECISION=6']
json_lyr = json_ds.CreateLayer('africa',
                               shp_lyr.GetSpatialRef(),
                               ogr.wkbMultiPolygon,
                               lyr_options)

# 写入数据
shp_lyr.ResetReading()
json_feat = ogr.Feature(json_lyr.GetLayerDefn())
for shp_feat in shp_lyr:
    if shp_feat.GetField('CONTINENT') == 'Africa':
        json_feat.SetGeometry(shp_feat.geometry())
        json_lyr.CreateFeature(json_feat)
json_ds

  结果：

  3.Bounding box

import sys
import os
from osgeo import ogr
data_dir = r'E:\Google chrome\Download'

shp_fn = os.path.join(data_dir, 'global', 'ne_50m_admin_0_countries.shp')
shp_ds = ogr.Open(shp_fn, 0)
if shp_ds is None:
    sys.exit('Could not open {0}'.format(shp_fn))
shp_lyr = shp_ds.GetLayer(0)
json_driver = ogr.GetDriverByName('GeoJSON')

# 使用选择项COORDINATE_PRECISION 和 WRITE_BBOX 创建一个json文件
json_fn = os.path.join(data_dir, 'GIS with python', 'africa-bbox.geojson')
if os.path.exists(json_fn):
    json_driver.DeleteDataSource(json_fn)
json_ds = json_driver.CreateDataSource(json_fn)
if json_ds is None:
    sys.exit('Could not create {0}.'.format(json_fn))

lyr_options = ['COORDINATE_PRECISION=6', 'WRITE_BBOX=YES']
json_lyr = json_ds.CreateLayer('africa',
                               shp_lyr.GetSpatialRef(),
                               ogr.wkbMultiPolygon,
                               lyr_options)

# Write some data.
shp_lyr.ResetReading()
json_feat = ogr.Feature(json_lyr.GetLayerDefn())
for shp_feat in shp_lyr:
    if shp_feat.GetField('CONTINENT') == 'Africa':
        json_feat.SetGeometry(shp_feat.geometry())
        json_lyr.CreateFeature(json_feat)
del json_ds

  结果：

  4.创建一个新的例子

import sys
import os
from osgeo import ogr
data_dir = r'E:\Google chrome\Download'

shp_fn = os.path.join(data_dir, 'global', 'ne_50m_admin_0_countries.shp')
json_fn = os.path.join(data_dir, 'GIS with python', 'africa.geojson')

# 打开输入文件
shp_ds = ogr.Open(shp_fn, 0)
if shp_ds is None:
    sys.exit('Could not open {0}'.format(shp_fn))
shp_lyr = shp_ds.GetLayer(0)

# 创建输出文件
json_driver = ogr.GetDriverByName('GeoJSON')
if os.path.exists(json_fn):
    json_driver.DeleteDataSource(json_fn)
json_ds = json_driver.CreateDataSource(json_fn)
if json_ds is None:
    sys.exit('Could not create {0}.'.format(json_fn))
lyr_options = ['COORDINATE_PRECISION=6']
json_lyr = json_ds.CreateLayer('africa',
                               shp_lyr.GetSpatialRef(),
                               ogr.wkbMultiPolygon,
                               lyr_options)

# 在图层中添加属性文件
name_fld = ogr.FieldDefn('Name', ogr.OFTString)
json_lyr.CreateField(name_fld)
pop_fld = ogr.FieldDefn('Population', ogr.OFTInteger)
json_lyr.CreateField(pop_fld)


# 添加整型文件，及错误的添加字符串
test_fld = ogr.FieldDefn('Test_field', ogr.OFTInteger)
json_lyr.CreateField(test_fld)

# 获取属性定义
feat_defn = json_lyr.GetLayerDefn()

# 创建输出特征并反复使用
json_feat = ogr.Feature(feat_defn)

for shp_feat in shp_lyr:
    if shp_feat.GetField('CONTINENT') == 'Africa':

        # 复制属性值到文件中
        name = shp_feat.GetField('NAME')
        pop = shp_feat.GetField('POP_EST')
        json_feat.SetField('Name', name)
        json_feat.SetField('Population', pop)

        # 在整型文件中放置字符串
        json_feat.SetField('Test_field', name)

        # 赋值几何形状
        json_feat.SetGeometry(shp_feat.geometry())

        # 将数据输入到geojson文件中
        json_lyr.CreateFeature(json_feat)

del json_ds, shp_ds

  结果：