class GeoDataFrame(pd.DataFrame):
    """
    一个包含几何列的 pandas DataFrame，
    支持空间操作和地理空间数据处理。
    """
    pass

import geopandas as gpd
from shapely.geometry import Point

# 创建一个 GeoDataFrame
gdf = gpd.GeoDataFrame({
    '名称': ['北京', '上海', '广州'],
    '人口': [2189, 2487, 1868],
    'GDP': [41610, 44652, 28231],
    'geometry': [
        Point(116.40, 39.90),
        Point(121.47, 31.23),
        Point(113.26, 23.13)
    ]
}, crs="EPSG:4326")

# 查看结构
print("类型:", type(gdf))
print("列:", gdf.columns.tolist())
print("几何列名:", gdf.geometry.name)
print("CRS:", gdf.crs)
print("空间索引:", type(gdf.sindex))

import geopandas as gpd
from shapely.geometry import Point

gdf = gpd.GeoDataFrame({
    '名称': ['A', 'B', 'C'],
    'geometry': [Point(0, 0), Point(1, 1), Point(2, 2)],
    'centroid_geom': [Point(0.5, 0.5), Point(1.5, 1.5), Point(2.5, 2.5)]
})

# 当前活跃几何列
print("活跃几何列:", gdf.geometry.name)  # 'geometry'

# 切换活跃几何列
gdf = gdf.set_geometry('centroid_geom')
print("切换后活跃几何列:", gdf.geometry.name)  # 'centroid_geom'

# 空间操作现在作用于 centroid_geom
print(gdf.total_bounds)

class GeoSeries(pd.Series):
    """
    存储几何对象的 pandas Series，
    支持空间操作和属性访问。
    """
    pass

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

# 方法1：从 Shapely 对象列表创建
gs = gpd.GeoSeries([
    Point(0, 0),
    Point(1, 1),
    Point(2, 2)
], crs="EPSG:4326")
print("方法1:", gs)

# 方法2：从 WKT 创建
gs_wkt = gpd.GeoSeries.from_wkt([
    'POINT (116.40 39.90)',
    'POINT (121.47 31.23)',
    'POINT (113.26 23.13)'
], crs="EPSG:4326")
print("\n方法2:", gs_wkt)

# 方法3：从 WKB 创建
wkb_data = gs.to_wkb()
gs_wkb = gpd.GeoSeries.from_wkb(wkb_data, crs="EPSG:4326")
print("\n方法3:", gs_wkb)

# 方法4：从坐标数组创建点
gs_xy = gpd.GeoSeries.from_xy(
    x=[116.40, 121.47, 113.26],
    y=[39.90, 31.23, 23.13],
    crs="EPSG:4326"
)
print("\n方法4:", gs_xy)

# 方法5：包含不同类型的几何对象
gs_mixed = gpd.GeoSeries([
    Point(0, 0),
    LineString([(0, 0), (1, 1)]),
    Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
])
print("\n方法5 (混合类型):", gs_mixed)

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

gs = gpd.GeoSeries([Point(0, 0), Point(1, 1)], crs="EPSG:4326")

# GeoSeries 是 pandas Series 的子类
print(isinstance(gs, pd.Series))     # True
print(isinstance(gs, gpd.GeoSeries)) # True

# pandas Series 的方法都可以使用
print(gs.index)        # RangeIndex
print(len(gs))         # 2
print(gs.is_empty)     # GeoSeries 特有方法
print(gs.dtype)        # geometry

import geopandas as gpd
from shapely.geometry import Polygon

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

# 基本属性
print("几何类型:", gs.geom_type.tolist())    # ['Polygon', 'Polygon']
print("是否为空:", gs.is_empty.tolist())      # [False, False]
print("是否有效:", gs.is_valid.tolist())      # [True, True]
print("是否简单:", gs.is_simple.tolist())     # [True, True]

# 度量属性
print("面积:", gs.area.tolist())              # [4.0, 6.0]
print("长度:", gs.length.tolist())            # [8.0, 10.0]
print("边界框:", gs.bounds)
print("总边界框:", gs.total_bounds)            # [0. 0. 4. 3.]

