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

# 创建包含不同几何类型的 GeoDataFrame
gdf = gpd.GeoDataFrame(
    {"name": ["点", "线", "面", "多点"]},
    geometry=[
        Point(0, 0),
        LineString([(0, 0), (1, 1), (2, 0)]),
        Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
        MultiPoint([(0, 0), (1, 1)])
    ]
)

# 查看每个要素的几何类型
print(gdf.geom_type)

# 判断是否为特定类型
is_point = gdf.geom_type == "Point"
is_polygon = gdf.geom_type.isin(["Polygon", "MultiPolygon"])
is_line = gdf.geom_type.isin(["LineString", "MultiLineString"])

print(f"点要素: {is_point.sum()} 个")
print(f"面要素: {is_polygon.sum()} 个")
print(f"线要素: {is_line.sum()} 个")

# 按类型筛选
gdf_points = gdf[gdf.geom_type == "Point"]
gdf_polygons = gdf[gdf.geom_type.isin(["Polygon", "MultiPolygon"])]

# 统计各类型数量
type_counts = gdf.geom_type.value_counts()
print("几何类型分布:")
print(type_counts)

# 检查是否为单一类型
if len(gdf.geom_type.unique()) == 1:
    print(f"数据集为单一类型: {gdf.geom_type.unique()[0]}")
else:
    print("数据集包含混合类型")

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

gdf = gpd.GeoDataFrame(
    {"name": ["正方形", "三角形", "点", "线"]},
    geometry=[
        Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
        Polygon([(0, 0), (2, 0), (1, 1)]),
        Point(0, 0),
        LineString([(0, 0), (1, 1)])
    ]
)

print(gdf.area)

import geopandas as gpd
from shapely.geometry import box

# 创建一个近似北京市区范围的矩形
beijing_approx = gpd.GeoDataFrame(
    {"name": ["北京近似"]},
    geometry=[box(116.0, 39.6, 116.8, 40.2)],
    crs="EPSG:4326"
)

# ❌ 错误：在地理坐标系下计算面积
area_degrees = beijing_approx.area.iloc[0]
print(f"地理坐标系面积: {area_degrees:.6f} 平方度（无意义）")

# ✅ 正确：转换到投影坐标系后计算
utm_crs = beijing_approx.estimate_utm_crs()
beijing_utm = beijing_approx.to_crs(utm_crs)
area_sqm = beijing_utm.area.iloc[0]
print(f"投影坐标系面积: {area_sqm:,.0f} 平方米")
print(f"投影坐标系面积: {area_sqm / 1e6:,.2f} 平方公里")

def calculate_area(gdf, unit="sqm"):
    """
    计算面积，自动处理坐标系转换

    参数:
        gdf: GeoDataFrame
        unit: 面积单位 - "sqm"(平方米), "sqkm"(平方公里),
              "ha"(公顷), "mu"(亩)
    """
    # 如果是地理坐标系，先转换到投影坐标系
    if gdf.crs and gdf.crs.is_geographic:
        utm_crs = gdf.estimate_utm_crs()
        gdf_proj = gdf.to_crs(utm_crs)
    else:
        gdf_proj = gdf

    area_sqm = gdf_proj.area

    # 单位转换
    conversions = {
        "sqm": 1,
        "sqkm": 1e-6,
        "ha": 1e-4,
        "mu": 1 / 666.667,  # 1亩 ≈ 666.667 平方米
        "acre": 1 / 4046.86
    }

    factor = conversions.get(unit, 1)
    return area_sqm * factor

# 使用示例
gdf["area_sqkm"] = calculate_area(gdf, unit="sqkm")
gdf["area_mu"] = calculate_area(gdf, unit="mu")
print(gdf[["name", "area_sqkm", "area_mu"]])

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

gdf = gpd.GeoDataFrame(
    {"name": ["线段", "折线", "多边形", "点"]},
    geometry=[
        LineString([(0, 0), (3, 4)]),              # 直线段
        LineString([(0, 0), (1, 0), (1, 1)]),      # 折线
        Polygon([(0, 0), (4, 0), (4, 3), (0, 3)]), # 矩形
        Point(0, 0)                                  # 点
    ]
)

print(gdf.length)

from shapely.geometry import Polygon

# 带洞的多边形
polygon_with_hole = Polygon(
    [(0, 0), (10, 0), (10, 10), (0, 10)],  # 外环
    [[(2, 2), (4, 2), (4, 4), (2, 4)]]     # 内环（洞）
)

# 注意：Polygon 的 length 返回的是外环的周长
print(f"多边形周长（外环）: {polygon_with_hole.length}")
# 40.0

