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

# 计算 DE-9IM 矩阵
matrix = poly_a.relate(poly_b)
print(f"DE-9IM 矩阵字符串: {matrix}")
# 输出: 212101212

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

# 创建 GeoDataFrame
gdf = gpd.GeoDataFrame(geometry=[
    Polygon([(0, 0), (2, 0), (2, 2), (0, 2)]),
    Polygon([(3, 3), (5, 3), (5, 5), (3, 5)]),
    Polygon([(1, 1), (3, 1), (3, 3), (1, 3)]),
])

# 创建查询几何
query_point = Point(1, 1)

# 判断每个几何对象是否与查询点相交
result = gdf.intersects(query_point)
print(result)
# 0     True
# 1    False
# 2     True
# dtype: bool

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

gs2 = gpd.GeoSeries([
    Point(0.5, 0.5),        # 在第一个多边形内部
    Point(4, 4),             # 不在第二个多边形内
    LineString([(4, 4), (6, 6)]),  # 与第三个点相交
])

# 逐元素判断
result = gs1.intersects(gs2)
print(result)
# 0     True
# 1    False
# 2     True
# dtype: bool

import geopandas as gpd
from shapely.geometry import box

# 读取数据
cities = gpd.read_file("cities.shp")

# 定义查询区域
query_area = box(100, 20, 120, 40)  # 经度100-120, 纬度20-40

# 查找与查询区域相交的城市
intersecting_cities = cities[cities.intersects(query_area)]
print(f"与查询区域相交的城市数量: {len(intersecting_cities)}")

from shapely.geometry import Point, Polygon

# 创建几何对象
poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
point_inside = Point(1, 1)    # 在多边形内部
point_outside = Point(5, 5)   # 在多边形外部
point_on_boundary = Point(2, 1)  # 在多边形边界上

gdf = gpd.GeoDataFrame(geometry=[poly, poly, poly])

# 分别判断
print(gpd.GeoSeries([poly]).disjoint(point_inside)[0])    # False - 点在内部，不分离
print(gpd.GeoSeries([poly]).disjoint(point_outside)[0])   # True  - 完全分离
print(gpd.GeoSeries([poly]).disjoint(point_on_boundary)[0])  # False - 接触边界，不分离

# 查找不在禁区内的设施
facilities = gpd.read_file("facilities.shp")
restricted_area = Polygon([(...)])  # 禁区多边形

# 筛选完全在禁区外的设施
safe_facilities = facilities[facilities.disjoint(restricted_area)]

from shapely.geometry import Polygon, Point, LineString

# 两个相邻但不重叠的多边形
poly_a = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
poly_b = Polygon([(1, 0), (2, 0), (2, 1), (1, 1)])  # 共享边 (1,0)-(1,1)

gs = gpd.GeoSeries([poly_a])
print(gs.touches(poly_b)[0])  # True - 边界接触

# 两个重叠的多边形
poly_c = Polygon([(0.5, 0), (1.5, 0), (1.5, 1), (0.5, 1)])
print(gs.touches(poly_c)[0])  # False - 内部有交集

# 点在多边形边界上
point_on_edge = Point(1, 0.5)
print(gs.touches(point_on_edge)[0])  # True - 点在边界上

# 点在多边形内部
point_inside = Point(0.5, 0.5)
print(gs.touches(point_inside)[0])  # False - 点在内部

# 查找与某省相邻的所有省份
provinces = gpd.read_file("provinces.shp")
target_province = provinces[provinces['name'] == '湖北省'].geometry.iloc[0]

# 使用 touches 找到所有相邻省份
adjacent = provinces[provinces.touches(target_province)]
print(f"与湖北省相邻的省份: {adjacent['name'].tolist()}")

from shapely.geometry import Polygon, Point

# 创建几何对象
outer_poly = Polygon([(0, 0), (4, 0), (4, 4), (0, 4)])
inner_poly = Polygon([(1, 1), (3, 1), (3, 3), (1, 3)])
point_inside = Point(2, 2)
point_on_boundary = Point(0, 0)
point_outside = Point(5, 5)

gs = gpd.GeoSeries([outer_poly])

# contains 测试
print(f"外部多边形包含内部多边形: {gs.contains(inner_poly)[0]}")      # True
print(f"外部多边形包含内部点: {gs.contains(point_inside)[0]}")         # True
print(f"外部多边形包含边界点: {gs.contains(point_on_boundary)[0]}")    # True（注意！）
print(f"外部多边形包含外部点: {gs.contains(point_outside)[0]}")       # False