# 派生几何
print("质心:", gs.centroid.tolist())
print("边界:", gs.boundary.tolist())
print("包络框:", gs.envelope.tolist())
print("凸包:", gs.convex_hull.tolist())

# GeometryArray 的位置
from geopandas.array import GeometryArray

# 它实现了 pandas ExtensionArray 接口
import pandas as pd
print(issubclass(GeometryArray, pd.api.extensions.ExtensionArray))  # True

import geopandas as gpd
from shapely.geometry import Point
import numpy as np

gs = gpd.GeoSeries([Point(0, 0), Point(1, 1), Point(2, 2)])

# 获取底层 GeometryArray
ga = gs.values
print("类型:", type(ga))                    # GeometryArray
print("底层数据类型:", type(ga._data))       # numpy.ndarray
print("元素类型:", type(ga._data[0]))        # shapely.geometry.point.Point

# GeometryArray 的 dtype
print("dtype:", ga.dtype)                    # geometry
print("dtype name:", ga.dtype.name)          # geometry

import geopandas as gpd
from geopandas.array import GeometryDtype

# dtype 信息
dtype = GeometryDtype()
print("名称:", dtype.name)           # geometry
print("类型:", dtype.type)           # <class 'shapely.geometry.base.BaseGeometry'>
print("numpy 类型:", dtype.numpy_dtype)  # object

# 在 DataFrame 中识别
import pandas as pd
from shapely.geometry import Point

gdf = gpd.GeoDataFrame({'geometry': [Point(0, 0)]})
print(gdf.dtypes)
# geometry    geometry
# dtype: object

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

# 创建 GeoDataFrame
gdf = gpd.GeoDataFrame({
    '名称': ['A', 'B'],
    'geometry': [
        Polygon([(0,0), (2,0), (2,2), (0,2)]),
        Polygon([(1,1), (3,1), (3,3), (1,3)])
    ]
})

# 在 GeoDataFrame 上调用 - 操作活跃几何列
print("GeoDataFrame.area:")
print(gdf.area)

# 在 GeoSeries 上调用 - 直接操作
gs = gdf.geometry
print("\nGeoSeries.area:")
print(gs.area)

# 两者结果相同
print("\n结果相同:", (gdf.area == gs.area).all())  # True

# 缓冲区 - 返回 GeoSeries
buffered_gdf = gdf.buffer(0.5)
buffered_gs = gs.buffer(0.5)
print("\nbuffer 返回类型 (GeoDataFrame):", type(buffered_gdf))  # GeoSeries
print("buffer 返回类型 (GeoSeries):", type(buffered_gs))       # GeoSeries

from geopandas.array import GeometryDtype
import pandas as pd

# GeometryDtype 继承自 pandas ExtensionDtype
print(issubclass(GeometryDtype, pd.api.extensions.ExtensionDtype))  # True

# 注册为 "geometry"
dtype = GeometryDtype()
print(dtype.name)  # "geometry"

# 可以通过字符串识别
print(pd.api.types.pandas_dtype("geometry"))  # geometry

import geopandas as gpd
from shapely.geometry import Point

gdf = gpd.GeoDataFrame({
    '名称': ['A', 'B'],
    'x': [1.0, 2.0],
    'geometry': [Point(0, 0), Point(1, 1)]
})

# dtype 识别
print(gdf.dtypes)
# 名称          object
# x           float64
# geometry    geometry    ← GeometryDtype
# dtype: object

# 通过 dtype 识别几何列
for col in gdf.columns:
    if gdf[col].dtype.name == 'geometry':
        print(f"'{col}' 是几何列")

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

# GeoSeries 的 dtype 始终是 geometry
gs = gpd.GeoSeries([Point(0, 0), Point(1, 1)])
print(gs.dtype)  # geometry

# 转换为普通 Series 会丢失 geometry dtype
s = pd.Series(gs.values)
print(s.dtype)  # object

# 但仍然包含 Shapely 对象
print(type(s.iloc[0]))  # <class 'shapely.geometry.point.Point'>