# 如果需要包括内环的总长度
exterior_length = polygon_with_hole.exterior.length
interior_length = sum(ring.length for ring in polygon_with_hole.interiors)
total_length = exterior_length + interior_length
print(f"外环周长: {exterior_length}")
print(f"内环周长: {interior_length}")
print(f"总周长: {total_length}")

# ❌ 地理坐标系下的长度单位是"度"，无意义
gdf_geo = gpd.GeoDataFrame(
    geometry=[LineString([(116.0, 39.0), (117.0, 39.0)])],
    crs="EPSG:4326"
)
print(f"地理坐标系长度: {gdf_geo.length.iloc[0]:.4f} 度")

# ✅ 转换到投影坐标系后计算
gdf_utm = gdf_geo.to_crs(gdf_geo.estimate_utm_crs())
print(f"投影坐标系长度: {gdf_utm.length.iloc[0]:,.0f} 米")
# 约 85,394 米（在纬度 39° 处，1° 经度约 85.4 公里）

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

gdf = gpd.GeoDataFrame(
    {"name": ["三角形", "折线", "点"]},
    geometry=[
        Polygon([(0, 0), (4, 0), (2, 3)]),
        LineString([(1, 1), (5, 3), (3, 5)]),
        Point(2, 2)
    ]
)

print(gdf.bounds)

# 将边界框信息添加到 GeoDataFrame
bounds = gdf.bounds
gdf["bbox_width"] = bounds["maxx"] - bounds["minx"]
gdf["bbox_height"] = bounds["maxy"] - bounds["miny"]
gdf["bbox_area"] = gdf["bbox_width"] * gdf["bbox_height"]

print(gdf[["name", "bbox_width", "bbox_height", "bbox_area"]])

# 所有要素的总范围
total = gdf.total_bounds
print(f"总范围: {total}")
# [minx, miny, maxx, maxy]

print(f"X 范围: {total[0]:.4f} - {total[2]:.4f}")
print(f"Y 范围: {total[1]:.4f} - {total[3]:.4f}")

# 用途示例：创建覆盖所有数据的边界矩形
from shapely.geometry import box
bbox_polygon = box(*total)
print(f"边界矩形: {bbox_polygon}")

import geopandas as gpd

# 读取中国省级数据
# gdf = gpd.read_file("china_provinces.shp")

# 获取每个省份的边界范围
# province_bounds = gdf.bounds
# province_bounds["province"] = gdf["name"]
# print(province_bounds.head())

# 获取中国的整体范围
# china_extent = gdf.total_bounds
# print(f"中国范围:")
# print(f"  经度: {china_extent[0]:.2f}° - {china_extent[2]:.2f}°")
# print(f"  纬度: {china_extent[1]:.2f}° - {china_extent[3]:.2f}°")
# 约为: 经度 73.5° - 135.1°, 纬度 18.2° - 53.6°

# 使用总范围进行空间查询
# other_gdf = gpd.read_file("data.shp", bbox=tuple(gdf.total_bounds))

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

gdf = gpd.GeoDataFrame(
    {"name": ["正方形", "三角形", "线段"]},
    geometry=[
        Polygon([(0, 0), (4, 0), (4, 4), (0, 4)]),
        Polygon([(0, 0), (6, 0), (3, 4)]),
        LineString([(0, 0), (10, 0)])
    ]
)

centroids = gdf.centroid
print("质心坐标:")
for name, centroid in zip(gdf["name"], centroids):
    print(f"  {name}: ({centroid.x:.2f}, {centroid.y:.2f})")

from shapely.geometry import Polygon

# 注意：质心可能不在几何内部！
# 例如凹多边形（L 形）
l_shape = Polygon([
    (0, 0), (2, 0), (2, 1), (1, 1),
    (1, 3), (0, 3)
])

centroid = l_shape.centroid
print(f"L 形质心: ({centroid.x:.2f}, {centroid.y:.2f})")
print(f"质心在几何内部: {l_shape.contains(centroid)}")
# 质心可能落在 L 形的凹陷处，不在几何内部

# 常见用途：在地图上为面要素添加标注点
gdf["centroid_x"] = gdf.centroid.x
gdf["centroid_y"] = gdf.centroid.y

# 创建质心点图层
gdf_centroids = gdf.copy()
gdf_centroids = gdf_centroids.set_geometry(gdf.centroid)
print(gdf_centroids.geom_type.unique())
# ['Point']

import geopandas as gpd
from shapely.geometry import Polygon