# within 测试
gs_inner = gpd.GeoSeries([inner_poly])
print(f"内部多边形在外部多边形内: {gs_inner.within(outer_poly)[0]}")    # True

from shapely.geometry import Polygon, Point

poly = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])

# 边界上的点
point_on_edge = Point(0.5, 0)  # 在底边中点

gs = gpd.GeoSeries([poly])
print(f"contains 边界点: {gs.contains(point_on_edge)[0]}")  # True（边界也算包含）

# 判断哪些点落在哪个行政区内
points = gpd.read_file("sample_points.shp")
districts = gpd.read_file("districts.shp")

# 方法一：逐个判断
for idx, district in districts.iterrows():
    mask = points.within(district.geometry)
    points_in_district = points[mask]
    print(f"{district['name']}: {len(points_in_district)} 个点")

# 方法二：使用空间连接（推荐，更高效）
result = gpd.sjoin(points, districts, predicate='within')

from shapely.geometry import Polygon, Point, LineString

poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])

# 内部点
point_inside = Point(1, 1)

# 边界上的点
point_on_boundary = Point(0, 1)

# 接触边界的线
line_touching = LineString([(0.5, 0.5), (2, 1)])  # 一端在内部，一端在边界上

gs = gpd.GeoSeries([poly])

print(f"contains(内部点): {gs.contains(point_inside)[0]}")                # True
print(f"contains_properly(内部点): {gs.contains_properly(point_inside)[0]}")  # True

print(f"contains(边界点): {gs.contains(point_on_boundary)[0]}")               # True
print(f"contains_properly(边界点): {gs.contains_properly(point_on_boundary)[0]}")  # False

print(f"contains(接触边界的线): {gs.contains(line_touching)[0]}")              # True
print(f"contains_properly(接触边界的线): {gs.contains_properly(line_touching)[0]}")  # False

# 只选择严格位于保护区内部的监测站
stations = gpd.read_file("monitoring_stations.shp")
reserve = gpd.read_file("nature_reserve.shp").geometry.iloc[0]

# 严格在内部的站点
strict_inside = stations[gpd.GeoSeries([reserve] * len(stations), index=stations.index)
                         .contains_properly(stations.geometry)]

from shapely.geometry import LineString

# 十字交叉的两条线
line_a = LineString([(0, 0), (2, 2)])
line_b = LineString([(0, 2), (2, 0)])

gs = gpd.GeoSeries([line_a])
print(f"两条线交叉: {gs.crosses(line_b)[0]}")  # True

# 不交叉的两条平行线
line_c = LineString([(0, 1), (2, 3)])
print(f"平行线交叉: {gs.crosses(line_c)[0]}")  # False

# 部分重叠的线
line_d = LineString([(1, 1), (3, 3)])
print(f"重叠线交叉: {gs.crosses(line_d)[0]}")  # False（重叠不算交叉）

from shapely.geometry import LineString, Polygon

poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])

# 穿过多边形的线
line_through = LineString([(-1, 1), (3, 1)])  # 线穿过多边形
gs = gpd.GeoSeries([line_through])
print(f"线穿过多边形: {gs.crosses(poly)[0]}")  # True

# 完全在多边形内的线
line_inside = LineString([(0.5, 0.5), (1.5, 1.5)])
gs2 = gpd.GeoSeries([line_inside])
print(f"线在多边形内: {gs2.crosses(poly)[0]}")  # False（完全包含不算交叉）

from shapely.geometry import Polygon, Point, 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])
print(f"部分重叠: {gs.overlaps(poly_b)[0]}")  # True

# 完全包含的多边形（不算重叠）
poly_c = Polygon([(0.5, 0.5), (1.5, 0.5), (1.5, 1.5), (0.5, 1.5)])
print(f"完全包含: {gs.overlaps(poly_c)[0]}")  # False

# 相同的多边形（不算重叠）
print(f"完全相同: {gs.overlaps(poly_a)[0]}")  # False

# 仅接触的多边形（不算重叠）
poly_d = Polygon([(2, 0), (4, 0), (4, 2), (2, 2)])
print(f"仅接触: {gs.overlaps(poly_d)[0]}")  # False

# 点与多边形不能 overlap
poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
point = Point(1, 1)
gs = gpd.GeoSeries([poly])
print(f"点与多边形 overlaps: {gs.overlaps(point)[0]}")  # False

# 两条部分重叠的线可以 overlap
line_a = LineString([(0, 0), (2, 0)])
line_b = LineString([(1, 0), (3, 0)])
gs2 = gpd.GeoSeries([line_a])
print(f"两条线 overlaps: {gs2.overlaps(line_b)[0]}")  # True

