第05章 - 矢量数据处理基础

第05章：矢量数据处理基础

5.1 矢量数据概述

5.1.1 什么是矢量数据

矢量数据使用点、线、面等几何对象来表示地理要素，每个几何对象可以关联多个属性信息。与栅格数据相比，矢量数据具有精确的位置表达和清晰的边界定义。

┌─────────────────────────────────────────────────────────────┐
│                     矢量数据结构                              │
├─────────────────────────────────────────────────────────────┤
│                                                              │
│   点 (Point)          线 (LineString)      面 (Polygon)     │
│       •                    ●───●              ●───●         │
│                           /     \            /     \        │
│                          ●       ●          ●       ●       │
│                                              \     /        │
│                                               ●───●         │
│                                                              │
│   多点 (MultiPoint)   多线 (MultiLineString) 多面(MultiPolygon)│
│     • • •              ●──●  ●──●           ┌─┐  ┌─┐       │
│     • •                ●──●                 └─┘  └─┘       │
│                                                              │
└─────────────────────────────────────────────────────────────┘

5.1.2 矢量数据组成

组成部分 描述 示例
数据源 存储矢量数据的容器 Shapefile、GeoJSON文件、PostGIS数据库
图层 包含同类要素的集合 道路图层、建筑物图层
要素 单个地理对象 一条河流、一栋建筑
几何 要素的空间形状 点坐标、线段、多边形
属性 要素的描述信息 名称、面积、类型

5.1.3 常见矢量格式

格式 扩展名 特点 适用场景
Shapefile .shp 行业标准，兼容性最好 通用交换
GeoJSON .json/.geojson 文本格式，Web友好 Web应用
GeoPackage .gpkg 现代开放格式，SQLite基础 移动/桌面应用
FileGDB .gdb Esri格式，功能丰富 ArcGIS生态
KML/KMZ .kml/.kmz Google Earth格式 可视化展示
PostGIS - 空间数据库 企业级应用
GML .gml OGC标准XML格式 数据交换
FlatGeobuf .fgb 高性能二进制格式 大数据集

5.2 读取矢量数据

5.2.1 打开数据源

from osgeo import ogr

# 启用异常
ogr.UseExceptions()

# 打开数据源（只读）
ds = ogr.Open('data.shp')

# 打开数据源（读写）
ds = ogr.Open('data.shp', 1)  # 1 表示可写

# 检查是否成功打开
if ds is None:
    print("无法打开数据源")
else:
    print(f"驱动: {ds.GetDriver().GetName()}")
    print(f"图层数: {ds.GetLayerCount()}")

# 关闭数据源
ds = None

5.2.2 访问图层

from osgeo import ogr

ds = ogr.Open('data.shp')

# 获取图层 - 方法1：按索引
layer = ds.GetLayer(0)

# 获取图层 - 方法2：按名称
layer = ds.GetLayerByName('layer_name')

# 图层基本信息
print(f"图层名称: {layer.GetName()}")
print(f"要素数量: {layer.GetFeatureCount()}")
print(f"几何类型: {ogr.GeometryTypeToName(layer.GetGeomType())}")

# 获取图层范围
extent = layer.GetExtent()
print(f"范围: X({extent[0]:.4f}, {extent[1]:.4f}), Y({extent[2]:.4f}, {extent[3]:.4f})")

# 获取空间参考
srs = layer.GetSpatialRef()
if srs:
    print(f"坐标系: {srs.GetName()}")

ds = None

5.2.3 获取图层定义

from osgeo import ogr

ds = ogr.Open('data.shp')
layer = ds.GetLayer()

# 获取图层定义
layer_defn = layer.GetLayerDefn()

# 获取字段信息
field_count = layer_defn.GetFieldCount()
print(f"字段数量: {field_count}")

print("\n字段列表:")
for i in range(field_count):
    field_defn = layer_defn.GetFieldDefn(i)
    field_name = field_defn.GetName()
    field_type = field_defn.GetTypeName()
    field_width = field_defn.GetWidth()
    field_precision = field_defn.GetPrecision()
    
    print(f"  {i+1}. {field_name}: {field_type}({field_width},{field_precision})")

# 获取几何字段信息
geom_field_count = layer_defn.GetGeomFieldCount()
for i in range(geom_field_count):
    geom_field_defn = layer_defn.GetGeomFieldDefn(i)
    print(f"\n几何字段: {geom_field_defn.GetName()}")
    print(f"几何类型: {ogr.GeometryTypeToName(geom_field_defn.GetType())}")

