import geopandas as gpd
from shapely.geometry import Polygon, Point

# 创建两个 GeoSeries
gs1 = gpd.GeoSeries([
    Polygon([(0, 0), (2, 0), (2, 2), (0, 2)]),
    Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
])

gs2 = gpd.GeoSeries([
    Polygon([(1, 1), (3, 1), (3, 3), (1, 3)]),
    Polygon([(0.5, 0.5), (1.5, 0.5), (1.5, 1.5), (0.5, 1.5)]),
])

# 逐元素合并
result = gs1.union(gs2)
print(result)
print(f"第一个合并结果面积: {result.iloc[0].area}")
print(f"第一个合并结果类型: {result.iloc[0].geom_type}")

# 将所有几何对象与一个公共区域合并
common_area = Polygon([(1, 0), (3, 0), (3, 2), (1, 2)])

result = gs1.union(common_area)
print(result)

from shapely.geometry import Polygon

# 重叠的多边形 → Polygon
poly_a = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
poly_b = Polygon([(1, 1), (3, 1), (3, 3), (1, 3)])
gs = gpd.GeoSeries([poly_a])
result = gs.union(poly_b)
print(f"重叠多边形合并: {result.iloc[0].geom_type}")  # Polygon

# 不重叠的多边形 → MultiPolygon
poly_c = Polygon([(5, 5), (7, 5), (7, 7), (5, 7)])
result2 = gs.union(poly_c)
print(f"分离多边形合并: {result2.iloc[0].geom_type}")  # MultiPolygon

import geopandas as gpd
from shapely.geometry import Polygon, LineString

# 两个重叠的多边形
poly_a = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
poly_b = Polygon([(1, 1), (3, 1), (3, 3), (1, 3)])

gs = gpd.GeoSeries([poly_a])
result = gs.intersection(poly_b)

print(f"交集类型: {result.iloc[0].geom_type}")  # Polygon
print(f"交集面积: {result.iloc[0].area}")        # 1.0
# 交集是 (1,1)-(2,1)-(2,2)-(1,2) 的正方形

# 多边形与线的交集
poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
line = LineString([(-1, 1), (3, 1)])  # 穿过多边形的线

gs = gpd.GeoSeries([poly])
result = gs.intersection(line)
print(f"多边形与线交集: {result.iloc[0]}")
# LineString [(0, 1), (2, 1)] - 线在多边形内的部分

# 不相交的几何对象 → 空几何
poly2 = Polygon([(5, 5), (6, 5), (6, 6), (5, 6)])
result2 = gs.intersection(poly2)
print(f"空交集: {result2.iloc[0].is_empty}")  # True

import geopandas as gpd
from shapely.geometry import box

# 读取全国数据
land_use = gpd.read_file("land_use.shp")

# 定义研究区域
study_area = box(113.5, 22.0, 114.5, 23.0)  # 珠三角区域

# 裁剪到研究区域
clipped = land_use.copy()
clipped['geometry'] = land_use.intersection(study_area)

# 移除空几何
clipped = clipped[~clipped.is_empty]
print(f"裁剪后要素数: {len(clipped)}")

import geopandas as gpd
from shapely.geometry import Polygon

# 两个重叠的多边形
poly_a = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
poly_b = Polygon([(1, 1), (3, 1), (3, 3), (1, 3)])

gs = gpd.GeoSeries([poly_a])

# A - B
diff_ab = gs.difference(poly_b)
print(f"A - B 面积: {diff_ab.iloc[0].area}")  # 3.0

# B - A
gs2 = gpd.GeoSeries([poly_b])
diff_ba = gs2.difference(poly_a)
print(f"B - A 面积: {diff_ba.iloc[0].area}")  # 3.0

# 验证：A的面积 = (A-B)面积 + (A∩B)面积
intersection_area = gs.intersection(poly_b).iloc[0].area
print(f"A面积 = {poly_a.area}, (A-B) + (A∩B) = {diff_ab.iloc[0].area + intersection_area}")

import geopandas as gpd

# 从城市区域中排除公园和水域
city_area = gpd.read_file("city_boundary.shp").geometry.iloc[0]
parks = gpd.read_file("parks.shp")
water = gpd.read_file("water_bodies.shp")