from shapely.geometry import Polygon, Point, LineString

poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])

# 内部点
point_inside = Point(1, 1)
# 边界点
point_on_boundary = Point(0, 1)
# 外部点
point_outside = Point(3, 3)

gs = gpd.GeoSeries([poly])

# covers 测试
print(f"covers 内部点: {gs.covers(point_inside)[0]}")       # True
print(f"covers 边界点: {gs.covers(point_on_boundary)[0]}")  # True
print(f"covers 外部点: {gs.covers(point_outside)[0]}")      # False

# covered_by 测试
gs_point = gpd.GeoSeries([point_inside])
print(f"内部点 covered_by 多边形: {gs_point.covered_by(poly)[0]}")  # True

from shapely.geometry import LineString, Point

# 线段
line = LineString([(0, 0), (1, 0)])
# 线段的端点
endpoint = Point(0, 0)

gs_line = gpd.GeoSeries([line])

# 线段 covers 其端点
print(f"线 covers 端点: {gs_line.covers(endpoint)[0]}")      # True
# 线段 contains 其端点
print(f"线 contains 端点: {gs_line.contains(endpoint)[0]}")  # True（注意：Shapely 2.0 行为）

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

# 创建 GeoSeries
gs = gpd.GeoSeries([
    Point(0, 0),
    Point(3, 0),
    Point(5, 0),
    Point(10, 0),
])

# 查询点
query = Point(0, 0)

# 判断哪些点在距离 4 以内
result = gs.dwithin(query, 4)
print(result)
# 0     True   (距离 0)
# 1     True   (距离 3)
# 2    False   (距离 5)
# 3    False   (距离 10)
# dtype: bool

# 查找某个位置 5 公里范围内的所有学校
schools = gpd.read_file("schools.shp")
# 注意：确保使用投影坐标系（单位为米）
schools_proj = schools.to_crs(epsg=32650)  # UTM Zone 50N

target = Point(500000, 3500000)  # 投影坐标
nearby = schools_proj[schools_proj.dwithin(target, 5000)]  # 5000 米
print(f"5公里范围内的学校: {len(nearby)} 所")

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

# 点到点的距离
gs = gpd.GeoSeries([Point(0, 0), Point(3, 0)])
other = Point(0, 4)
print(gs.distance(other))
# 0    4.0
# 1    5.0
# dtype: float64

# 点到多边形的距离
poly = Polygon([(2, 2), (4, 2), (4, 4), (2, 4)])
gs2 = gpd.GeoSeries([Point(0, 0)])
print(gs2.distance(poly))
# 0    2.828427  (约 2√2)
# dtype: float64

# 点在多边形内部，距离为 0
gs3 = gpd.GeoSeries([Point(3, 3)])
print(gs3.distance(poly))
# 0    0.0
# dtype: float64

from shapely.geometry import LineString
import geopandas as gpd

# 两条相似但不完全相同的线
line_a = LineString([(0, 0), (1, 0), (2, 0)])
line_b = LineString([(0, 0.5), (1, 0.3), (2, 0.5)])

gs = gpd.GeoSeries([line_a])
print(f"Hausdorff 距离: {gs.hausdorff_distance(line_b)[0]:.4f}")
# 输出: Hausdorff 距离: 0.5000

from shapely.geometry import LineString
import geopandas as gpd

# 两条路径
path_a = LineString([(0, 0), (1, 1), (2, 0)])
path_b = LineString([(0, 0), (1, 2), (2, 0)])

gs = gpd.GeoSeries([path_a])
print(f"Fréchet 距离: {gs.frechet_distance(path_b)[0]:.4f}")

# 正确的距离计算流程
import geopandas as gpd

# 1. 读取数据（通常为 WGS84）
points = gpd.read_file("points.shp")

# 2. 投影到合适的坐标系
points_proj = points.to_crs(epsg=32650)

# 3. 计算距离（单位为米）
distances = points_proj.distance(Point(500000, 3500000))
print(f"平均距离: {distances.mean():.2f} 米")

from shapely.geometry import Polygon
import geopandas as gpd

# 相同的多边形，但顶点起始位置不同
poly_a = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
poly_b = Polygon([(1, 0), (1, 1), (0, 1), (0, 0)])  # 不同的起始点

gs = gpd.GeoSeries([poly_a])
print(f"拓扑相等: {gs.geom_equals(poly_b)[0]}")  # True

from shapely.geometry import Point
import geopandas as gpd

point_a = Point(0, 0)
point_b = Point(0.0000001, 0.0000001)  # 微小差异