ds = None

5.2.4 遍历要素

from osgeo import ogr

ds = ogr.Open('data.shp')
layer = ds.GetLayer()

# 方法1：使用 GetNextFeature 遍历
layer.ResetReading()  # 重置读取位置
feature = layer.GetNextFeature()
while feature:
    fid = feature.GetFID()
    geom = feature.GetGeometryRef()
    
    print(f"FID: {fid}, 几何类型: {geom.GetGeometryName()}")
    
    feature = layer.GetNextFeature()

# 方法2：使用 for 循环遍历（Python 风格）
layer.ResetReading()
for feature in layer:
    print(f"FID: {feature.GetFID()}")

# 方法3：按 FID 获取特定要素
feature = layer.GetFeature(5)  # 获取 FID=5 的要素

# 方法4：使用索引获取
layer.SetNextByIndex(10)  # 设置下一个要读取的位置
feature = layer.GetNextFeature()

ds = None

5.2.5 读取属性值

from osgeo import ogr

ds = ogr.Open('data.shp')
layer = ds.GetLayer()

for feature in layer:
    # 获取 FID
    fid = feature.GetFID()
    
    # 方法1：按字段名获取
    name = feature.GetField('name')
    population = feature.GetField('population')
    area = feature.GetField('area')
    
    # 方法2：按字段索引获取
    value = feature.GetField(0)
    
    # 方法3：获取特定类型
    int_value = feature.GetFieldAsInteger('count')
    float_value = feature.GetFieldAsDouble('area')
    str_value = feature.GetFieldAsString('name')
    
    # 方法4：获取所有字段为字典
    fields = {}
    layer_defn = layer.GetLayerDefn()
    for i in range(layer_defn.GetFieldCount()):
        field_name = layer_defn.GetFieldDefn(i).GetName()
        fields[field_name] = feature.GetField(i)
    
    # 检查字段是否为空
    if feature.IsFieldNull('name'):
        print("name 字段为空")
    
    print(f"FID={fid}: name={name}, population={population}")

ds = None

5.2.6 读取几何

from osgeo import ogr

ds = ogr.Open('data.shp')
layer = ds.GetLayer()

for feature in layer:
    geom = feature.GetGeometryRef()
    
    if geom is None:
        print("几何为空")
        continue
    
    # 几何类型
    geom_type = geom.GetGeometryType()
    geom_name = geom.GetGeometryName()
    print(f"几何类型: {geom_name}")
    
    # 导出为各种格式
    wkt = geom.ExportToWkt()
    geojson = geom.ExportToJson()
    wkb = geom.ExportToWkb()
    
    # 获取坐标
    if geom_type == ogr.wkbPoint:
        x, y = geom.GetX(), geom.GetY()
        print(f"点坐标: ({x}, {y})")
    
    elif geom_type == ogr.wkbLineString:
        point_count = geom.GetPointCount()
        print(f"线点数: {point_count}")
        for i in range(point_count):
            x, y = geom.GetPoint(i)[:2]
            print(f"  点{i}: ({x}, {y})")
    
    elif geom_type == ogr.wkbPolygon:
        ring_count = geom.GetGeometryCount()
        print(f"环数量: {ring_count}")
        for i in range(ring_count):
            ring = geom.GetGeometryRef(i)
            point_count = ring.GetPointCount()
            print(f"  环{i}有 {point_count} 个点")
    
    # 几何属性
    print(f"面积: {geom.GetArea()}")
    print(f"长度: {geom.Length()}")
    
    # 获取范围
    envelope = geom.GetEnvelope()
    print(f"范围: {envelope}")
    
    # 获取质心
    centroid = geom.Centroid()
    print(f"质心: ({centroid.GetX()}, {centroid.GetY()})")

ds = None

5.3 过滤和查询

5.3.1 属性过滤

from osgeo import ogr

ds = ogr.Open('cities.shp')
layer = ds.GetLayer()

# 设置属性过滤器（SQL WHERE 子句语法）
layer.SetAttributeFilter("population > 1000000")

print(f"人口超过100万的城市数量: {layer.GetFeatureCount()}")

for feature in layer:
    name = feature.GetField('name')
    pop = feature.GetField('population')
    print(f"  {name}: {pop}")

