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)

geopandas.overlay(df1, df2, how='intersection', keep_geom_type=False, make_valid=True)

df1.overlay(df2, how='intersection')

# 执行交集分析
intersection = gpd.overlay(layer_a, layer_b, how='intersection')
print("交集分析结果:")
print(intersection)
print(f"\n结果要素数: {len(intersection)}")
print(f"包含列: {intersection.columns.tolist()}")

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']}")

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)

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)}")

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())

difference = gpd.overlay(layer_a, layer_b, how='difference')

print("差集结果（A - B）:")
print(difference)
print(f"\n结果要素数: {len(difference)}")

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}")

# 差集是非对称的
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}")

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)}")

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}")

# 两个时期的城市边界
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}")

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)}")

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']])

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

# gdf1 和 gdf2 都有 'id' 列
result = gpd.overlay(gdf1, gdf2, how='intersection')
print(result.columns.tolist())
# 输出: ['id_1', 'name', 'value', 'id_2', 'type', 'score', 'geometry']

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)

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)

# 默认行为 - 可能产生点、线、面混合结果
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)} 个要素")

# 自动修复无效几何
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()

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} 秒")

# 过滤微小碎片
result = gpd.overlay(layer_a, layer_b, how='intersection')
min_area = 0.001  # 最小面积阈值
result = result[result.geometry.area > min_area]
print(f"过滤碎片后要素数: {len(result)}")

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)

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}")

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']])