gs = gpd.GeoSeries([point_a])

# 精确比较
print(f"精确相等: {gs.geom_equals(point_b)[0]}")  # False

# 带容差比较
print(f"近似相等(容差1e-6): {gs.geom_equals_exact(point_b, tolerance=1e-6)[0]}")  # True
print(f"近似相等(容差1e-8): {gs.geom_equals_exact(point_b, tolerance=1e-8)[0]}")  # False

from shapely.geometry import Polygon
import geopandas as gpd

poly_a = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
poly_b = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])  # 完全相同
poly_c = Polygon([(1, 0), (1, 1), (0, 1), (0, 0)])  # 起始点不同

gs = gpd.GeoSeries([poly_a])

print(f"identical(完全相同): {gs.geom_equals_identical(poly_b)[0]}")  # True
print(f"identical(起始点不同): {gs.geom_equals_identical(poly_c)[0]}")  # False

# 对比 geom_equals
print(f"geom_equals(起始点不同): {gs.geom_equals(poly_c)[0]}")  # True

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

# 创建几何对象
poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])

# 内部点
point_inside = Point(1, 1)
gs = gpd.GeoSeries([poly])
matrix = gs.relate(point_inside)
print(f"内部点的 DE-9IM: {matrix[0]}")  # 0F2FF1FF2

# 边界点
point_boundary = Point(0, 0)
matrix = gs.relate(point_boundary)
print(f"边界点的 DE-9IM: {matrix[0]}")  # FF20F1FF2

# 外部点
point_outside = Point(5, 5)
matrix = gs.relate(point_outside)
print(f"外部点的 DE-9IM: {matrix[0]}")  # FF2FF10F2

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

poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
point = Point(1, 1)

gs = gpd.GeoSeries([poly])

# 使用 contains 的 DE-9IM 模式
contains_pattern = "T*****FF*"
print(f"matches contains: {gs.relate_pattern(point, contains_pattern)[0]}")  # True

# 使用 within 的 DE-9IM 模式
within_pattern = "T*F**F***"
gs_point = gpd.GeoSeries([point])
print(f"matches within: {gs_point.relate_pattern(poly, within_pattern)[0]}")  # True

# 自定义模式：内部与内部相交，维度为 2（面）
custom_pattern = "2********"
poly2 = Polygon([(1, 1), (3, 1), (3, 3), (1, 3)])
print(f"内部相交为面: {gs.relate_pattern(poly2, custom_pattern)[0]}")  # True

# 例如：判断两个多边形是否共享一条边（而非仅一个点）
# 条件：边界的交集维度为 1（线）
shared_edge_pattern = "****1****"

poly_a = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
poly_b = Polygon([(1, 0), (2, 0), (2, 1), (1, 1)])  # 共享边

gs = gpd.GeoSeries([poly_a])
print(f"共享边: {gs.relate_pattern(poly_b, shared_edge_pattern)[0]}")  # True

# 仅共享一个点的情况
poly_c = Polygon([(1, 1), (2, 1), (2, 2), (1, 2)])  # 仅共享点(1,1)
print(f"仅共享点: {gs.relate_pattern(poly_c, shared_edge_pattern)[0]}")  # False

# align=True（默认）：根据索引对齐后逐元素比较
# align=False：按位置逐元素比较，忽略索引

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

# 创建两个索引不同的 GeoSeries
gs1 = gpd.GeoSeries(
    [Polygon([(0,0),(1,0),(1,1),(0,1)]), Polygon([(2,2),(3,2),(3,3),(2,3)])],
    index=[0, 1]
)
gs2 = gpd.GeoSeries(
    [Point(2.5, 2.5), Point(0.5, 0.5)],
    index=[1, 0]
)

# align=True：按索引对齐
result = gs1.intersects(gs2, align=True)
print(result)
# 0    True   (gs1[0] vs gs2[0]，即第一个多边形 vs Point(0.5,0.5))
# 1    True   (gs1[1] vs gs2[1]，即第二个多边形 vs Point(2.5,2.5))

# align=False：按位置对齐
result = gs1.intersects(gs2, align=False)
print(result)
# 0    False  (gs1位置0 vs gs2位置0，即第一个多边形 vs Point(2.5,2.5))
# 1    False  (gs1位置1 vs gs2位置1，即第二个多边形 vs Point(0.5,0.5))

import geopandas as gpd
from shapely.geometry import Point
import pandas as pd