# 清除过滤器
layer.SetAttributeFilter(None)

# 更复杂的过滤条件
layer.SetAttributeFilter("""
    population > 500000 AND 
    country = 'China' AND 
    name LIKE '北%'
""")

# 支持的运算符
# =, <>, <, >, <=, >=
# AND, OR, NOT
# LIKE (通配符 % _)
# IN (值列表)
# BETWEEN
# IS NULL, IS NOT NULL

ds = None

5.3.2 空间过滤

from osgeo import ogr

ds = ogr.Open('data.shp')
layer = ds.GetLayer()

# 方法1：矩形范围过滤
min_x, min_y, max_x, max_y = 116.0, 39.0, 117.0, 40.0
layer.SetSpatialFilterRect(min_x, min_y, max_x, max_y)

print(f"范围内要素数: {layer.GetFeatureCount()}")

# 方法2：几何对象过滤
# 创建过滤几何（圆形缓冲区）
center = ogr.Geometry(ogr.wkbPoint)
center.AddPoint(116.4, 39.9)
buffer = center.Buffer(0.5)  # 0.5度缓冲区

layer.SetSpatialFilter(buffer)
print(f"缓冲区内要素数: {layer.GetFeatureCount()}")

# 方法3：使用多边形过滤
filter_polygon = ogr.CreateGeometryFromWkt("""
    POLYGON((116 39, 117 39, 117 40, 116 40, 116 39))
""")
layer.SetSpatialFilter(filter_polygon)

# 清除空间过滤器
layer.SetSpatialFilter(None)

ds = None

5.3.3 SQL 查询

from osgeo import ogr

ds = ogr.Open('data.gpkg')

# 执行 SQL 查询
sql = """
    SELECT name, population, ST_Area(geom) as area
    FROM cities
    WHERE population > 1000000
    ORDER BY population DESC
    LIMIT 10
"""

result_layer = ds.ExecuteSQL(sql)

if result_layer:
    print("查询结果:")
    for feature in result_layer:
        name = feature.GetField('name')
        pop = feature.GetField('population')
        area = feature.GetField('area')
        print(f"  {name}: {pop}, {area:.2f}")
    
    # 释放结果图层
    ds.ReleaseResultSet(result_layer)

# 空间 SQL 查询（需要支持空间函数的格式，如 GeoPackage、PostGIS）
spatial_sql = """
    SELECT a.name, a.population
    FROM cities a, provinces b
    WHERE ST_Within(a.geom, b.geom) AND b.name = '北京'
"""

# 使用 SQLite 方言
ds.ExecuteSQL(sql, dialect='SQLITE')

# 使用 OGRSQL 方言（默认）
ds.ExecuteSQL(sql, dialect='OGRSQL')

ds = None

5.4 创建矢量数据

5.4.1 创建数据源和图层

from osgeo import ogr, osr

def create_shapefile(filepath, layer_name, geom_type, srs_epsg=4326, fields=None):
    """创建 Shapefile"""
    
    # 获取驱动
    driver = ogr.GetDriverByName('ESRI Shapefile')
    
    # 如果文件存在则删除
    if driver.DeleteDataSource(filepath) != 0:
        pass  # 文件不存在，忽略错误
    
    # 创建数据源
    ds = driver.CreateDataSource(filepath)
    
    # 创建空间参考
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(srs_epsg)
    
    # 创建图层
    layer = ds.CreateLayer(layer_name, srs, geom_type)
    
    # 添加字段
    if fields:
        for field_name, field_type, *extra in fields:
            field_defn = ogr.FieldDefn(field_name, field_type)
            
            if field_type == ogr.OFTString and extra:
                field_defn.SetWidth(extra[0])
            elif field_type == ogr.OFTReal and len(extra) >= 2:
                field_defn.SetWidth(extra[0])
                field_defn.SetPrecision(extra[1])
            
            layer.CreateField(field_defn)
    
    return ds, layer

# 创建点图层
fields = [
    ('name', ogr.OFTString, 50),
    ('population', ogr.OFTInteger64),
    ('area', ogr.OFTReal, 12, 2),
    ('create_date', ogr.OFTDate),
]

ds, layer = create_shapefile(
    'cities.shp', 
    'cities', 
    ogr.wkbPoint,
    srs_epsg=4326,
    fields=fields
)

# 不要忘记关闭
ds = None