# 合并所有需要排除的区域
excluded = parks.union_all().union(water.union_all())

# 从城市区域中减去排除区域
buildable_area = gpd.GeoSeries([city_area]).difference(excluded)
print(f"可建设面积: {buildable_area.iloc[0].area:.2f}")

# 使用 overlay 进行擦除
result = gpd.overlay(gdf1, gdf2, how='difference')

import geopandas as gpd
from shapely.geometry import Polygon

poly_a = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
poly_b = Polygon([(1, 1), (3, 1), (3, 3), (1, 3)])

gs = gpd.GeoSeries([poly_a])
result = gs.symmetric_difference(poly_b)

print(f"对称差面积: {result.iloc[0].area}")  # 6.0
# A面积(4) + B面积(4) - 2×交集面积(1) = 6

print(f"对称差类型: {result.iloc[0].geom_type}")  # MultiPolygon（两块不相连的区域）

# 对比两个时期的城市建成区，找出变化区域
urban_2010 = gpd.read_file("urban_2010.shp").union_all()
urban_2020 = gpd.read_file("urban_2020.shp").union_all()

# 对称差即为变化区域（新增 + 减少）
changed_area = gpd.GeoSeries([urban_2010]).symmetric_difference(urban_2020)

# 分别计算新增和减少
newly_built = gpd.GeoSeries([urban_2020]).difference(urban_2010)  # 新增区域
demolished = gpd.GeoSeries([urban_2010]).difference(urban_2020)   # 减少区域

print(f"变化总面积: {changed_area.iloc[0].area:.2f}")
print(f"新增面积: {newly_built.iloc[0].area:.2f}")
print(f"减少面积: {demolished.iloc[0].area:.2f}")

import geopandas as gpd
from shapely.geometry import LineString, Point

# 两条几乎相连但有微小间隙的线
line_a = LineString([(0, 0), (1, 0)])
line_b = LineString([(1.01, 0), (2, 0)])  # 起点距离 line_a 终点 0.01

gs = gpd.GeoSeries([line_a])

# 将 line_a 捕捉到 line_b，容差 0.05
snapped = gs.snap(line_b, tolerance=0.05)
print(f"捕捉前终点: {list(line_a.coords)[-1]}")
print(f"捕捉后终点: {list(snapped.iloc[0].coords)[-1]}")
# 捕捉后 line_a 的终点会被移到 (1.01, 0)

# 修复拓扑不一致：确保相邻多边形共享精确的边界
import geopandas as gpd

polygons = gpd.read_file("parcels.shp")

# 对于每对相邻多边形，将其捕捉到彼此
tolerance = 0.001  # 容差（与坐标系单位一致）

for i in range(len(polygons)):
    for j in range(i + 1, len(polygons)):
        if polygons.geometry.iloc[i].distance(polygons.geometry.iloc[j]) < tolerance:
            # 将几何 i 捕捉到几何 j
            polygons.geometry.iloc[i] = polygons.geometry.iloc[i].snap(
                polygons.geometry.iloc[j], tolerance
            )

import geopandas as gpd
from shapely.geometry import Polygon

# 创建多个多边形
gs = gpd.GeoSeries([
    Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
    Polygon([(0.5, 0.5), (1.5, 0.5), (1.5, 1.5), (0.5, 1.5)]),
    Polygon([(1, 1), (2, 1), (2, 2), (1, 2)]),
])

# 合并所有几何对象
merged = gs.unary_union
print(f"合并结果类型: {merged.geom_type}")  # Polygon
print(f"合并结果面积: {merged.area}")

import geopandas as gpd

# 读取县级行政区数据
counties = gpd.read_file("counties.shp")

# 按省份名称分组合并，得到省级边界
provinces = counties.dissolve(by='province_name')
print(f"县级区域数: {len(counties)}")
print(f"省级区域数: {len(provinces)}")

import geopandas as gpd
from shapely.geometry import Polygon

gs = gpd.GeoSeries([
    Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
    Polygon([(0.5, 0), (1.5, 0), (1.5, 1), (0.5, 1)]),
    Polygon([(1, 0), (2, 0), (2, 1), (1, 1)]),
])

# 使用 union_all
merged = gs.union_all()
print(f"合并面积: {merged.area}")  # 2.0
print(f"合并类型: {merged.geom_type}")  # Polygon