# 创建 L 形（凹多边形）
l_shape = Polygon([
    (0, 0), (4, 0), (4, 2), (2, 2),
    (2, 4), (0, 4)
])

gdf = gpd.GeoDataFrame(
    {"name": ["L形"]},
    geometry=[l_shape]
)

# 对比质心和代表点
centroid = gdf.centroid.iloc[0]
rep_point = gdf.representative_point().iloc[0]

print(f"质心: ({centroid.x:.2f}, {centroid.y:.2f})")
print(f"代表点: ({rep_point.x:.2f}, {rep_point.y:.2f})")
print(f"质心在内部: {l_shape.contains(centroid)}")
print(f"代表点在内部: {l_shape.contains(rep_point)}")

# 实际建议：
# - 统计分析（如加权平均中心）→ 使用 centroid
# - 地图标注 → 使用 representative_point()
# - 空间索引点 → 使用 representative_point()

# 为每个面要素生成标注点
gdf["label_point"] = gdf.representative_point()

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

# 创建包含空几何的数据
gdf = gpd.GeoDataFrame(
    {"name": ["正常点", "空几何", "正常线"]},
    geometry=[
        Point(1, 2),
        wkt.loads("GEOMETRYCOLLECTION EMPTY"),
        LineString([(0, 0), (1, 1)])
    ]
)

print(gdf.is_empty)

# 过滤掉空几何
gdf_valid = gdf[~gdf.is_empty]
print(f"移除空几何后: {len(gdf_valid)} 条记录")

# 统计空几何数量
empty_count = gdf.is_empty.sum()
print(f"空几何数量: {empty_count}")

# 替换空几何
# gdf.loc[gdf.is_empty, "geometry"] = Point(0, 0)  # 用默认点替代

import geopandas as gpd
from shapely.geometry import Polygon

# 创建有效和无效的几何
valid_polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])

# "蝴蝶结"（自相交多边形）— 无效
invalid_polygon = Polygon([(0, 0), (1, 1), (1, 0), (0, 1)])

gdf = gpd.GeoDataFrame(
    {"name": ["有效", "无效"]},
    geometry=[valid_polygon, invalid_polygon]
)

print(gdf.is_valid)

from shapely.validation import explain_validity

# 检查无效原因
for idx, row in gdf.iterrows():
    if not row.geometry.is_valid:
        reason = explain_validity(row.geometry)
        print(f"要素 {row['name']}: {reason}")

# 在 GeoSeries 级别
# GeoPandas 提供了 is_valid_reason 属性（取决于版本）
# print(gdf.geometry.is_valid_reason)

from shapely.validation import make_valid

# 方法1：使用 make_valid()
gdf["geometry"] = gdf.geometry.make_valid()

# 方法2：使用 buffer(0) 修复（经典方法）
gdf["geometry"] = gdf.geometry.buffer(0)

# 方法3：逐个检查和修复
def fix_geometry(geom):
    """尝试修复无效几何"""
    if geom is None or geom.is_empty:
        return geom
    if not geom.is_valid:
        # 先尝试 make_valid
        fixed = make_valid(geom)
        if fixed.is_valid:
            return fixed
        # 再尝试 buffer(0)
        fixed = geom.buffer(0)
        return fixed
    return geom

gdf["geometry"] = gdf.geometry.apply(fix_geometry)

# 验证修复结果
print(f"所有几何有效: {gdf.is_valid.all()}")

import geopandas as gpd
from shapely.geometry import LineString

# 简单线（不自相交）
simple_line = LineString([(0, 0), (1, 1), (2, 0)])

# 非简单线（自相交，形成8字形）
complex_line = LineString([(0, 0), (1, 1), (1, 0), (0, 1)])

gdf = gpd.GeoDataFrame(
    {"name": ["简单线", "自相交线"]},
    geometry=[simple_line, complex_line]
)

print(gdf.is_simple)

# 对于多边形：有效则一定简单
# 对于线：简单不一定有效（线没有有效性概念）
# 点始终是简单的

import geopandas as gpd
from shapely.geometry import LineString

# 闭合环（首尾相连）
ring = LineString([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])

# 非闭合线
open_line = LineString([(0, 0), (1, 0), (1, 1)])

# 自相交的闭合线（不是环，因为不简单）
self_intersecting = LineString([(0, 0), (1, 1), (1, 0), (0, 1), (0, 0)])

gdf = gpd.GeoDataFrame(
    {"name": ["闭合环", "非闭合", "自相交闭合"]},
    geometry=[ring, open_line, self_intersecting]
)