5.4.2 创建 GeoJSON

from osgeo import ogr, osr

def create_geojson(filepath, layer_name, geom_type, srs_epsg=4326, fields=None):
    """创建 GeoJSON 文件"""
    
    driver = ogr.GetDriverByName('GeoJSON')
    ds = driver.CreateDataSource(filepath)
    
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(srs_epsg)
    
    layer = ds.CreateLayer(layer_name, srs, geom_type)
    
    if fields:
        for field_name, field_type, *extra in fields:
            field_defn = ogr.FieldDefn(field_name, field_type)
            layer.CreateField(field_defn)
    
    return ds, layer

# 创建示例
ds, layer = create_geojson('output.geojson', 'points', ogr.wkbPoint)
ds = None

5.4.3 创建 GeoPackage

from osgeo import ogr, osr

def create_geopackage(filepath, layer_name, geom_type, srs_epsg=4326, fields=None):
    """创建 GeoPackage 文件"""
    
    driver = ogr.GetDriverByName('GPKG')
    
    # GeoPackage 可以包含多个图层，所以检查文件是否存在
    import os
    if os.path.exists(filepath):
        ds = ogr.Open(filepath, 1)  # 打开现有文件
    else:
        ds = driver.CreateDataSource(filepath)
    
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(srs_epsg)
    
    # 创建选项
    options = ['FID=fid', 'GEOMETRY_NAME=geom']
    
    layer = ds.CreateLayer(layer_name, srs, geom_type, options=options)
    
    if fields:
        for field_name, field_type, *extra in fields:
            field_defn = ogr.FieldDefn(field_name, field_type)
            layer.CreateField(field_defn)
    
    return ds, layer

# 创建多图层 GeoPackage
ds, points_layer = create_geopackage('data.gpkg', 'points', ogr.wkbPoint)
ds, lines_layer = create_geopackage('data.gpkg', 'roads', ogr.wkbLineString)
ds, polygons_layer = create_geopackage('data.gpkg', 'buildings', ogr.wkbPolygon)
ds = None

5.4.4 添加要素

from osgeo import ogr, osr

# 创建 Shapefile
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.CreateDataSource('cities.shp')

srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

layer = ds.CreateLayer('cities', srs, ogr.wkbPoint)

# 添加字段
layer.CreateField(ogr.FieldDefn('name', ogr.OFTString))
layer.CreateField(ogr.FieldDefn('population', ogr.OFTInteger64))

# 获取图层定义
layer_defn = layer.GetLayerDefn()

# 添加要素 - 方法1：完整方式
feature = ogr.Feature(layer_defn)
feature.SetField('name', '北京')
feature.SetField('population', 21540000)

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(116.4, 39.9)
feature.SetGeometry(point)

layer.CreateFeature(feature)
feature = None

# 添加要素 - 方法2：使用字典方式设置属性
cities = [
    {'name': '上海', 'population': 24280000, 'coords': (121.5, 31.2)},
    {'name': '广州', 'population': 18680000, 'coords': (113.3, 23.1)},
    {'name': '深圳', 'population': 17560000, 'coords': (114.1, 22.5)},
]

for city in cities:
    feature = ogr.Feature(layer_defn)
    feature.SetField('name', city['name'])
    feature.SetField('population', city['population'])
    
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(*city['coords'])
    feature.SetGeometry(point)
    
    layer.CreateFeature(feature)
    feature = None

# 保存并关闭
ds = None

5.4.5 创建不同类型的几何

from osgeo import ogr

# 创建点
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(116.4, 39.9)

# 创建带高程的点
point_z = ogr.Geometry(ogr.wkbPoint25D)
point_z.AddPoint(116.4, 39.9, 100)

# 创建线
line = ogr.Geometry(ogr.wkbLineString)
line.AddPoint(116.0, 39.0)
line.AddPoint(116.5, 39.5)
line.AddPoint(117.0, 40.0)

# 创建多边形
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(116.0, 39.0)
ring.AddPoint(117.0, 39.0)
ring.AddPoint(117.0, 40.0)
ring.AddPoint(116.0, 40.0)
ring.AddPoint(116.0, 39.0)  # 闭合

polygon = ogr.Geometry(ogr.wkbPolygon)
polygon.AddGeometry(ring)

