# Shapely 的基本用法一览
from shapely import Point, Polygon

# 创建一个点
point = Point(1.0, 2.0)
print(f"点: {point}")           # POINT (1 2)
print(f"类型: {point.geom_type}")  # Point

# 创建一个多边形
polygon = Polygon([(0, 0), (4, 0), (4, 4), (0, 4)])
print(f"面积: {polygon.area}")     # 16.0
print(f"周长: {polygon.length}")   # 16.0

# 空间关系判断
print(f"点在多边形内: {polygon.contains(point)}")  # True

from shapely import Point

p = Point(1, 2)
# Shapely 2.0 中，几何对象是不可变的
# 不能修改已有几何体，只能创建新的
# p.coords = [(3, 4)]  # 这在 2.0 中会抛出异常

# 正确方式：创建新的几何体
p_new = Point(3, 4)

import numpy as np
import shapely

# 批量创建点——向量化操作，性能远优于循环
coords = np.array([[0, 0], [1, 1], [2, 2], [3, 3]])
points = shapely.points(coords)
print(points)  # 返回 NumPy 数组，包含 4 个 Point 对象

# 批量计算距离
origin = Point(0, 0)
distances = shapely.distance(points, origin)
print(distances)  # [0.  1.41421356  2.82842712  4.24264069]

from shapely import Point

p1 = Point(1, 2)
p2 = Point(3, 4)

# 2.0 中几何对象可用作字典键和集合元素
point_set = {p1, p2}
point_dict = {p1: "A", p2: "B"}
print(point_dict[p1])  # "A"

import numpy as np
import shapely

# NumPy 数组可直接存储几何对象
geom_array = np.array([Point(0, 0), Point(1, 1), Point(2, 2)])

# 所有 shapely 顶级函数都支持数组输入
areas = shapely.area(geom_array)
is_valid = shapely.is_valid(geom_array)

# Shapely 专注于几何计算的示例
from shapely import Point, Polygon

# 创建几何体（不需要指定坐标系）
park = Polygon([(0, 0), (100, 0), (100, 80), (0, 80)])
school = Point(50, 40)

# 几何分析
buffer_zone = school.buffer(30)          # 缓冲区分析
overlap = park.intersection(buffer_zone)  # 交集运算
print(f"重叠面积: {overlap.area:.2f}")    # 几何计算结果

from shapely import Point, Polygon, box

# 像使用普通 Python 对象一样使用几何体
p = Point(1, 2)
print(f"x={p.x}, y={p.y}")        # 属性访问
print(f"WKT: {p}")                 # 字符串表示即 WKT
print(f"repr: {repr(p)}")          # POINT (1 2)

# 运算符重载
poly1 = box(0, 0, 2, 2)
poly2 = box(1, 1, 3, 3)

union = poly1 | poly2              # | 运算符 = union
inter = poly1 & poly2              # & 运算符 = intersection
diff = poly1 - poly2               # - 运算符 = difference
sym_diff = poly1 ^ poly2           # ^ 运算符 = symmetric_difference

# in 运算符
print(Point(1, 1) in poly1)       # 等价于 poly1.contains(Point(1, 1))

import shapely
import numpy as np
from concurrent.futures import ThreadPoolExecutor

# 生成大量几何体
polygons = shapely.box(
    np.random.uniform(0, 100, 1000),
    np.random.uniform(0, 100, 1000),
    np.random.uniform(100, 200, 1000),
    np.random.uniform(100, 200, 1000)
)

# GEOS 操作在多线程中可以并行执行
# 因为 GEOS 计算时释放了 GIL
areas = shapely.area(polygons)       # 向量化操作，内部已并行
buffers = shapely.buffer(polygons, 5.0)

import shapely
# 查看 GEOS 版本
print(f"GEOS 版本: {shapely.geos_version}")         # (3, 11, 2)
print(f"Shapely 版本: {shapely.__version__}")        # 2.0.x

import geopandas as gpd
from shapely import Point

# GeoPandas 内部使用 Shapely 几何对象
gdf = gpd.GeoDataFrame({
    'name': ['北京', '上海', '广州'],
    'geometry': [Point(116.4, 39.9), Point(121.5, 31.2), Point(113.3, 23.1)]
})

# 每个元素都是 Shapely 对象
print(type(gdf.geometry[0]))  # <class 'shapely.geometry.point.Point'>

import fiona
from shapely.geometry import shape

# Fiona 读取，Shapely 分析
with fiona.open("buildings.shp") as src:
    for feature in src:
        geom = shape(feature['geometry'])  # 转为 Shapely 对象
        print(f"面积: {geom.area}")

from shapely import Polygon, MultiPolygon
from shapely.ops import unary_union
from shapely.validation import make_valid

# 修复无效的几何体
invalid_poly = Polygon([(0, 0), (2, 2), (2, 0), (0, 2)])  # 自交叉
valid_poly = make_valid(invalid_poly)
print(f"修复前有效: {invalid_poly.is_valid}")  # False
print(f"修复后有效: {valid_poly.is_valid}")    # True

# 合并多个几何体
polygons = [
    Polygon([(0, 0), (2, 0), (2, 2), (0, 2)]),
    Polygon([(1, 1), (3, 1), (3, 3), (1, 3)]),
    Polygon([(2, 2), (4, 2), (4, 4), (2, 4)])
]
merged = unary_union(polygons)
print(f"合并后类型: {merged.geom_type}")

