几何叠加分析在工程实践中如何应用？

第17章：几何叠加分析（Overlay） 叠加分析是 GIS 空间分析的核心操作之一。与空间连接仅关联属性不同，叠加分析会同时对几何对象进行集合运算（交集、并集、差集等），并合并两个数据集的属性信息。本章将详细介绍 GeoPandas 中 overlay() 函数的使用方法及各种叠加类型的原理与应用。 17.1 叠加分析概述 17.1.1 什么是叠加分析 叠加分析（Overlay Analysis）是将两个面图层进行几何集合运算，生成新的几何要素并合并双方属性的过程。它是传统 GIS 中最经典的空间分析方法之一。 与空间连接的区别： 操作 空间连接（sjoin） 叠加分析（overlay） 几何变化 不修改几何 生成新的几何 属性处理 附加匹配要素的属性 合并两层的属性 输入限制 任意几何类型 主要用于面数据 输出几何 保留左表几何 新生成的几何 典型用途 属性关联 空间分区分析 17.1.2 叠加类型概览 GeoPandas 的 overlay() 支持以下叠加类型： 类型 关键字 说明 交集 'intersection' 两层重叠的部分 并集 'union' 两层的全部区域 差集 'difference' 左层减去右层的部分 对称差集 'symmetric_difference' 两层不重叠的部分 身份叠加 'identity' 保留左层全部，叠加右层 import geopandas as gpd from shapely.geometry import Polygon # 创建两个重叠的多边形层 layer_a = gpd.GeoDataFrame({ 'name_a': ['A1', 'A2'], 'value_a': [100, 200], 'geometry': [ Polygon([(0, 0), (3, 0), (3, 3), (0, 3)]), Polygon([(3, 0), (6, 0), (6, 3), (3, 3)]) ] }) layer_b = gpd.GeoDataFrame({ 'name_b': ['B1', 'B2'], 'value_b': [10, 20], 'geometry': [ Polygon([(1, 1), (4, 1), (4, 4), (1, 4)]), Polygon([(4, 1), (7, 1), (7, 4), (4, 4)]) ] }) print("图层 A:") print(layer_a) print("\n图层 B:") print(layer_b) 17.1.3 几何集合运算示意 图层 A 图层 B 交集(intersection) +---+---+ +---+---+ +--+--+ | A1| A2| | B1| B2| |AB|AB| | | | | | | |11|21| +---+---+ +---+---+ +--+--+ 并集(union) 差集(difference) 身份(identity) +--+--+--+--+ +--+ +--+--+--+ |A | A| A| | |A1| |A1|AB|A1| |1 |AB|AB|B2| | | | |11| | +--+--+--+--+ +--+ +--+--+--+ | |B1|B2| | +--+--+--+--+ 17.2 overlay() 基础用法 17.2.1 基本语法 geopandas.overlay(df1, df2, how='intersection', keep_geom_type=False, make_valid=True) 也可以使用 GeoDataFrame 的方法形式： df1.overlay(df2, how='intersection') 17.2.2 参数说明 参数 类型 说明 df1 GeoDataFrame 第一个输入图层（左层） df2 GeoDataFrame 第二个输入图层（右层） how str 叠加类型 keep_geom_type bool 是否只保留与输入相同的几何类型 make_valid bool 是否自动修复无效几何 17.2.3 第一个叠加分析 # 执行交集分析 intersection = gpd.overlay(layer_a, layer_b, how='intersection') print("交集分析结果:") print(intersection) print(f"\n结果要素数: {len(intersection)}") print(f"包含列: {intersection.columns.tolist()}") 输出结果将包含： 两层重叠部分的 新几何 来自两层的 所有属性列 17.3 交集分析（intersection） 17.3.1 基本原理 交集分析提取两个图层共同覆盖的区域，结果只包含两层重叠的部分。每个结果要素拥有来自两层的属性。 intersection = gpd.overlay(layer_a, layer_b, how='intersection') print("交集结果:") for idx, row in intersection.iterrows(): print(f" {row['name_a']}-{row['name_b']}: " f"面积={row.geometry.area:.2f}, " f"value_a={row['value_a']}, value_b={row['value_b']}") 17.3.2 实际应用 - 土地利用与洪水风险交叉分析 import geopandas as gpd from shapely.geometry import Polygon # 土地利用图层 landuse = gpd.GeoDataFrame({ 'landuse': ['居住用地', '商业用地', '农业用地'], 'population': [50000, 20000, 5000], 'geometry': [ Polygon([(0, 0), (4, 0), (4, 3), (0, 3)]), Polygon([(4, 0), (8, 0), (8, 3), (4, 3)]), Polygon([(0, 3), (8, 3), (8, 6), (0, 6)]) ] }) # 洪水风险区 flood_risk = gpd.GeoDataFrame({ 'risk_level': ['高风险', '中风险'], 'geometry': [ Polygon([(0, 0), (5, 0), (5, 4), (0, 4)]), Polygon([(5, 0), (8, 0), (8, 4), (5, 4)]) ] }) # 交叉分析 risk_landuse = gpd.overlay(landuse, flood_risk, how='intersection') print("土地利用与洪水风险交叉分析:") print(risk_landuse[['landuse', 'risk_level', 'population']]) # 计算各风险等级下受影响人口的面积加权估算 risk_landuse['overlap_area'] = risk_landuse.geometry.area for idx, row in risk_landuse.iterrows(): original = landuse[landuse['landuse'] == row['landuse']].geometry.iloc[0] ratio = row['overlap_area'] / original.area risk_landuse.loc[idx, 'estimated_pop'] = row['population'] * ratio print("\n各风险等级受影响人口估算:") pop_by_risk = risk_landuse.groupby('risk_level')['estimated_pop'].sum() print(pop_by_risk) 17.4 并集分析（union） 17.4.1 基本原理 并集分析将两层的全部区域合并，重叠区域拆分为独立要素并保留双方属性，非重叠区域保留原始属性（另一方属性为 NaN）。 union = gpd.overlay(layer_a, layer_b, how='union') print("并集结果:") print(union[['name_a', 'name_b', 'value_a', 'value_b']]) print(f"\n结果要素数: {len(union)}") 17.4.2 并集结果的组成 区域类型 name_a name_b 说明 仅 A 覆盖 有值 NaN 不在 B 范围内的 A 区域 A 与 B 重叠 有值 有值 两层共同覆盖的区域 仅 B 覆盖 NaN 有值 不在 A 范围内的 B 区域 17.4.3 实际应用 - 合并两个行政区划 import geopandas as gpd from shapely.geometry import Polygon # 两个不同来源的区域划分 zoning_a = gpd.GeoDataFrame({ 'zone_a': ['住宅区', '工业区'], 'a_code': ['R1', 'I1'], 'geometry': [ Polygon([(0, 0), (5, 0), (5, 5), (0, 5)]), Polygon([(5, 0), (10, 0), (10, 5), (5, 5)]) ] }) zoning_b = gpd.GeoDataFrame({ 'zone_b': ['禁建区', '可建区'], 'b_code': ['N', 'Y'], 'geometry': [ Polygon([(0, 0), (3, 0), (3, 5), (0, 5)]), Polygon([(3, 0), (10, 0), (10, 5), (3, 5)]) ] }) # 并集叠加 union_result = gpd.overlay(zoning_a, zoning_b, how='union') print("并集叠加结果:") print(union_result[['zone_a', 'zone_b', 'a_code', 'b_code']]) # 创建综合管理分区 union_result['管理分区'] = ( union_result['zone_a'].fillna('') + '-' + union_result['zone_b'].fillna('') ) print("\n综合管理分区:") print(union_result[['管理分区', 'geometry']].to_string()) 17.5 差集分析（difference） 17.5.1 基本原理 差集分析从左层中减去与右层重叠的部分，只保留左层中不被右层覆盖的区域。结果只包含左层的属性。 difference = gpd.overlay(layer_a, layer_b, how='difference') print("差集结果（A - B）:") print(difference) print(f"\n结果要素数: {len(difference)}") 17.5.2 差集与擦除 差集分析在 GIS 中也称为"擦除"（Erase）操作。它常用于从数据中排除特定区域： import geopandas as gpd from shapely.geometry import Polygon # 研究区域 study_area = gpd.GeoDataFrame({ 'name': ['研究区'], 'geometry': [Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])] }) # 需要排除的区域（如水体、保护区） excluded = gpd.GeoDataFrame({ 'type': ['湖泊', '自然保护区'], 'geometry': [ Polygon([(2, 2), (4, 2), (4, 4), (2, 4)]), Polygon([(7, 6), (9, 6), (9, 9), (7, 9)]) ] }) # 从研究区中排除特定区域 valid_area = gpd.overlay(study_area, excluded, how='difference') print(f"原始面积: {study_area.geometry.area.sum():.1f}") print(f"排除后面积: {valid_area.geometry.area.sum():.1f}") print(f"排除面积: {study_area.geometry.area.sum() - valid_area.geometry.area.sum():.1f}") 17.5.3 注意事项 差集操作是非对称的：A - B ≠ B - A 结果只包含左层的属性，不包含右层属性 如果右层完全覆盖左层某个要素，该要素将从结果中消失 # 差集是非对称的 diff_ab = gpd.overlay(layer_a, layer_b, how='difference') # A - B diff_ba = gpd.overlay(layer_b, layer_a, how='difference') # B - A print(f"A - B 面积: {diff_ab.geometry.area.sum():.2f}") print(f"B - A 面积: {diff_ba.geometry.area.sum():.2f}") 17.6 对称差集（symmetric_difference） 17.6.1 基本原理 对称差集返回两层中不重叠的区域，即 (A - B) + (B - A)。重叠部分被排除。 sym_diff = gpd.overlay(layer_a, layer_b, how='symmetric_difference') print("对称差集结果:") print(sym_diff[['name_a', 'name_b']]) print(f"\n结果要素数: {len(sym_diff)}") 17.6.2 与其他叠加类型的关系 对称差集 = 并集 - 交集 验证这个关系： union_result = gpd.overlay(layer_a, layer_b, how='union') intersection_result = gpd.overlay(layer_a, layer_b, how='intersection') sym_diff_result = gpd.overlay(layer_a, layer_b, how='symmetric_difference') union_area = union_result.geometry.area.sum() intersection_area = intersection_result.geometry.area.sum() sym_diff_area = sym_diff_result.geometry.area.sum() print(f"并集面积: {union_area:.2f}") print(f"交集面积: {intersection_area:.2f}") print(f"对称差集面积: {sym_diff_area:.2f}") print(f"并集 - 交集 = {union_area - intersection_area:.2f}") print(f"等式验证: {abs(sym_diff_area - (union_area - intersection_area)) < 1e-10}") 17.6.3 应用场景 对称差集常用于： 变化检测：比较两个时期的用地边界，找出变化区域 差异分析：找出两个规划方案之间不同的区域 # 两个时期的城市边界 urban_2020 = gpd.GeoDataFrame({ 'year': [2020], 'geometry': [Polygon([(0, 0), (6, 0), (6, 6), (0, 6)])] }) urban_2023 = gpd.GeoDataFrame({ 'year': [2023], 'geometry': [Polygon([(1, 1), (8, 1), (8, 7), (1, 7)])] }) # 对称差集 - 找出变化区域 changed = gpd.overlay(urban_2020, urban_2023, how='symmetric_difference') print(f"变化区域面积: {changed.geometry.area.sum():.2f}") print(f"其中 2020 独有: {changed[changed['year'].notna()].geometry.area.sum():.2f}") 17.7 身份叠加（identity） 17.7.1 基本原理 身份叠加保留左层的全部区域，但将与右层重叠的部分进行拆分，并附加右层的属性。非重叠区域保留原始几何和属性。 identity = gpd.overlay(layer_a, layer_b, how='identity') print("身份叠加结果:") print(identity[['name_a', 'name_b', 'value_a', 'value_b']]) print(f"\n结果要素数: {len(identity)}") 17.7.2 身份叠加与交集的区别 特征 交集（intersection） 身份叠加（identity） 保留区域 仅重叠区域 左层全部区域 非重叠左层 丢弃 保留（右层属性为 NaN） 非重叠右层 丢弃 丢弃 右层属性 始终有值 重叠区域有值，其他为 NaN 17.7.3 实际应用 身份叠加适用于需要保持左层完整性的场景： import geopandas as gpd from shapely.geometry import Polygon # 地块数据（需要保持完整） parcels = gpd.GeoDataFrame({ 'parcel_id': ['P001', 'P002', 'P003'], 'owner': ['张三', '李四', '王五'], 'geometry': [ Polygon([(0, 0), (3, 0), (3, 3), (0, 3)]), Polygon([(3, 0), (6, 0), (6, 3), (3, 3)]), Polygon([(6, 0), (9, 0), (9, 3), (6, 3)]) ] }) # 规划限制区域（部分覆盖） restrictions = gpd.GeoDataFrame({ 'restriction': ['限高区', '绿化带'], 'max_height': [15, 5], 'geometry': [ Polygon([(1, 0), (4, 0), (4, 3), (1, 3)]), Polygon([(7, 0), (10, 0), (10, 3), (7, 3)]) ] }) # 身份叠加 - 保留所有地块，添加规划限制信息 result = gpd.overlay(parcels, restrictions, how='identity') print("地块与规划限制叠加:") print(result[['parcel_id', 'owner', 'restriction', 'max_height']]) 17.8 叠加分析的属性处理 17.8.1 属性列的合并规则 叠加分析自动合并两层的所有非几何列： import geopandas as gpd from shapely.geometry import Polygon gdf1 = gpd.GeoDataFrame({ 'id': [1], 'name': ['区域A'], 'value': [100], 'geometry': [Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])] }) gdf2 = gpd.GeoDataFrame({ 'id': [1], 'type': ['类型B'], 'score': [0.8], 'geometry': [Polygon([(1, 1), (3, 1), (3, 3), (1, 3)])] }) result = gpd.overlay(gdf1, gdf2, how='union') print("合并后的列:") print(result.columns.tolist()) # 同名列自动添加后缀: id_1, id_2 17.8.2 同名列处理 当两层存在同名的非几何列时，GeoPandas 会自动添加数字后缀来区分： # gdf1 和 gdf2 都有 'id' 列 result = gpd.overlay(gdf1, gdf2, how='intersection') print(result.columns.tolist()) # 输出: ['id_1', 'name', 'value', 'id_2', 'type', 'score', 'geometry'] 17.8.3 NaN 值的含义 在并集、身份叠加等操作中，非重叠区域的对方属性为 NaN： union = gpd.overlay(gdf1, gdf2, how='union') # 检查 NaN 值 print("含 NaN 的行:") print(union[union.isna().any(axis=1)]) # 填充 NaN 值 union['value'] = union['value'].fillna(0) union['score'] = union['score'].fillna(0) print("\n填充后:") print(union) 17.8.4 面积加权属性计算 叠加分析后，常需要对数值属性进行面积加权重新分配： import geopandas as gpd from shapely.geometry import Polygon # 人口普查区 census = gpd.GeoDataFrame({ 'census_id': ['C1', 'C2'], 'population': [10000, 8000], 'geometry': [ Polygon([(0, 0), (5, 0), (5, 5), (0, 5)]), Polygon([(5, 0), (10, 0), (10, 5), (5, 5)]) ] }) # 学区 school_districts = gpd.GeoDataFrame({ 'school': ['学校A学区', '学校B学区'], 'geometry': [ Polygon([(0, 0), (7, 0), (7, 5), (0, 5)]), Polygon([(7, 0), (10, 0), (10, 5), (7, 5)]) ] }) # 交集分析 overlay_result = gpd.overlay(census, school_districts, how='intersection') # 计算面积加权人口 for idx, row in overlay_result.iterrows(): original = census[census['census_id'] == row['census_id']].geometry.iloc[0] area_ratio = row.geometry.area / original.area overlay_result.loc[idx, 'weighted_pop'] = row['population'] * area_ratio print("各学区人口估算:") school_pop = overlay_result.groupby('school')['weighted_pop'].sum() print(school_pop) 17.9 性能优化与注意事项 17.9.1 keep_geom_type 参数 叠加运算可能产生混合几何类型。keep_geom_type=True 只保留与输入相同类型的几何： # 默认行为 - 可能产生点、线、面混合结果 result_all = gpd.overlay(layer_a, layer_b, how='intersection', keep_geom_type=False) # 仅保留面类型 result_poly = gpd.overlay(layer_a, layer_b, how='intersection', keep_geom_type=True) print(f"保留所有类型: {len(result_all)} 个要素") print(f"仅保留面类型: {len(result_poly)} 个要素") 17.9.2 make_valid 参数 无效几何会导致叠加分析失败或产生错误结果。make_valid=True（默认）自动修复无效几何： # 自动修复无效几何 result = gpd.overlay(gdf1, gdf2, how='intersection', make_valid=True) # 手动检查和修复 print(f"gdf1 有效几何数: {gdf1.is_valid.sum()}/{len(gdf1)}") print(f"gdf2 有效几何数: {gdf2.is_valid.sum()}/{len(gdf2)}") # 手动修复 gdf1['geometry'] = gdf1.geometry.make_valid() gdf2['geometry'] = gdf2.geometry.make_valid() 17.9.3 性能优化策略 import geopandas as gpd import time # 1. 先裁剪到共同范围 bounds_a = layer_a.total_bounds bounds_b = layer_b.total_bounds common_bounds = [ max(bounds_a[0], bounds_b[0]), max(bounds_a[1], bounds_b[1]), min(bounds_a[2], bounds_b[2]), min(bounds_a[3], bounds_b[3]) ] # 2. 简化几何（在容许误差范围内） layer_a_simple = layer_a.copy() layer_a_simple['geometry'] = layer_a_simple.geometry.simplify(tolerance=0.01) # 3. 移除不需要的列（减少内存占用） cols_needed = ['name_a', 'value_a', 'geometry'] layer_a_clean = layer_a[cols_needed] # 4. 执行叠加分析 start = time.time() result = gpd.overlay(layer_a_clean, layer_b, how='intersection') print(f"叠加分析耗时: {time.time() - start:.3f} 秒") 17.9.4 常见问题与解决方案 问题 原因 解决方案 TopologyException 无效几何 使用 make_valid=True 或预处理 结果包含碎片 精度问题 设置 keep_geom_type=True，过滤微小面积 内存溢出 数据量太大 分块处理或简化几何 结果属性为 NaN 非重叠区域 根据叠加类型预期行为处理 CRS 不一致 坐标系不同 统一 CRS 后再叠加 # 过滤微小碎片 result = gpd.overlay(layer_a, layer_b, how='intersection') min_area = 0.001 # 最小面积阈值 result = result[result.geometry.area > min_area] print(f"过滤碎片后要素数: {len(result)}") 17.10 实际应用案例 17.10.1 案例一：城市规划用地冲突分析 import geopandas as gpd from shapely.geometry import Polygon # 城市规划用地 planning = gpd.GeoDataFrame({ 'plan_type': ['商业区', '住宅区', '工业区', '绿化区'], 'plan_code': ['C', 'R', 'I', 'G'], 'geometry': [ Polygon([(0, 0), (4, 0), (4, 4), (0, 4)]), Polygon([(4, 0), (8, 0), (8, 4), (4, 4)]), Polygon([(0, 4), (4, 4), (4, 8), (0, 8)]), Polygon([(4, 4), (8, 4), (8, 8), (4, 8)]) ] }) # 生态保护红线 redline = gpd.GeoDataFrame({ 'protect_level': ['一级保护', '二级保护'], 'geometry': [ Polygon([(3, 3), (6, 3), (6, 6), (3, 6)]), Polygon([(6, 5), (9, 5), (9, 8), (6, 8)]) ] }) # 交集分析 - 找出规划用地与保护红线的冲突区域 conflict = gpd.overlay(planning, redline, how='intersection') conflict['conflict_area'] = conflict.geometry.area print("用地规划与生态红线冲突分析:") print(conflict[['plan_type', 'protect_level', 'conflict_area']]) print(f"\n总冲突面积: {conflict['conflict_area'].sum():.2f}") # 冲突汇总 conflict_summary = conflict.groupby('plan_type')['conflict_area'].sum() print("\n各类用地冲突面积:") print(conflict_summary) 17.10.2 案例二：多规合一分析 import geopandas as gpd from shapely.geometry import Polygon # 城市总体规划 urban_plan = gpd.GeoDataFrame({ 'urban_use': ['建设用地', '建设用地', '非建设用地'], 'geometry': [ Polygon([(0, 0), (5, 0), (5, 3), (0, 3)]), Polygon([(5, 0), (10, 0), (10, 3), (5, 3)]), Polygon([(0, 3), (10, 3), (10, 6), (0, 6)]) ] }) # 土地利用规划 land_plan = gpd.GeoDataFrame({ 'land_use': ['耕地', '建设用地', '林地'], 'geometry': [ Polygon([(0, 0), (3, 0), (3, 6), (0, 6)]), Polygon([(3, 0), (7, 0), (7, 6), (3, 6)]), Polygon([(7, 0), (10, 0), (10, 6), (7, 6)]) ] }) # 并集叠加 - 生成综合分区 combined = gpd.overlay(urban_plan, land_plan, how='union') combined['combo'] = combined['urban_use'].fillna('未覆盖') + ' | ' + combined['land_use'].fillna('未覆盖') combined['area'] = combined.geometry.area print("多规合一分析结果:") for _, row in combined.iterrows(): print(f" {row['combo']}: 面积 = {row['area']:.2f}") # 识别矛盾区域 combined['conflict'] = ( (combined['urban_use'] == '建设用地') & (combined['land_use'].isin(['耕地', '林地'])) ) conflicts = combined[combined['conflict']] print(f"\n矛盾区域数: {len(conflicts)}") print(f"矛盾总面积: {conflicts['area'].sum():.2f}") 17.10.3 案例三：缓冲区叠加分析 import geopandas as gpd from shapely.geometry import Point, Polygon # 污染源 sources = gpd.GeoDataFrame({ 'source': ['工厂A', '工厂B'], 'pollution_level': ['高', '中'], 'geometry': [Point(3, 3), Point(7, 5)] }) # 生成缓冲区（假设已投影到平面坐标系） buffer_1km = sources.copy() buffer_1km['geometry'] = sources.geometry.buffer(2) buffer_1km['buffer_dist'] = '2km' # 敏感区域 sensitive = gpd.GeoDataFrame({ 'area_type': ['学校', '医院', '居民区'], 'geometry': [ Polygon([(1, 4), (3, 4), (3, 6), (1, 6)]), Polygon([(5, 2), (7, 2), (7, 4), (5, 4)]), Polygon([(6, 5), (9, 5), (9, 8), (6, 8)]) ] }) # 叠加分析 - 找出受影响的敏感区域 affected = gpd.overlay(sensitive, buffer_1km, how='intersection') affected['affected_area'] = affected.geometry.area print("受污染影响的敏感区域:") print(affected[['area_type', 'source', 'pollution_level', 'affected_area']]) 17.11 本章小结 本章全面介绍了 GeoPandas 中的几何叠加分析。主要内容回顾： 叠加类型对比 类型 关键字 保留区域 属性来源 交集 intersection 仅重叠 两层 并集 union 全部 两层（含 NaN） 差集 difference 左层非重叠 仅左层 对称差集 symmetric_difference 两层非重叠 两层（含 NaN） 身份叠加 identity 左层全部 两层（含 NaN） 关键参数 参数 说明 默认值 how 叠加类型 'intersection' keep_geom_type 保留相同几何类型 False make_valid 自动修复无效几何 True 使用建议 选择正确的叠加类型：根据分析目标选择合适的 how 参数 预处理数据：确保 CRS 一致、几何有效 处理碎片：使用 keep_geom_type=True 和面积过滤处理碎片 属性重分配：叠加后常需对数值属性进行面积加权计算 性能优化：对大数据集先裁剪到感兴趣区域再叠加 在下一章中，我们将学习裁剪与掩膜操作，这是叠加分析的一种特殊应用，用于将数据限制在特定区域范围内。