print(gdf.is_ring)

# is_closed 只检查首尾是否相连
print("is_closed:")
for name, geom in zip(gdf["name"], gdf.geometry):
    if hasattr(geom, 'is_closed'):
        print(f"  {name}: {geom.is_closed}")

# is_ring = is_closed AND is_simple

import geopandas as gpd
from shapely.geometry import Point

# 二维点
point_2d = Point(116.4, 39.9)

# 三维点（含高程）
point_3d = Point(116.4, 39.9, 50.0)

gdf = gpd.GeoDataFrame(
    {"name": ["2D点", "3D点"]},
    geometry=[point_2d, point_3d]
)

print(gdf.has_z)

# Z 值常用于：
# - 地形高程数据
# - 建筑物高度
# - 地下管线深度
# - 三维城市模型

# 获取 Z 值
point_3d = Point(116.4, 39.9, 50.0)
print(f"X: {point_3d.x}")
print(f"Y: {point_3d.y}")
print(f"Z: {point_3d.z}")

# 注意：大多数二维空间分析会忽略 Z 值
# 面积、距离等计算通常只使用 X、Y 坐标

# 注意：Shapely 对 M 值的支持有限
# has_m 属性在较新版本的 Shapely 中可用

# M 值常用于：
# - 道路里程桩
# - 管线测量
# - 时间序列轨迹

# 检查是否支持
# print(gdf.has_m)  # 取决于 Shapely 版本

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

# 创建不同复杂度的几何
gdf = gpd.GeoDataFrame(
    {"name": ["点", "线(3点)", "三角形", "矩形"]},
    geometry=[
        Point(0, 0),                                    # 1 个坐标
        LineString([(0, 0), (1, 1), (2, 0)]),           # 3 个坐标
        Polygon([(0, 0), (1, 0), (0, 1)]),              # 4 个坐标（闭合）
        Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])       # 5 个坐标（闭合）
    ]
)

# 使用 shapely 的 count_coordinates 函数
coord_counts = shapely.count_coordinates(gdf.geometry)
print("坐标数量:")
for name, count in zip(gdf["name"], coord_counts):
    print(f"  {name}: {count}")

from shapely.geometry import MultiPoint, MultiPolygon, Point, Polygon
import shapely

gdf = gpd.GeoDataFrame(
    {"name": ["单点", "3个多点", "2个多面"]},
    geometry=[
        Point(0, 0),
        MultiPoint([(0, 0), (1, 1), (2, 2)]),
        MultiPolygon([
            Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
            Polygon([(2, 2), (3, 2), (3, 3), (2, 3)])
        ])
    ]
)

geom_counts = shapely.count_geometries(gdf.geometry)
print("子几何数量:")
for name, count in zip(gdf["name"], geom_counts):
    print(f"  {name}: {count}")

from shapely.geometry import Polygon
import shapely

# 无洞多边形
simple_poly = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])

# 有 1 个洞的多边形
one_hole = Polygon(
    [(0, 0), (10, 0), (10, 10), (0, 10)],
    [[(2, 2), (4, 2), (4, 4), (2, 4)]]
)

# 有 2 个洞的多边形
two_holes = Polygon(
    [(0, 0), (10, 0), (10, 10), (0, 10)],
    [
        [(1, 1), (3, 1), (3, 3), (1, 3)],
        [(5, 5), (7, 5), (7, 7), (5, 7)]
    ]
)

gdf = gpd.GeoDataFrame(
    {"name": ["无洞", "1个洞", "2个洞"]},
    geometry=[simple_poly, one_hole, two_holes]
)

ring_counts = shapely.count_interior_rings(gdf.geometry)
print("内环数量:")
for name, count in zip(gdf["name"], ring_counts):
    print(f"  {name}: {count}")

import shapely
from shapely.geometry import Point

# 默认精度为 0（无限精度/浮点精度）
point = Point(116.40723456789, 39.90456789123)
precision = shapely.get_precision(point)
print(f"默认精度: {precision}")
# 0.0

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

# 降低精度（四舍五入坐标）
point = Point(116.40723456789, 39.90456789123)

# 精度为 0.001（约 100 米级别）
point_low = shapely.set_precision(point, grid_size=0.001)
print(f"原始: ({point.x}, {point.y})")
print(f"降精: ({point_low.x}, {point_low.y})")

# 对整个 GeoSeries 设置精度
gdf = gpd.GeoDataFrame(
    {"name": ["A", "B"]},
    geometry=[
        Point(116.40723456789, 39.90456789123),
        Point(121.47456789012, 31.23056789012)
    ],
    crs="EPSG:4326"
)