# 创建带孔的多边形
outer_ring = ogr.Geometry(ogr.wkbLinearRing)
outer_ring.AddPoint(0, 0)
outer_ring.AddPoint(10, 0)
outer_ring.AddPoint(10, 10)
outer_ring.AddPoint(0, 10)
outer_ring.AddPoint(0, 0)

inner_ring = ogr.Geometry(ogr.wkbLinearRing)
inner_ring.AddPoint(2, 2)
inner_ring.AddPoint(8, 2)
inner_ring.AddPoint(8, 8)
inner_ring.AddPoint(2, 8)
inner_ring.AddPoint(2, 2)

polygon_with_hole = ogr.Geometry(ogr.wkbPolygon)
polygon_with_hole.AddGeometry(outer_ring)
polygon_with_hole.AddGeometry(inner_ring)

# 创建多点
multipoint = ogr.Geometry(ogr.wkbMultiPoint)
for x, y in [(1, 1), (2, 2), (3, 3)]:
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(x, y)
    multipoint.AddGeometry(point)

# 创建多线
multiline = ogr.Geometry(ogr.wkbMultiLineString)
for coords in [[(0, 0), (1, 1)], [(2, 2), (3, 3)]]:
    line = ogr.Geometry(ogr.wkbLineString)
    for x, y in coords:
        line.AddPoint(x, y)
    multiline.AddGeometry(line)

# 创建多多边形
multipolygon = ogr.Geometry(ogr.wkbMultiPolygon)
multipolygon.AddGeometry(polygon)

# 从 WKT 创建
geom_from_wkt = ogr.CreateGeometryFromWkt(
    'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))'
)

# 从 GeoJSON 创建
geom_from_json = ogr.CreateGeometryFromJson(
    '{"type":"Point","coordinates":[116.4,39.9]}'
)

# 从 WKB 创建
# geom_from_wkb = ogr.CreateGeometryFromWkb(wkb_bytes)

5.5 编辑矢量数据

5.5.1 更新要素

from osgeo import ogr

# 打开数据源（可写模式）
ds = ogr.Open('data.shp', 1)
layer = ds.GetLayer()

# 更新单个要素
feature = layer.GetFeature(5)  # 获取 FID=5 的要素
if feature:
    # 更新属性
    feature.SetField('name', '新名称')
    feature.SetField('population', 1000000)
    
    # 更新几何
    new_geom = ogr.CreateGeometryFromWkt('POINT(116.5 40.0)')
    feature.SetGeometry(new_geom)
    
    # 保存更新
    layer.SetFeature(feature)

# 批量更新
layer.ResetReading()
for feature in layer:
    old_value = feature.GetField('population')
    new_value = int(old_value * 1.1)  # 增加10%
    feature.SetField('population', new_value)
    layer.SetFeature(feature)

ds = None

5.5.2 删除要素

from osgeo import ogr

ds = ogr.Open('data.shp', 1)
layer = ds.GetLayer()

# 删除单个要素
layer.DeleteFeature(5)  # 删除 FID=5 的要素

# 条件删除
layer.SetAttributeFilter("population < 10000")
for feature in layer:
    layer.DeleteFeature(feature.GetFID())

# 清除过滤器
layer.SetAttributeFilter(None)

# 使用 SQL 删除（某些格式支持）
# ds.ExecuteSQL("DELETE FROM layer_name WHERE condition")

ds = None

# 注意：Shapefile 删除要素后会留下空洞
# 需要重新打包来真正删除
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.Open('data.shp', 1)
ds.ExecuteSQL('REPACK data')  # 重新打包
ds = None

5.5.3 添加和删除字段

from osgeo import ogr

ds = ogr.Open('data.shp', 1)
layer = ds.GetLayer()

# 添加新字段
new_field = ogr.FieldDefn('new_column', ogr.OFTString)
new_field.SetWidth(50)
layer.CreateField(new_field)

# 删除字段（某些格式支持）
layer_defn = layer.GetLayerDefn()
field_index = layer_defn.GetFieldIndex('old_column')
if field_index >= 0:
    layer.DeleteField(field_index)

# 修改字段（某些格式支持）
# layer.AlterFieldDefn(field_index, new_field_defn, ogr.ALTER_NAME_FLAG)

ds = None

5.6 几何操作

5.6.1 空间关系判断

from osgeo import ogr

# 创建两个几何对象
geom1 = ogr.CreateGeometryFromWkt('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))')
geom2 = ogr.CreateGeometryFromWkt('POLYGON((5 5, 15 5, 15 15, 5 15, 5 5))')
point = ogr.CreateGeometryFromWkt('POINT(5 5)')