from shapely import Point, Polygon, LineString

# 点在面内判断（地理围栏）
fence = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
location = Point(5, 5)
print(f"在围栏内: {fence.contains(location)}")  # True

# 路径与区域的交集（道路穿过公园的部分）
road = LineString([(0, 5), (15, 5)])
park = Polygon([(2, 2), (8, 2), (8, 8), (2, 8)])
road_in_park = road.intersection(park)
print(f"穿过公园的路段长度: {road_in_park.length}")  # 6.0

from shapely import Point, box
from shapely.geometry import mapping
import json

# GeoJSON 输出（Web API 返回）
polygon = box(0, 0, 10, 10)
geojson = mapping(polygon)
print(json.dumps(geojson, indent=2))
# {
#   "type": "Polygon",
#   "coordinates": [[[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0]]]
# }

# 空间查询：找到视口范围内的点
viewport = box(100, 30, 120, 40)  # 经度 100-120, 纬度 30-40
cities = [Point(116.4, 39.9), Point(121.5, 31.2), Point(113.3, 23.1)]
visible = [c for c in cities if viewport.contains(c)]
print(f"可见城市数: {len(visible)}")

import numpy as np
import shapely

# 批量生成随机多边形并计算面积分布
np.random.seed(42)
n = 10000
x = np.random.uniform(0, 100, n)
y = np.random.uniform(0, 100, n)
points = shapely.points(x, y)
buffers = shapely.buffer(points, np.random.uniform(1, 10, n))
areas = shapely.area(buffers)

print(f"平均面积: {areas.mean():.2f}")
print(f"最大面积: {areas.max():.2f}")
print(f"面积标准差: {areas.std():.2f}")

from shapely import Point, LineString, Polygon

# 创建点
p1 = Point(0, 0)
p2 = Point(3, 4)
print(f"点 p1: {p1}")               # POINT (0 0)
print(f"点 p2: {p2}")               # POINT (3 4)
print(f"两点距离: {p1.distance(p2)}")  # 5.0

# 创建线
line = LineString([(0, 0), (1, 1), (2, 0), (3, 1)])
print(f"线: {line}")
print(f"线长度: {line.length:.4f}")   # 4.2426

# 创建多边形
poly = Polygon([(0, 0), (4, 0), (4, 3), (0, 3)])
print(f"多边形: {poly}")
print(f"面积: {poly.area}")           # 12.0
print(f"周长: {poly.length}")         # 14.0
print(f"质心: {poly.centroid}")       # POINT (2 1.5)

from shapely import Point

# 以点为圆心创建缓冲区（近似圆）
center = Point(0, 0)
circle = center.buffer(10)
print(f"缓冲区面积: {circle.area:.2f}")  # 约 314.16（≈ π * 10²）
print(f"缓冲区顶点数: {len(circle.exterior.coords)}")  # 默认 65 个顶点

# 控制缓冲区精度
circle_low = center.buffer(10, quad_segs=4)   # 低精度
circle_high = center.buffer(10, quad_segs=64)  # 高精度
print(f"低精度面积: {circle_low.area:.2f}")
print(f"高精度面积: {circle_high.area:.2f}")

from shapely import Polygon

# 两个重叠的矩形
rect1 = Polygon([(0, 0), (3, 0), (3, 3), (0, 3)])
rect2 = Polygon([(1, 1), (4, 1), (4, 4), (1, 4)])

# 交集
intersection = rect1.intersection(rect2)
print(f"交集面积: {intersection.area}")  # 4.0
print(f"交集形状: {intersection}")
# POLYGON ((1 1, 3 1, 3 3, 1 3, 1 1))

# 并集
union = rect1.union(rect2)
print(f"并集面积: {union.area}")  # 14.0

# 差集
difference = rect1.difference(rect2)
print(f"差集面积: {difference.area}")  # 5.0

# 对称差集
sym_diff = rect1.symmetric_difference(rect2)
print(f"对称差集面积: {sym_diff.area}")  # 10.0

from shapely import Point, Polygon

# 定义围栏区域（模拟校园边界）
campus = Polygon([
    (116.32, 39.98), (116.34, 39.98),
    (116.34, 40.00), (116.32, 40.00)
])

# 模拟 GPS 定位点
locations = [
    ("小明", Point(116.33, 39.99)),   # 校园内
    ("小红", Point(116.35, 39.99)),   # 校园外
    ("小刚", Point(116.34, 39.98)),   # 边界上
]

for name, loc in locations:
    if campus.contains(loc):
        status = "在校园内 ✓"
    elif campus.touches(loc) or campus.boundary.distance(loc) < 0.0001:
        status = "在校园边界 ~"
    else:
        status = "在校园外 ✗"
    print(f"{name}: {status}")

# 典型的 Python GIS 工作流
import fiona          # 读取空间数据文件
from shapely.geometry import shape, mapping  # 几何操作
import pyproj         # 坐标转换
import geopandas as gpd  # 空间数据框分析

# 1. 用 Fiona 读取数据
# 2. 用 Shapely 做几何计算
# 3. 用 pyproj 做坐标转换
# 4. 用 GeoPandas 管理数据和批量操作
# 5. 用 matplotlib 可视化结果