import geopandas as gpd
from shapely.geometry import Polygon

gs = gpd.GeoSeries([
    Polygon([(0, 0), (2, 0), (2, 2), (0, 2)]),
    Polygon([(0.5, 0.5), (2.5, 0.5), (2.5, 2.5), (0.5, 2.5)]),
    Polygon([(1, 1), (3, 1), (3, 3), (1, 3)]),
])

# 计算所有几何对象的交集
common = gs.intersection_all()
print(f"公共区域面积: {common.area}")      # 1.0
print(f"公共区域类型: {common.geom_type}")  # Polygon
# 公共区域是 (1,1)-(2,1)-(2,2)-(1,2)

import geopandas as gpd
from shapely.geometry import LineString, MultiLineString

# 两条部分重叠的线
line_a = LineString([(0, 0), (1, 0), (2, 0), (3, 0)])
line_b = LineString([(1, 0), (2, 0), (3, 0), (4, 0)])

gs = gpd.GeoSeries([line_a])
result = gs.shared_paths(line_b)

print(f"共享路径: {result.iloc[0]}")
# GeometryCollection，包含同方向和反方向的共享段

from shapely.geometry import LineString
from shapely import shared_paths

line_a = LineString([(0, 0), (1, 0), (2, 0)])
line_b = LineString([(1, 0), (2, 0), (3, 0)])

result = shared_paths(line_a, line_b)
# result 是 GeometryCollection
# result.geoms[0]: 同方向共享段（MultiLineString）
# result.geoms[1]: 反方向共享段（MultiLineString）

print(f"同方向共享: {result.geoms[0]}")
print(f"反方向共享: {result.geoms[1]}")

# 找出两条道路的共享路段
road_a = gpd.read_file("road_a.shp").geometry.iloc[0]
road_b = gpd.read_file("road_b.shp").geometry.iloc[0]

gs = gpd.GeoSeries([road_a])
shared = gs.shared_paths(road_b)

if not shared.iloc[0].is_empty:
    print("两条道路有共享路段")
    # 共享路段可能是同一条物理道路被两个数据源重复记录

import geopandas as gpd
from shapely.geometry import Point, Polygon, LineString

# 点到多边形的最短连线
point = Point(5, 5)
poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])

gs = gpd.GeoSeries([poly])
shortest = gs.shortest_line(point)

print(f"最短连线: {shortest.iloc[0]}")
print(f"最短距离: {shortest.iloc[0].length:.4f}")  # 等于 distance()

# 两条不相交的线之间的最短连线
line_a = LineString([(0, 0), (2, 0)])
line_b = LineString([(1, 1), (3, 1)])

gs = gpd.GeoSeries([line_a])
shortest = gs.shortest_line(line_b)
print(f"线到线最短连线: {shortest.iloc[0]}")
# LINESTRING (1 0, 1 1) - 垂直于两条线的连线

# 为每个消防站找到最近的河流取水点
stations = gpd.read_file("fire_stations.shp")
rivers = gpd.read_file("rivers.shp")

# 合并所有河流为一个几何
all_rivers = rivers.union_all()

# 计算每个消防站到最近河流的最短连线
shortest_lines = stations.geometry.apply(
    lambda pt: pt.shortest_line(all_rivers) if hasattr(pt, 'shortest_line')
    else gpd.GeoSeries([pt]).shortest_line(all_rivers).iloc[0]
)

# 获取取水点（最短连线的终点）
water_access_points = shortest_lines.apply(
    lambda line: Point(line.coords[-1])
)

import geopandas as gpd
from shapely.geometry import Point, LineString

# 创建线段
line = LineString([(0, 0), (10, 0)])

# 要投影的点
points = gpd.GeoSeries([
    Point(3, 5),    # 在线上方
    Point(7, -3),   # 在线下方
    Point(-2, 1),   # 在线起点之前
    Point(12, 2),   # 在线终点之后
])

gs_line = gpd.GeoSeries([line] * 4)

# 计算沿线距离
distances = gs_line.project(points)
print(distances)
# 0     3.0  (投影到 (3, 0)，沿线距离 3)
# 1     7.0  (投影到 (7, 0)，沿线距离 7)
# 2     0.0  (投影到起点 (0, 0)，限制为 0)
# 3    10.0  (投影到终点 (10, 0)，限制为 10)