# 将坐标精度降低到小数点后 4 位
gdf["geometry"] = shapely.set_precision(gdf.geometry.values, grid_size=0.0001)
print(gdf.geometry)

# 1. 减小文件大小（较低精度 → 更少的存储空间）
# 地理坐标系下，不同精度对应的近似空间精度：
# 0.1      → 约 11 km
# 0.01     → 约 1.1 km
# 0.001    → 约 110 m
# 0.0001   → 约 11 m
# 0.00001  → 约 1.1 m
# 0.000001 → 约 0.11 m

# 2. 避免浮点数精度问题导致的拓扑错误
# 3. 匹配不同数据源的精度级别

import geopandas as gpd
from shapely.geometry import Polygon

# 创建数据（WGS84 地理坐标系）
gdf = gpd.GeoDataFrame(
    {"name": ["区域A", "区域B"]},
    geometry=[
        Polygon([(116.0, 39.5), (116.5, 39.5), (116.5, 40.0), (116.0, 40.0)]),
        Polygon([(121.0, 31.0), (121.5, 31.0), (121.5, 31.5), (121.0, 31.5)])
    ],
    crs="EPSG:4326"
)

# 步骤1：检查当前 CRS
print(f"CRS: {gdf.crs}")
print(f"是地理坐标系: {gdf.crs.is_geographic}")

# 步骤2：转换到投影坐标系
# 方法A：使用 estimate_utm_crs()
utm_crs = gdf.estimate_utm_crs()
gdf_utm = gdf.to_crs(utm_crs)

# 方法B：使用等面积投影（更适合面积计算）
# gdf_aea = gdf.to_crs("+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +datum=WGS84 +units=m")

# 步骤3：计算面积和周长
gdf["area_sqm"] = gdf_utm.area
gdf["area_sqkm"] = gdf_utm.area / 1e6
gdf["perimeter_m"] = gdf_utm.length
gdf["perimeter_km"] = gdf_utm.length / 1000

print(gdf[["name", "area_sqkm", "perimeter_km"]])

import geopandas as gpd
from shapely.geometry import Polygon

# 创建一个近似北京市区的矩形
beijing = gpd.GeoDataFrame(
    {"name": ["北京"]},
    geometry=[Polygon([
        (116.0, 39.6), (116.8, 39.6),
        (116.8, 40.2), (116.0, 40.2)
    ])],
    crs="EPSG:4326"
)

# 使用不同投影计算面积
projections = {
    "UTM 50N": "EPSG:32650",
    "Web Mercator": "EPSG:3857",
    "Albers 等面积": "+proj=aea +lat_1=25 +lat_2=47 +lon_0=105 +datum=WGS84 +units=m"
}

print("不同投影下的面积对比:")
for name, crs in projections.items():
    area = beijing.to_crs(crs).area.iloc[0] / 1e6
    print(f"  {name}: {area:,.2f} 平方公里")

# 建议：
# - 面积计算优先使用等面积投影
# - UTM 在小范围内也足够准确
# - Web Mercator 不适合面积计算（面积变形严重）

# 单位转换工具函数
class AreaConverter:
    """面积单位转换器"""
    FACTORS = {
        "sqm": 1,
        "sqkm": 1e-6,
        "ha": 1e-4,
        "mu": 1 / 666.667,
        "acre": 1 / 4046.86,
        "sqmi": 1 / 2589988
    }

    @staticmethod
    def convert(value_sqm, to_unit):
        """将平方米转换为目标单位"""
        factor = AreaConverter.FACTORS.get(to_unit)
        if factor is None:
            raise ValueError(f"不支持的单位: {to_unit}")
        return value_sqm * factor

class LengthConverter:
    """长度单位转换器"""
    FACTORS = {
        "m": 1,
        "km": 0.001,
        "mi": 1 / 1609.34,
        "nmi": 1 / 1852,
        "ft": 3.28084
    }

    @staticmethod
    def convert(value_m, to_unit):
        factor = LengthConverter.FACTORS.get(to_unit)
        if factor is None:
            raise ValueError(f"不支持的单位: {to_unit}")
        return value_m * factor

# 使用示例
area_sqm = 16400000  # 16.4 平方公里
print(f"面积: {AreaConverter.convert(area_sqm, 'sqkm'):.2f} 平方公里")
print(f"面积: {AreaConverter.convert(area_sqm, 'ha'):.0f} 公顷")
print(f"面积: {AreaConverter.convert(area_sqm, 'mu'):.0f} 亩")