# 空间关系判断
print(f"geom1 包含 point: {geom1.Contains(point)}")
print(f"point 在 geom1 内: {point.Within(geom1)}")
print(f"geom1 与 geom2 相交: {geom1.Intersects(geom2)}")
print(f"geom1 与 geom2 相离: {geom1.Disjoint(geom2)}")
print(f"geom1 与 geom2 接触: {geom1.Touches(geom2)}")
print(f"geom1 与 geom2 交叉: {geom1.Crosses(geom2)}")
print(f"geom1 与 geom2 重叠: {geom1.Overlaps(geom2)}")
print(f"geom1 等于 geom2: {geom1.Equals(geom2)}")

# 空间关系详情（DE-9IM）
relation = geom1.Relate(geom2)
print(f"空间关系矩阵: {relation}")

5.6.2 空间操作

from osgeo import ogr

geom1 = ogr.CreateGeometryFromWkt('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))')
geom2 = ogr.CreateGeometryFromWkt('POLYGON((5 5, 15 5, 15 15, 5 15, 5 5))')

# 缓冲区
buffer = geom1.Buffer(1.0)
print(f"缓冲区面积: {buffer.GetArea()}")

# 负缓冲（收缩）
negative_buffer = geom1.Buffer(-1.0)

# 交集
intersection = geom1.Intersection(geom2)
print(f"交集: {intersection.ExportToWkt()}")

# 并集
union = geom1.Union(geom2)
print(f"并集面积: {union.GetArea()}")

# 差集
difference = geom1.Difference(geom2)
print(f"差集: {difference.ExportToWkt()}")

# 对称差
sym_difference = geom1.SymDifference(geom2)
print(f"对称差: {sym_difference.ExportToWkt()}")

# 凸包
convex_hull = geom1.ConvexHull()

# 简化
simplified = geom1.Simplify(0.1)
simplified_preserve = geom1.SimplifyPreserveTopology(0.1)

# 质心
centroid = geom1.Centroid()
print(f"质心: {centroid.ExportToWkt()}")

# 外包矩形
envelope = geom1.GetEnvelope()
print(f"外包矩形: MinX={envelope[0]}, MaxX={envelope[1]}, MinY={envelope[2]}, MaxY={envelope[3]}")

5.6.3 几何属性计算

from osgeo import ogr

polygon = ogr.CreateGeometryFromWkt(
    'POLYGON((0 0, 100 0, 100 100, 0 100, 0 0))'
)
line = ogr.CreateGeometryFromWkt(
    'LINESTRING(0 0, 100 0, 100 100)'
)

# 面积（平面坐标系下的单位面积）
area = polygon.GetArea()
print(f"面积: {area}")

# 长度/周长
perimeter = polygon.Length()  # 多边形周长
length = line.Length()        # 线长度
print(f"周长: {perimeter}, 线长度: {length}")

# 几何验证
is_valid = polygon.IsValid()
is_simple = polygon.IsSimple()
is_ring = polygon.IsRing()
is_empty = polygon.IsEmpty()

print(f"有效: {is_valid}, 简单: {is_simple}, 闭合: {is_ring}, 空: {is_empty}")

# 点数量
point_count = polygon.GetGeometryRef(0).GetPointCount()  # 外环点数
print(f"外环点数: {point_count}")

# 几何维度
dim = polygon.GetDimension()  # 0=点, 1=线, 2=面
coord_dim = polygon.GetCoordinateDimension()  # 2 或 3
print(f"几何维度: {dim}, 坐标维度: {coord_dim}")

5.6.4 坐标转换

from osgeo import ogr, osr

# 创建几何
geom = ogr.CreateGeometryFromWkt('POINT(116.4 39.9)')

# 创建坐标转换
src_srs = osr.SpatialReference()
src_srs.ImportFromEPSG(4326)  # WGS84

dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(3857)  # Web Mercator

# 创建转换器
transform = osr.CoordinateTransformation(src_srs, dst_srs)

# 转换几何（就地修改）
geom.Transform(transform)
print(f"转换后: {geom.ExportToWkt()}")

# 克隆并转换
original = ogr.CreateGeometryFromWkt('POINT(116.4 39.9)')
transformed = original.Clone()
transformed.Transform(transform)
print(f"原始: {original.ExportToWkt()}")
print(f"转换后: {transformed.ExportToWkt()}")