# 归一化距离（0-1范围）
normalized = gs_line.project(points, normalized=True)
print(normalized)
# 0    0.3
# 1    0.7
# 2    0.0
# 3    1.0

import geopandas as gpd
from shapely.geometry import LineString

# 创建线段
line = LineString([(0, 0), (5, 0), (5, 5)])

gs = gpd.GeoSeries([line, line, line])
distances = gpd.array.GeometryDtype  # 不同距离

# 沿线不同距离处插值
result = gs.interpolate([2.0, 5.0, 8.0])
print(result)
# 0    POINT (2 0)      距离 2: 在第一段上
# 1    POINT (5 0)      距离 5: 在转折点
# 2    POINT (5 3)      距离 8: 在第二段上

# 使用归一化距离
total_length = line.length  # 10.0
result_norm = gs.interpolate([0.2, 0.5, 0.8], normalized=True)
print(result_norm)

from shapely.geometry import LineString, Point

line = LineString([(0, 0), (5, 0), (5, 5)])

# 从点投影得到距离
point = Point(3, 2)
distance = line.project(point)
print(f"投影距离: {distance}")

# 从距离反插值得到线上的点
projected_point = line.interpolate(distance)
print(f"线上投影点: {projected_point}")
# 这是 point 在 line 上的最近点

import geopandas as gpd
from shapely.geometry import LineString, Point

# 道路中心线
road = LineString([(0, 0), (100, 0), (100, 50), (200, 50)])

# 每隔 50 米设置一个里程标记
gs = gpd.GeoSeries([road])
total_length = road.length
markers = []
distance = 0

while distance <= total_length:
    point = road.interpolate(distance)
    markers.append({
        'distance': distance,
        'geometry': point
    })
    distance += 50

markers_gdf = gpd.GeoDataFrame(markers)
print(f"共设置 {len(markers_gdf)} 个里程标记")

import geopandas as gpd
from shapely.geometry import Polygon

# 不同索引的 GeoSeries
gs1 = gpd.GeoSeries(
    [Polygon([(0,0),(2,0),(2,2),(0,2)])],
    index=[0]
)
gs2 = gpd.GeoSeries(
    [Polygon([(1,1),(3,1),(3,3),(1,3)])],
    index=[1]
)

# align=True（默认）
result_aligned = gs1.intersection(gs2, align=True)
print(f"对齐后结果长度: {len(result_aligned)}")
# 索引 0 和 1 都会出现，不匹配的变为 NaN

# align=False
result_positional = gs1.intersection(gs2, align=False)
print(f"按位置结果: {result_positional.iloc[0]}")
# 按位置配对：gs1[位置0] ∩ gs2[位置0]

import geopandas as gpd
from shapely.geometry import Polygon

# 而不是对所有要素进行集合运算，先用空间索引筛选
gdf = gpd.read_file("large_dataset.shp")
clip_area = Polygon([(...)]) 

# 先筛选可能相交的要素
possible_idx = gdf.sindex.query(clip_area, predicate='intersects')
candidates = gdf.iloc[possible_idx]

# 只对候选要素进行精确运算
result = candidates.intersection(clip_area)

import numpy as np

# 对于大型 unary_union，可以分批处理
def batch_union(geometries, batch_size=1000):
    """分批合并几何对象"""
    batches = [geometries[i:i+batch_size]
               for i in range(0, len(geometries), batch_size)]

    partial_unions = [gpd.GeoSeries(batch).union_all()
                      for batch in batches]

    return gpd.GeoSeries(partial_unions).union_all()

import geopandas as gpd

gdf1 = gpd.read_file("layer1.shp")
gdf2 = gpd.read_file("layer2.shp")

# overlay 支持的运算类型
result_union = gpd.overlay(gdf1, gdf2, how='union')
result_intersection = gpd.overlay(gdf1, gdf2, how='intersection')
result_difference = gpd.overlay(gdf1, gdf2, how='difference')
result_sym_diff = gpd.overlay(gdf1, gdf2, how='symmetric_difference')
result_identity = gpd.overlay(gdf1, gdf2, how='identity')