import geopandas as gpd
from shapely.geometry import Point, box
import numpy as np

# 创建大量随机点
np.random.seed(42)
n = 10000
points = gpd.GeoSeries(
    gpd.points_from_xy(np.random.rand(n) * 100, np.random.rand(n) * 100)
)
gdf = gpd.GeoDataFrame({'id': range(n), 'geometry': points})

# 空间索引自动创建（延迟创建）
sindex = gdf.sindex
print("索引类型:", type(sindex))
print("索引大小:", sindex.size)

# 查询范围内的点
query_box = box(10, 10, 20, 20)  # 查询范围
candidates_idx = sindex.query(query_box)
print(f"\n范围查询结果: 在 {n} 个点中找到 {len(candidates_idx)} 个候选")

# 精确过滤
result = gdf.iloc[candidates_idx]
result = result[result.within(query_box)]
print(f"精确过滤后: {len(result)} 个点在范围内")

import geopandas as gpd
from shapely.geometry import Point, box
import numpy as np

np.random.seed(42)
points = gpd.GeoSeries(gpd.points_from_xy(np.random.rand(100), np.random.rand(100)))

sindex = points.sindex

# 1. 单个几何对象查询
query_geom = box(0.2, 0.2, 0.4, 0.4)
idx = sindex.query(query_geom)
print(f"方法1 - query(): 找到 {len(idx)} 个候选")

# 2. 批量查询（向量化，高性能）
query_geoms = gpd.GeoSeries([
    box(0, 0, 0.3, 0.3),
    box(0.5, 0.5, 0.8, 0.8)
])
tree_idx, input_idx = sindex.query(query_geoms, predicate='intersects')
print(f"方法2 - 批量 query(): 找到 {len(tree_idx)} 个结果对")

# 3. 最近邻查询
nearest_idx = sindex.nearest(Point(0.5, 0.5))
print(f"方法3 - nearest(): 最近点索引 = {nearest_idx}")

import geopandas as gpd
from shapely.geometry import box
import numpy as np
import time

# 创建测试数据
np.random.seed(42)
n = 100000
gdf = gpd.GeoDataFrame({
    'geometry': gpd.points_from_xy(np.random.rand(n), np.random.rand(n))
})
query = box(0.3, 0.3, 0.5, 0.5)

# 不使用空间索引（暴力搜索）
start = time.time()
result1 = gdf[gdf.within(query)]
time1 = time.time() - start

# 使用空间索引（两步过滤）
start = time.time()
candidates = gdf.iloc[gdf.sindex.query(query)]
result2 = candidates[candidates.within(query)]
time2 = time.time() - start

print(f"数据量: {n} 条")
print(f"暴力搜索: {time1:.4f} 秒, 结果: {len(result1)} 条")
print(f"空间索引: {time2:.4f} 秒, 结果: {len(result2)} 条")
print(f"加速比: {time1/time2:.1f}x")

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

# 继承关系验证
gdf = gpd.GeoDataFrame({'geometry': [Point(0, 0)]})
gs = gdf.geometry

# GeoDataFrame 的继承链
print("GeoDataFrame MRO:")
for cls in type(gdf).__mro__[:5]:
    print(f"  → {cls.__name__}")
# GeoDataFrame → DataFrame → NDFrame → PandasObject → object

# GeoSeries 的继承链
print("\nGeoSeries MRO:")
for cls in type(gs).__mro__[:5]:
    print(f"  → {cls.__name__}")
# GeoSeries → Series → NDFrame → PandasObject → object

# 类型检查
print("\n类型检查:")
print(f"isinstance(gdf, pd.DataFrame): {isinstance(gdf, pd.DataFrame)}")  # True
print(f"isinstance(gs, pd.Series): {isinstance(gs, pd.Series)}")          # True
print(f"isinstance(gdf, gpd.GeoDataFrame): {isinstance(gdf, gpd.GeoDataFrame)}")  # True

import geopandas as gpd

# 从文件读取
gdf = gpd.read_file("input.shp")

# 从数据库读取
gdf = gpd.read_postgis("SELECT * FROM table", engine)