5.7 格式转换

5.7.1 使用 ogr2ogr 函数

from osgeo import gdal, ogr

def convert_vector(src_path, dst_path, dst_format='ESRI Shapefile', options=None):
    """转换矢量格式"""
    
    if options is None:
        options = []
    
    vector_options = gdal.VectorTranslateOptions(
        format=dst_format,
        layerCreationOptions=options
    )
    
    ds = gdal.VectorTranslate(dst_path, src_path, options=vector_options)
    ds = None

# Shapefile 转 GeoJSON
convert_vector('input.shp', 'output.geojson', 'GeoJSON')

# Shapefile 转 GeoPackage
convert_vector('input.shp', 'output.gpkg', 'GPKG')

# GeoJSON 转 Shapefile
convert_vector('input.geojson', 'output.shp', 'ESRI Shapefile')

5.7.2 手动复制转换

from osgeo import ogr, osr

def copy_layer(src_path, dst_path, dst_format='ESRI Shapefile'):
    """复制图层到新格式"""
    
    # 打开源数据
    src_ds = ogr.Open(src_path)
    src_layer = src_ds.GetLayer()
    
    # 创建目标数据源
    driver = ogr.GetDriverByName(dst_format)
    
    # 删除已存在的目标文件
    if driver.DeleteDataSource(dst_path) != 0:
        pass
    
    dst_ds = driver.CreateDataSource(dst_path)
    
    # 复制图层定义
    src_layer_defn = src_layer.GetLayerDefn()
    srs = src_layer.GetSpatialRef()
    geom_type = src_layer.GetGeomType()
    
    dst_layer = dst_ds.CreateLayer(
        src_layer.GetName(),
        srs,
        geom_type
    )
    
    # 复制字段
    for i in range(src_layer_defn.GetFieldCount()):
        field_defn = src_layer_defn.GetFieldDefn(i)
        dst_layer.CreateField(field_defn)
    
    # 复制要素
    dst_layer_defn = dst_layer.GetLayerDefn()
    
    for src_feature in src_layer:
        dst_feature = ogr.Feature(dst_layer_defn)
        
        # 复制几何
        geom = src_feature.GetGeometryRef()
        if geom:
            dst_feature.SetGeometry(geom.Clone())
        
        # 复制属性
        for i in range(dst_layer_defn.GetFieldCount()):
            dst_feature.SetField(i, src_feature.GetField(i))
        
        dst_layer.CreateFeature(dst_feature)
        dst_feature = None
    
    src_ds = None
    dst_ds = None

# 使用示例
copy_layer('input.shp', 'output.geojson', 'GeoJSON')

5.7.3 转换时进行投影

from osgeo import gdal

def convert_with_reprojection(src_path, dst_path, dst_srs='EPSG:3857', dst_format='GeoJSON'):
    """转换格式并重投影"""
    
    options = gdal.VectorTranslateOptions(
        format=dst_format,
        dstSRS=dst_srs,
        reproject=True
    )
    
    ds = gdal.VectorTranslate(dst_path, src_path, options=options)
    ds = None

# WGS84 转 Web Mercator
convert_with_reprojection('wgs84.shp', 'webmercator.geojson', 'EPSG:3857')

5.7.4 转换时进行过滤

from osgeo import gdal

def convert_with_filter(src_path, dst_path, where_clause=None, 
                        bbox=None, dst_format='GeoJSON'):
    """转换格式并过滤"""
    
    options_dict = {
        'format': dst_format,
    }
    
    if where_clause:
        options_dict['where'] = where_clause
    
    if bbox:
        options_dict['spatFilter'] = bbox  # [minx, miny, maxx, maxy]
    
    options = gdal.VectorTranslateOptions(**options_dict)
    
    ds = gdal.VectorTranslate(dst_path, src_path, options=options)
    ds = None

# 按属性过滤
convert_with_filter(
    'cities.shp', 
    'big_cities.geojson',
    where_clause="population > 1000000"
)

# 按范围过滤
convert_with_filter(
    'cities.shp',
    'beijing_cities.geojson',
    bbox=[116.0, 39.0, 117.0, 40.0]
)

5.8 数据库操作

5.8.1 连接 PostGIS

from osgeo import ogr

