第12章 - 矢量数据高级处理
第12章：矢量数据高级处理
12.1 空间分析
12.1.1 缓冲区分析
from osgeo import ogr, osr
def buffer_analysis(input_path, output_path, distance, dissolve=False):
"""缓冲区分析"""
ogr.UseExceptions()
src_ds = ogr.Open(input_path)
src_layer = src_ds.GetLayer()
# 创建输出
driver = ogr.GetDriverByName('ESRI Shapefile')
if driver.DeleteDataSource(output_path) != 0:
pass
dst_ds = driver.CreateDataSource(output_path)
dst_layer = dst_ds.CreateLayer(
'buffer',
src_layer.GetSpatialRef(),
ogr.wkbPolygon
)
# 复制字段
src_defn = src_layer.GetLayerDefn()
for i in range(src_defn.GetFieldCount()):
dst_layer.CreateField(src_defn.GetFieldDefn(i))
dst_defn = dst_layer.GetLayerDefn()
if dissolve:
# 溶解所有缓冲区
union_geom = None
for feature in src_layer:
geom = feature.GetGeometryRef()
if geom:
buffer_geom = geom.Buffer(distance)
if union_geom is None:
union_geom = buffer_geom
else:
union_geom = union_geom.Union(buffer_geom)
if union_geom:
out_feature = ogr.Feature(dst_defn)
out_feature.SetGeometry(union_geom)
dst_layer.CreateFeature(out_feature)
else:
# 单独的缓冲区
for feature in src_layer:
geom = feature.GetGeometryRef()
if geom:
buffer_geom = geom.Buffer(distance)
out_feature = ogr.Feature(dst_defn)
out_feature.SetGeometry(buffer_geom)
for i in range(src_defn.GetFieldCount()):
out_feature.SetField(i, feature.GetField(i))
dst_layer.CreateFeature(out_feature)
src_ds = None
dst_ds = None
12.1.2 叠加分析
from osgeo import ogr
def intersection_analysis(layer1_path, layer2_path, output_path):
"""交集分析"""
ogr.UseExceptions()
ds1 = ogr.Open(layer1_path)
layer1 = ds1.GetLayer()
ds2 = ogr.Open(layer2_path)
layer2 = ds2.GetLayer()
# 创建输出
driver = ogr.GetDriverByName('ESRI Shapefile')
out_ds = driver.CreateDataSource(output_path)
out_layer = out_ds.CreateLayer(
'intersection',
layer1.GetSpatialRef(),
ogr.wkbPolygon
)
# 添加两个图层的字段
for layer in [layer1, layer2]:
defn = layer.GetLayerDefn()
for i in range(defn.GetFieldCount()):
field_defn = defn.GetFieldDefn(i)
try:
out_layer.CreateField(field_defn)
except:
pass # 字段已存在
out_defn = out_layer.GetLayerDefn()
# 执行交集
for feat1 in layer1:
geom1 = feat1.GetGeometryRef()
if not geom1:
continue
layer2.SetSpatialFilter(geom1)
for feat2 in layer2:
geom2 = feat2.GetGeometryRef()
if not geom2:
continue
if geom1.Intersects(geom2):
intersection = geom1.Intersection(geom2)
if not intersection.IsEmpty():
out_feat = ogr.Feature(out_defn)
out_feat.SetGeometry(intersection)
# 复制属性
for i in range(feat1.GetFieldCount()):
out_feat.SetField(i, feat1.GetField(i))
out_layer.CreateFeature(out_feat)
layer2.ResetReading()
ds1 = None
ds2 = None
out_ds = None
def union_analysis(layer1_path, layer2_path, output_path):
"""并集分析"""
ogr.UseExceptions()
ds1 = ogr.Open(layer1_path)
layer1 = ds1.GetLayer()
ds2 = ogr.Open(layer2_path)
layer2 = ds2.GetLayer()
# 创建输出
driver = ogr.GetDriverByName('ESRI Shapefile')
out_ds = driver.CreateDataSource(output_path)
out_layer = out_ds.CreateLayer(
'union',
layer1.GetSpatialRef(),
ogr.wkbPolygon
)
out_defn = out_layer.GetLayerDefn()
# 收集所有几何
union_geom = None
for layer in [layer1, layer2]:
for feature in layer:
geom = feature.GetGeometryRef()
if geom:
if union_geom is None:
union_geom = geom.Clone()
else:
union_geom = union_geom.Union(geom)
layer.ResetReading()
if union_geom:
out_feat = ogr.Feature(out_defn)
out_feat.SetGeometry(union_geom)
out_layer.CreateFeature(out_feat)
ds1 = None
ds2 = None
out_ds = None
def difference_analysis(layer1_path, layer2_path, output_path):
"""差集分析"""
ogr.UseExceptions()
ds1 = ogr.Open(layer1_path)
layer1 = ds1.GetLayer()
ds2 = ogr.Open(layer2_path)
layer2 = ds2.GetLayer()
# 合并 layer2 的所有几何
union2 = None
for feature in layer2:
geom = feature.GetGeometryRef()
if geom:
if union2 is None:
union2 = geom.Clone()
else:
union2 = union2.Union(geom)
layer2.ResetReading()
# 创建输出
driver = ogr.GetDriverByName('ESRI Shapefile')
out_ds = driver.CreateDataSource(output_path)
out_layer = out_ds.CreateLayer(
'difference',
layer1.GetSpatialRef(),
layer1.GetGeomType()
)
# 复制字段
src_defn = layer1.GetLayerDefn()
for i in range(src_defn.GetFieldCount()):
out_layer.CreateField(src_defn.GetFieldDefn(i))
out_defn = out_layer.GetLayerDefn()
# 执行差集
for feature in layer1:
geom = feature.GetGeometryRef()
if geom and union2:
diff = geom.Difference(union2)
if diff and not diff.IsEmpty():
out_feat = ogr.Feature(out_defn)
out_feat.SetGeometry(diff)
for i in range(src_defn.GetFieldCount()):
out_feat.SetField(i, feature.GetField(i))
out_layer.CreateFeature(out_feat)
ds1 = None
ds2 = None
out_ds = None
12.1.3 空间连接
from osgeo import ogr
def spatial_join(target_path, join_path, output_path, join_type='intersects'):
"""空间连接"""
ogr.UseExceptions()
target_ds = ogr.Open(target_path)
target_layer = target_ds.GetLayer()
join_ds = ogr.Open(join_path)
join_layer = join_ds.GetLayer()
# 创建输出
driver = ogr.GetDriverByName('ESRI Shapefile')
out_ds = driver.CreateDataSource(output_path)
out_layer = out_ds.CreateLayer(
'joined',
target_layer.GetSpatialRef(),
target_layer.GetGeomType()
)
# 复制目标图层字段
target_defn = target_layer.GetLayerDefn()
for i in range(target_defn.GetFieldCount()):
out_layer.CreateField(target_defn.GetFieldDefn(i))
# 添加连接图层字段（加前缀）
join_defn = join_layer.GetLayerDefn()
for i in range(join_defn.GetFieldCount()):
field = join_defn.GetFieldDefn(i)
new_field = ogr.FieldDefn(f'join_{field.GetName()}', field.GetType())
new_field.SetWidth(field.GetWidth())
out_layer.CreateField(new_field)
out_defn = out_layer.GetLayerDefn()
# 空间关系判断方法
spatial_methods = {
'intersects': lambda g1, g2: g1.Intersects(g2),
'contains': lambda g1, g2: g1.Contains(g2),
'within': lambda g1, g2: g1.Within(g2),
'touches': lambda g1, g2: g1.Touches(g2),
}
check_spatial = spatial_methods.get(join_type, spatial_methods['intersects'])
# 执行连接
for target_feat in target_layer:
target_geom = target_feat.GetGeometryRef()
if not target_geom:
continue
join_layer.SetSpatialFilter(target_geom)
for join_feat in join_layer:
join_geom = join_feat.GetGeometryRef()
if join_geom and check_spatial(target_geom, join_geom):
out_feat = ogr.Feature(out_defn)
out_feat.SetGeometry(target_geom.Clone())
# 复制目标属性
for i in range(target_defn.GetFieldCount()):
out_feat.SetField(i, target_feat.GetField(i))
# 复制连接属性
for i in range(join_defn.GetFieldCount()):
out_feat.SetField(
target_defn.GetFieldCount() + i,
join_feat.GetField(i)
)
out_layer.CreateFeature(out_feat)
break # 只取第一个匹配
join_layer.ResetReading()
target_ds = None
join_ds = None
out_ds = None
12.2 几何操作
12.2.1 几何简化
from osgeo import ogr
def simplify_layer(input_path, output_path, tolerance, preserve_topology=True):
"""简化几何"""
ogr.UseExceptions()
src_ds = ogr.Open(input_path)
src_layer = src_ds.GetLayer()
driver = ogr.GetDriverByName('ESRI Shapefile')
dst_ds = driver.CreateDataSource(output_path)
dst_layer = dst_ds.CreateLayer(
src_layer.GetName(),
src_layer.GetSpatialRef(),
src_layer.GetGeomType()
)
# 复制字段
src_defn = src_layer.GetLayerDefn()
for i in range(src_defn.GetFieldCount()):
dst_layer.CreateField(src_defn.GetFieldDefn(i))
dst_defn = dst_layer.GetLayerDefn()
for feature in src_layer:
geom = feature.GetGeometryRef()
if geom:
if preserve_topology:
simplified = geom.SimplifyPreserveTopology(tolerance)
else:
simplified = geom.Simplify(tolerance)
out_feat = ogr.Feature(dst_defn)
out_feat.SetGeometry(simplified)
for i in range(src_defn.GetFieldCount()):
out_feat.SetField(i, feature.GetField(i))
dst_layer.CreateFeature(out_feat)
src_ds = None
dst_ds = None
12.2.2 凸包和最小边界
from osgeo import ogr
def create_convex_hull(input_path, output_path):
"""创建凸包"""
src_ds = ogr.Open(input_path)
src_layer = src_ds.GetLayer()
# 合并所有几何
multi = ogr.Geometry(ogr.wkbGeometryCollection)
for feature in src_layer:
geom = feature.GetGeometryRef()
if geom:
multi.AddGeometry(geom)
# 计算凸包
convex_hull = multi.ConvexHull()
# 保存
driver = ogr.GetDriverByName('ESRI Shapefile')
dst_ds = driver.CreateDataSource(output_path)
dst_layer = dst_ds.CreateLayer(
'convex_hull',
src_layer.GetSpatialRef(),
ogr.wkbPolygon
)
out_feat = ogr.Feature(dst_layer.GetLayerDefn())
out_feat.SetGeometry(convex_hull)
dst_layer.CreateFeature(out_feat)
src_ds = None
dst_ds = None
def create_minimum_bounding_rectangle(input_path, output_path):
"""创建最小外接矩形"""
src_ds = ogr.Open(input_path)
src_layer = src_ds.GetLayer()
# 获取范围
extent = src_layer.GetExtent()
# 创建矩形
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(extent[0], extent[2]) # 左下
ring.AddPoint(extent[1], extent[2]) # 右下
ring.AddPoint(extent[1], extent[3]) # 右上
ring.AddPoint(extent[0], extent[3]) # 左上
ring.AddPoint(extent[0], extent[2]) # 闭合
polygon = ogr.Geometry(ogr.wkbPolygon)
polygon.AddGeometry(ring)
# 保存
driver = ogr.GetDriverByName('ESRI Shapefile')
dst_ds = driver.CreateDataSource(output_path)
dst_layer = dst_ds.CreateLayer(
'mbr',
src_layer.GetSpatialRef(),
ogr.wkbPolygon
)
out_feat = ogr.Feature(dst_layer.GetLayerDefn())
out_feat.SetGeometry(polygon)
dst_layer.CreateFeature(out_feat)
src_ds = None
dst_ds = None
12.2.3 质心和点在面内
from osgeo import ogr
def create_centroids(input_path, output_path):
"""创建质心点"""
src_ds = ogr.Open(input_path)
src_layer = src_ds.GetLayer()
driver = ogr.GetDriverByName('ESRI Shapefile')
dst_ds = driver.CreateDataSource(output_path)
dst_layer = dst_ds.CreateLayer(
'centroids',
src_layer.GetSpatialRef(),
ogr.wkbPoint
)
# 复制字段
src_defn = src_layer.GetLayerDefn()
for i in range(src_defn.GetFieldCount()):
dst_layer.CreateField(src_defn.GetFieldDefn(i))
dst_defn = dst_layer.GetLayerDefn()
for feature in src_layer:
geom = feature.GetGeometryRef()
if geom:
centroid = geom.Centroid()
out_feat = ogr.Feature(dst_defn)
out_feat.SetGeometry(centroid)
for i in range(src_defn.GetFieldCount()):
out_feat.SetField(i, feature.GetField(i))
dst_layer.CreateFeature(out_feat)
src_ds = None
dst_ds = None
def create_point_on_surface(input_path, output_path):
"""创建面内点（保证点在多边形内部）"""
src_ds = ogr.Open(input_path)
src_layer = src_ds.GetLayer()
driver = ogr.GetDriverByName('ESRI Shapefile')
dst_ds = driver.CreateDataSource(output_path)
dst_layer = dst_ds.CreateLayer(
'points_on_surface',
src_layer.GetSpatialRef(),
ogr.wkbPoint
)
src_defn = src_layer.GetLayerDefn()
for i in range(src_defn.GetFieldCount()):
dst_layer.CreateField(src_defn.GetFieldDefn(i))
dst_defn = dst_layer.GetLayerDefn()
for feature in src_layer:
geom = feature.GetGeometryRef()
if geom:
point = geom.PointOnSurface()
out_feat = ogr.Feature(dst_defn)
out_feat.SetGeometry(point)
for i in range(src_defn.GetFieldCount()):
out_feat.SetField(i, feature.GetField(i))
dst_layer.CreateFeature(out_feat)
src_ds = None
dst_ds = None
12.3 拓扑处理
12.3.1 几何验证
from osgeo import ogr
def validate_geometries(input_path):
"""验证几何有效性"""
ds = ogr.Open(input_path)
layer = ds.GetLayer()
invalid_features = []
for feature in layer:
geom = feature.GetGeometryRef()
if geom:
if not geom.IsValid():
invalid_features.append({
'fid': feature.GetFID(),
'reason': geom.GetGeometryName(),
})
ds = None
print(f"无效几何数量: {len(invalid_features)}")
for item in invalid_features:
print(f" FID {item['fid']}: {item['reason']}")
return invalid_features
def fix_geometries(input_path, output_path):
"""修复几何"""
src_ds = ogr.Open(input_path)
src_layer = src_ds.GetLayer()
driver = ogr.GetDriverByName('ESRI Shapefile')
dst_ds = driver.CreateDataSource(output_path)
dst_layer = dst_ds.CreateLayer(
src_layer.GetName(),
src_layer.GetSpatialRef(),
src_layer.GetGeomType()
)
src_defn = src_layer.GetLayerDefn()
for i in range(src_defn.GetFieldCount()):
dst_layer.CreateField(src_defn.GetFieldDefn(i))
dst_defn = dst_layer.GetLayerDefn()
fixed_count = 0
for feature in src_layer:
geom = feature.GetGeometryRef()
if geom:
if not geom.IsValid():
# 尝试修复
fixed_geom = geom.MakeValid()
if fixed_geom and fixed_geom.IsValid():
geom = fixed_geom
fixed_count += 1
out_feat = ogr.Feature(dst_defn)
out_feat.SetGeometry(geom)
for i in range(src_defn.GetFieldCount()):
out_feat.SetField(i, feature.GetField(i))
dst_layer.CreateFeature(out_feat)
print(f"修复了 {fixed_count} 个几何")
src_ds = None
dst_ds = None
12.3.2 消除缝隙和重叠
from osgeo import ogr
def dissolve_layer(input_path, output_path, field=None):
"""溶解图层"""
src_ds = ogr.Open(input_path)
src_layer = src_ds.GetLayer()
driver = ogr.GetDriverByName('ESRI Shapefile')
dst_ds = driver.CreateDataSource(output_path)
dst_layer = dst_ds.CreateLayer(
'dissolved',
src_layer.GetSpatialRef(),
ogr.wkbMultiPolygon
)
if field:
# 按字段溶解
src_defn = src_layer.GetLayerDefn()
field_idx = src_defn.GetFieldIndex(field)
field_defn = src_defn.GetFieldDefn(field_idx)
dst_layer.CreateField(field_defn)
# 收集各组的几何
groups = {}
for feature in src_layer:
value = feature.GetField(field)
geom = feature.GetGeometryRef()
if geom:
if value not in groups:
groups[value] = []
groups[value].append(geom.Clone())
# 溶解各组
dst_defn = dst_layer.GetLayerDefn()
for value, geoms in groups.items():
union = geoms[0]
for g in geoms[1:]:
union = union.Union(g)
out_feat = ogr.Feature(dst_defn)
out_feat.SetField(field, value)
out_feat.SetGeometry(union)
dst_layer.CreateFeature(out_feat)
else:
# 溶解所有
union = None
for feature in src_layer:
geom = feature.GetGeometryRef()
if geom:
if union is None:
union = geom.Clone()
else:
union = union.Union(geom)
if union:
out_feat = ogr.Feature(dst_layer.GetLayerDefn())
out_feat.SetGeometry(union)
dst_layer.CreateFeature(out_feat)
src_ds = None
dst_ds = None
12.4 属性处理
12.4.1 属性计算
from osgeo import ogr
def calculate_field(input_path, field_name, expression, field_type=ogr.OFTReal):
"""计算字段值"""
ds = ogr.Open(input_path, 1) # 可写模式
layer = ds.GetLayer()
# 检查字段是否存在，不存在则创建
defn = layer.GetLayerDefn()
field_idx = defn.GetFieldIndex(field_name)
if field_idx < 0:
new_field = ogr.FieldDefn(field_name, field_type)
layer.CreateField(new_field)
field_idx = defn.GetFieldIndex(field_name)
# 计算值
for feature in layer:
geom = feature.GetGeometryRef()
# 内置变量
area = geom.GetArea() if geom else 0
length = geom.Length() if geom else 0
# 计算表达式
try:
value = eval(expression)
except:
value = 0
feature.SetField(field_name, value)
layer.SetFeature(feature)
ds = None
def calculate_geometry_fields(input_path):
"""计算几何属性字段"""
ds = ogr.Open(input_path, 1)
layer = ds.GetLayer()
# 添加字段
geom_type = layer.GetGeomType()
if geom_type in [ogr.wkbPolygon, ogr.wkbMultiPolygon]:
layer.CreateField(ogr.FieldDefn('area', ogr.OFTReal))
layer.CreateField(ogr.FieldDefn('perimeter', ogr.OFTReal))
if geom_type in [ogr.wkbLineString, ogr.wkbMultiLineString]:
layer.CreateField(ogr.FieldDefn('length', ogr.OFTReal))
layer.CreateField(ogr.FieldDefn('centroid_x', ogr.OFTReal))
layer.CreateField(ogr.FieldDefn('centroid_y', ogr.OFTReal))
# 计算值
for feature in layer:
geom = feature.GetGeometryRef()
if geom:
if geom_type in [ogr.wkbPolygon, ogr.wkbMultiPolygon]:
feature.SetField('area', geom.GetArea())
feature.SetField('perimeter', geom.Boundary().Length())
if geom_type in [ogr.wkbLineString, ogr.wkbMultiLineString]:
feature.SetField('length', geom.Length())
centroid = geom.Centroid()
feature.SetField('centroid_x', centroid.GetX())
feature.SetField('centroid_y', centroid.GetY())
layer.SetFeature(feature)
ds = None
12.4.2 属性连接
from osgeo import ogr
def attribute_join(target_path, join_table_path, target_key, join_key, output_path):
"""属性连接"""
target_ds = ogr.Open(target_path)
target_layer = target_ds.GetLayer()
join_ds = ogr.Open(join_table_path)
join_layer = join_ds.GetLayer()
# 创建输出
driver = ogr.GetDriverByName('ESRI Shapefile')
out_ds = driver.CreateDataSource(output_path)
out_layer = out_ds.CreateLayer(
target_layer.GetName(),
target_layer.GetSpatialRef(),
target_layer.GetGeomType()
)
# 复制目标字段
target_defn = target_layer.GetLayerDefn()
for i in range(target_defn.GetFieldCount()):
out_layer.CreateField(target_defn.GetFieldDefn(i))
# 添加连接字段
join_defn = join_layer.GetLayerDefn()
join_key_idx = join_defn.GetFieldIndex(join_key)
for i in range(join_defn.GetFieldCount()):
if i != join_key_idx:
field = join_defn.GetFieldDefn(i)
new_field = ogr.FieldDefn(f'join_{field.GetName()}', field.GetType())
out_layer.CreateField(new_field)
# 建立连接表索引
join_dict = {}
for feature in join_layer:
key_value = feature.GetField(join_key)
join_dict[key_value] = feature
# 执行连接
out_defn = out_layer.GetLayerDefn()
for target_feat in target_layer:
out_feat = ogr.Feature(out_defn)
# 复制几何
geom = target_feat.GetGeometryRef()
if geom:
out_feat.SetGeometry(geom.Clone())
# 复制目标属性
for i in range(target_defn.GetFieldCount()):
out_feat.SetField(i, target_feat.GetField(i))
# 连接属性
key_value = target_feat.GetField(target_key)
if key_value in join_dict:
join_feat = join_dict[key_value]
j = 0
for i in range(join_defn.GetFieldCount()):
if i != join_key_idx:
out_feat.SetField(
target_defn.GetFieldCount() + j,
join_feat.GetField(i)
)
j += 1
out_layer.CreateFeature(out_feat)
target_ds = None
join_ds = None
out_ds = None
12.5 本章小结
本章介绍了矢量数据高级处理：
- 空间分析：缓冲区、叠加分析、空间连接
- 几何操作：简化、凸包、质心
- 拓扑处理：验证、修复、溶解
- 属性处理：计算、连接
12.6 思考与练习
- 实现一个支持多种空间关系的空间查询函数。
- 编写代码检测并修复多边形的自相交问题。
- 实现泰森多边形（Voronoi）生成。
- 如何高效处理大量要素的空间连接？
- 实现支持多字段的溶解函数。
浙公网安备 33010602011771号