# 从 Parquet 读取
gdf = gpd.read_parquet("data.parquet")

# 从 pandas DataFrame 转换
import pandas as pd
from shapely.geometry import Point

df = pd.read_csv("data.csv")
gdf = gpd.GeoDataFrame(
    df,
    geometry=gpd.points_from_xy(df.longitude, df.latitude),
    crs="EPSG:4326"
)

# 1. 检查和设置 CRS
if gdf.crs is None:
    gdf = gdf.set_crs("EPSG:4326")

# 2. 投影转换（根据需要）
gdf = gdf.to_crs(epsg=32650)  # 转为 UTM 50N

# 3. 检查和修复无效几何
invalid = gdf[~gdf.is_valid]
print(f"无效几何数量: {len(invalid)}")
gdf['geometry'] = gdf.make_valid()

# 4. 删除空几何
gdf = gdf[~gdf.is_empty]

# 5. 删除重复几何
gdf = gdf.drop_duplicates(subset='geometry')

# 6. 处理缺失值
gdf = gdf.dropna(subset=['geometry'])

# 缓冲区分析
buffered = gdf.buffer(1000)

# 空间连接
joined = gpd.sjoin(points_gdf, polygon_gdf, predicate='within')

# 空间叠加
intersection = gpd.overlay(gdf1, gdf2, how='intersection')

# 融合
dissolved = gdf.dissolve(by='category', aggfunc='sum')

# 裁剪
clipped = gpd.clip(gdf, boundary)

# 静态地图
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(12, 8))
gdf.plot(ax=ax, column='value', legend=True, cmap='YlOrRd')
plt.savefig("output_map.png", dpi=300)

# 交互式地图
m = gdf.explore(column='value', cmap='YlOrRd')
m.save("interactive_map.html")

# 导出数据
gdf.to_file("result.gpkg", driver="GPKG")
gdf.to_parquet("result.parquet")
gdf.to_postgis("result_table", engine)

import geopandas as gpd
import matplotlib.pyplot as plt

def analyze_poi_coverage(poi_file, boundary_file, buffer_distance=500):
    """分析 POI（兴趣点）的服务覆盖范围"""

    # 1. 读取数据
    pois = gpd.read_file(poi_file)
    boundaries = gpd.read_file(boundary_file)

    # 2. 数据清洗
    pois = pois.set_crs("EPSG:4326", allow_override=True)
    boundaries = boundaries.set_crs("EPSG:4326", allow_override=True)

    # 转为投影坐标系（以便使用米为单位）
    pois = pois.to_crs(epsg=32650)
    boundaries = boundaries.to_crs(epsg=32650)

    # 修复无效几何
    boundaries['geometry'] = boundaries.make_valid()

    # 3. 空间分析
    # 创建服务范围（缓冲区）
    pois['service_area'] = pois.buffer(buffer_distance)

    # 计算每个区域内的 POI 数量
    poi_counts = gpd.sjoin(pois, boundaries, predicate='within')
    poi_summary = poi_counts.groupby('区域名称').size().reset_index(name='POI数量')

    # 合并统计结果
    boundaries = boundaries.merge(poi_summary, on='区域名称', how='left')
    boundaries['POI数量'] = boundaries['POI数量'].fillna(0)

    # 4. 可视化
    fig, axes = plt.subplots(1, 2, figsize=(16, 8))

    boundaries.plot(ax=axes[0], column='POI数量', legend=True, cmap='YlGn')
    pois.plot(ax=axes[0], color='red', markersize=5)
    axes[0].set_title('POI 分布与数量')

    # 服务覆盖率
    service_union = pois.set_geometry('service_area').union_all()
    boundaries.plot(ax=axes[1], color='lightgray', edgecolor='black')
    gpd.GeoSeries([service_union]).plot(ax=axes[1], alpha=0.5, color='blue')
    axes[1].set_title(f'POI 服务覆盖范围（{buffer_distance}m）')

    plt.tight_layout()
    plt.savefig("poi_analysis.png", dpi=150)

    # 5. 输出
    boundaries.to_file("poi_analysis_result.gpkg", driver="GPKG")

    return boundaries