# 连接 PostGIS 数据库
connection_string = "PG:host=localhost dbname=gisdb user=postgres password=secret"
ds = ogr.Open(connection_string)

if ds:
    print(f"连接成功，图层数: {ds.GetLayerCount()}")
    
    # 列出所有图层
    for i in range(ds.GetLayerCount()):
        layer = ds.GetLayer(i)
        print(f"  {layer.GetName()}: {layer.GetFeatureCount()} 个要素")
    
    ds = None
else:
    print("连接失败")

# 连接时只加载特定图层
connection_string = "PG:host=localhost dbname=gisdb user=postgres password=secret tables=cities,roads"
ds = ogr.Open(connection_string)

5.8.2 查询 PostGIS

from osgeo import ogr

ds = ogr.Open("PG:host=localhost dbname=gisdb user=postgres password=secret")

# 执行 SQL 查询
sql = """
    SELECT name, population, ST_AsText(geom) as geom_wkt
    FROM cities
    WHERE population > 1000000
"""

result = ds.ExecuteSQL(sql)

for feature in result:
    print(f"{feature.GetField('name')}: {feature.GetField('population')}")

ds.ReleaseResultSet(result)

# 空间查询
spatial_sql = """
    SELECT a.name, a.population
    FROM cities a, provinces b
    WHERE ST_Within(a.geom, b.geom)
    AND b.name = '广东省'
"""

result = ds.ExecuteSQL(spatial_sql)
# ... 处理结果

ds = None

5.8.3 写入 PostGIS

from osgeo import ogr, osr

# 连接数据库（可写）
ds = ogr.Open("PG:host=localhost dbname=gisdb user=postgres password=secret", 1)

# 创建新表
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

# 创建图层（表）
layer = ds.CreateLayer(
    'new_cities',
    srs,
    ogr.wkbPoint,
    ['GEOMETRY_NAME=geom', 'OVERWRITE=YES']
)

# 添加字段
layer.CreateField(ogr.FieldDefn('name', ogr.OFTString))
layer.CreateField(ogr.FieldDefn('population', ogr.OFTInteger64))

# 添加要素
layer_defn = layer.GetLayerDefn()

feature = ogr.Feature(layer_defn)
feature.SetField('name', '测试城市')
feature.SetField('population', 1000000)

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(116.4, 39.9)
feature.SetGeometry(point)

layer.CreateFeature(feature)

# 开始事务（提高批量插入性能）
layer.StartTransaction()

for i in range(1000):
    feature = ogr.Feature(layer_defn)
    feature.SetField('name', f'城市{i}')
    feature.SetField('population', i * 1000)
    
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(116.0 + i * 0.01, 39.0 + i * 0.01)
    feature.SetGeometry(point)
    
    layer.CreateFeature(feature)

layer.CommitTransaction()

ds = None

5.9 本章小结

本章详细介绍了 OGR 矢量数据处理的基础知识：

  1. 矢量数据概述：数据结构、组成和常见格式
  2. 读取矢量数据：
    • 打开数据源和访问图层
    • 遍历要素
    • 读取属性和几何
  3. 过滤和查询：
    • 属性过滤
    • 空间过滤
    • SQL 查询
  4. 创建矢量数据：
    • 创建不同格式的数据源
    • 定义字段和添加要素
    • 创建各种几何类型
  5. 编辑矢量数据：
    • 更新和删除要素
    • 修改字段定义
  6. 几何操作：
    • 空间关系判断
    • 空间运算（缓冲区、交集、并集等）
    • 几何属性计算
  7. 格式转换：
    • 使用 VectorTranslate
    • 手动复制
    • 带投影和过滤的转换
  8. 数据库操作：
    • PostGIS 连接和查询
    • 数据写入和事务

5.10 思考与练习

  1. 解释矢量数据中数据源、图层、要素、几何的层次关系。
  2. 如何高效遍历一个包含百万要素的 Shapefile？
  3. 比较 Shapefile 和 GeoPackage 格式的优缺点。
  4. 编写代码计算两个多边形图层的空间交集。
  5. 实现一个函数，将 PostGIS 中的数据导出为 GeoJSON。
  6. 如何在不加载所有数据的情况下获取图层的要素数量？
  7. 编写代码为多边形图层中的每个要素计算面积并存储到属性字段中。
posted @ 2025-12-29 11:40  我才是银古  阅读(11)  评论(0)    收藏  举报