# 模拟城市数据
# 创建医院位置
hospitals = gpd.GeoDataFrame({
    'name': ['第一医院', '第二医院', '第三医院', '第四医院'],
    'geometry': [
        Point(116.40, 39.90),
        Point(116.45, 39.92),
        Point(116.38, 39.88),
        Point(116.42, 39.95),
    ]
}, crs='EPSG:4326')

# 创建居民区位置
residential = gpd.GeoDataFrame({
    'name': ['小区A', '小区B', '小区C', '小区D', '小区E'],
    'geometry': [
        Point(116.41, 39.91),
        Point(116.46, 39.93),
        Point(116.39, 39.87),
        Point(116.50, 39.98),
        Point(116.35, 39.85),
    ]
}, crs='EPSG:4326')

# 转换到投影坐标系（米为单位）
hospitals_proj = hospitals.to_crs(epsg=32650)
residential_proj = residential.to_crs(epsg=32650)

# 分析每个医院 3000 米服务范围
service_radius = 3000  # 米
hospital_buffers = hospitals_proj.copy()
hospital_buffers['geometry'] = hospitals_proj.buffer(service_radius)

# 判断每个小区是否在某个医院的服务范围内
for idx, community in residential_proj.iterrows():
    covered = hospital_buffers[hospital_buffers.contains(community.geometry)]
    if len(covered) > 0:
        print(f"{residential.iloc[idx]['name']} 被 {covered['name'].tolist()} 覆盖")
    else:
        print(f"{residential.iloc[idx]['name']} 未被任何医院覆盖！")

import geopandas as gpd

# 读取数据
rivers = gpd.read_file("rivers.shp")
roads = gpd.read_file("roads.shp")

# 确保坐标系一致
roads = roads.to_crs(rivers.crs)

# 找到所有与河流交叉的道路
crossing_pairs = []
for r_idx, river in rivers.iterrows():
    crossing_roads = roads[roads.crosses(river.geometry)]
    for rd_idx, road in crossing_roads.iterrows():
        crossing_pairs.append({
            'river_name': river['name'],
            'road_name': road['name'],
            'intersection': river.geometry.intersection(road.geometry)
        })

crossing_df = pd.DataFrame(crossing_pairs)
print(f"共发现 {len(crossing_df)} 个河流-道路交叉点")
print("这些交叉点可能需要建设桥梁")

import geopandas as gpd
from shapely.geometry import Point

# 综合使用多种空间谓词进行复杂筛选
# 场景：选择新建学校的位置

# 条件1：必须在规划区域内（within）
# 条件2：不能在已有学校 1000 米范围内（disjoint + buffer）
# 条件3：必须在主干道 500 米范围内（dwithin 或 intersects + buffer）
# 条件4：不能在工厂 200 米范围内（disjoint + buffer）

def find_suitable_sites(candidate_sites, planning_area, existing_schools,
                        main_roads, factories):
    """综合筛选适合建校的位置"""

    # 条件1：在规划区域内
    mask1 = candidate_sites.within(planning_area)

    # 条件2：不在已有学校 1000 米范围内
    school_buffers = existing_schools.buffer(1000).union_all()
    mask2 = candidate_sites.disjoint(school_buffers)

    # 条件3：在主干道 500 米范围内
    road_buffers = main_roads.buffer(500).union_all()
    mask3 = candidate_sites.intersects(road_buffers)

    # 条件4：不在工厂 200 米范围内
    factory_buffers = factories.buffer(200).union_all()
    mask4 = candidate_sites.disjoint(factory_buffers)

    # 综合所有条件
    suitable = candidate_sites[mask1 & mask2 & mask3 & mask4]
    return suitable

import geopandas as gpd
import pandas as pd

# 计算所有区域两两之间的空间关系
regions = gpd.read_file("regions.shp")

n = len(regions)
relation_matrix = pd.DataFrame(index=regions['name'], columns=regions['name'])

for i in range(n):
    for j in range(n):
        if i == j:
            relation_matrix.iloc[i, j] = 'self'
        else:
            geom_i = regions.geometry.iloc[i]
            geom_j = regions.geometry.iloc[j]

            if geom_i.contains(geom_j):
                relation_matrix.iloc[i, j] = 'contains'
            elif geom_i.within(geom_j):
                relation_matrix.iloc[i, j] = 'within'
            elif geom_i.overlaps(geom_j):
                relation_matrix.iloc[i, j] = 'overlaps'
            elif geom_i.touches(geom_j):
                relation_matrix.iloc[i, j] = 'touches'
            elif geom_i.disjoint(geom_j):
                relation_matrix.iloc[i, j] = 'disjoint'
            else:
                relation_matrix.iloc[i, j] = geom_i.relate(geom_j)

print(relation_matrix)