import geopandas as gpd

# 1. 在 GeoPandas 中加载或创建数据
gdf = gpd.read_file("data.shp")

# 2. 提取所需信息传递给第三方库
geometries = gdf.geometry
attributes = gdf[["col1", "col2"]]

# 3. 调用第三方库进行处理
# result = third_party_lib.process(geometries, attributes)

# 4. 将结果整合回 GeoDataFrame
# gdf["new_col"] = result

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

# 创建 Shapely 几何对象
point = Point(116.4, 39.9)
polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])

# 直接用于构建 GeoDataFrame
gdf = gpd.GeoDataFrame(
    {"name": ["北京", "矩形区域"]},
    geometry=[point, polygon],
    crs="EPSG:4326"
)
print(type(gdf.geometry.iloc[0]))
# <class 'shapely.geometry.point.Point'>

from shapely.geometry import MultiPoint
from shapely.ops import nearest_points, split, snap
from shapely import wkt

# 最近点计算
multi = MultiPoint([(0, 0), (1, 1), (2, 2)])
target = Point(0.8, 0.8)
nearest = nearest_points(multi, target)
print(f"最近点: {nearest[0]}")

# WKT 互转
wkt_str = polygon.wkt
restored = wkt.loads(wkt_str)
print(f"WKT: {wkt_str}")

# 几何分割
line = LineString([(0, 0), (2, 2)])
splitter = LineString([(1, 0), (1, 2)])
result = split(line, splitter)
print(f"分割后片段数: {len(result.geoms)}")

import numpy as np

# 使用 apply 对每个几何对象执行 Shapely 操作
gdf = gpd.GeoDataFrame(
    geometry=[Point(x, y) for x, y in np.random.rand(100, 2) * 10],
    crs="EPSG:4326"
)

# 为每个点创建缓冲区并简化
gdf["buffered"] = gdf.geometry.buffer(0.5)
gdf["simplified"] = gdf["buffered"].simplify(tolerance=0.1)

# 提取坐标
gdf["x"] = gdf.geometry.x
gdf["y"] = gdf.geometry.y
print(gdf[["x", "y"]].head())

import fiona

# 直接使用 Fiona 查看数据源信息
with fiona.open("data.shp") as src:
    print(f"驱动: {src.driver}")
    print(f"CRS: {src.crs}")
    print(f"要素数量: {len(src)}")
    print(f"字段: {src.schema['properties']}")
    print(f"几何类型: {src.schema['geometry']}")

# 指定 pyogrio 作为引擎
gdf = gpd.read_file("large_dataset.gpkg", engine="pyogrio")

# pyogrio 支持读取部分列，减少内存占用
gdf_partial = gpd.read_file(
    "large_dataset.gpkg",
    engine="pyogrio",
    columns=["name", "population"]
)

# 写出时也可指定引擎
gdf.to_file("output.gpkg", driver="GPKG", engine="pyogrio")

# 性能对比示例
import time

for engine in ["fiona", "pyogrio"]:
    start = time.time()
    gdf = gpd.read_file("large_file.gpkg", engine=engine)
    elapsed = time.time() - start
    print(f"{engine}: {elapsed:.2f} 秒, {len(gdf)} 条记录")

import rasterio
from rasterio.plot import show

# 读取栅格数据基本信息
with rasterio.open("elevation.tif") as src:
    print(f"尺寸: {src.width} x {src.height}")
    print(f"波段数: {src.count}")
    print(f"CRS: {src.crs}")
    print(f"范围: {src.bounds}")
    data = src.read(1)  # 读取第一波段

import rasterio
from rasterio.mask import mask
import json

# 读取矢量边界
boundary = gpd.read_file("study_area.shp")

# 确保坐标系一致
with rasterio.open("elevation.tif") as src:
    boundary_reproj = boundary.to_crs(src.crs)
    shapes = [json.loads(boundary_reproj.to_json())["features"][0]["geometry"]]

    # 裁剪栅格
    out_image, out_transform = mask(src, shapes, crop=True)
    out_meta = src.meta.copy()
    out_meta.update({
        "height": out_image.shape[1],
        "width": out_image.shape[2],
        "transform": out_transform
    })

# 写出裁剪后的栅格
with rasterio.open("clipped.tif", "w", **out_meta) as dest:
    dest.write(out_image)

from rasterstats import zonal_stats

# 对每个多边形区域计算栅格统计量
zones = gpd.read_file("districts.shp")
stats = zonal_stats(
    zones,
    "elevation.tif",
    stats=["min", "max", "mean", "median", "std", "count"]
)

# 将统计结果合并到 GeoDataFrame
import pandas as pd
stats_df = pd.DataFrame(stats)
zones = pd.concat([zones, stats_df], axis=1)
print(zones[["name", "min", "max", "mean"]].head())

import rasterio
import numpy as np

points = gpd.read_file("sample_points.shp")

with rasterio.open("elevation.tif") as src:
    points_reproj = points.to_crs(src.crs)
    coords = [(geom.x, geom.y) for geom in points_reproj.geometry]
    values = [val[0] for val in src.sample(coords)]

points["elevation"] = values
print(points[["id", "elevation"]].head())

import geopandas as gpd
import numpy as np
from sklearn.cluster import KMeans, DBSCAN
from sklearn.preprocessing import StandardScaler

# 加载空间数据
gdf = gpd.read_file("buildings.shp")

# 提取特征：坐标 + 属性
coords = np.column_stack([gdf.geometry.x, gdf.geometry.y])
features = gdf[["area", "height", "year_built"]].values
X = np.hstack([coords, features])

# 标准化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# K-Means 聚类
kmeans = KMeans(n_clusters=5, random_state=42, n_init=10)
gdf["kmeans_cluster"] = kmeans.fit_predict(X_scaled)

# DBSCAN 密度聚类（适合发现不规则形状的簇）
dbscan = DBSCAN(eps=0.5, min_samples=10)
gdf["dbscan_cluster"] = dbscan.fit_predict(X_scaled)

# 可视化聚类结果
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
gdf.plot(column="kmeans_cluster", cmap="Set2", ax=axes[0], legend=True)
axes[0].set_title("K-Means 聚类")
gdf.plot(column="dbscan_cluster", cmap="Set2", ax=axes[1], legend=True)
axes[1].set_title("DBSCAN 聚类")
plt.tight_layout()
plt.show()

from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report

# 准备分类数据
gdf = gpd.read_file("land_use.shp")
feature_cols = ["ndvi", "slope", "elevation", "distance_to_road"]
X = gdf[feature_cols].values
y = gdf["land_type"].values

# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=42
)

# 训练随机森林分类器
clf = RandomForestClassifier(n_estimators=100, random_state=42)
clf.fit(X_train, y_train)

# 评估
y_pred = clf.predict(X_test)
print(classification_report(y_test, y_pred))

# 特征重要性
importance = dict(zip(feature_cols, clf.feature_importances_))
for feat, imp in sorted(importance.items(), key=lambda x: -x[1]):
    print(f"  {feat}: {imp:.4f}")

import networkx as nx
import geopandas as gpd

# 读取道路网络数据
roads = gpd.read_file("roads.shp")

# 构建无向图
G = nx.Graph()
for idx, row in roads.iterrows():
    coords = list(row.geometry.coords)
    start = coords[0]
    end = coords[-1]
    length = row.geometry.length
    G.add_edge(
        start, end,
        weight=length,
        name=row.get("road_name", ""),
        fid=idx
    )

print(f"节点数: {G.number_of_nodes()}")
print(f"边数: {G.number_of_edges()}")

# 寻找最近的图节点
from shapely.geometry import Point

def find_nearest_node(G, point):
    """找到距给定点最近的图节点"""
    min_dist = float("inf")
    nearest = None
    for node in G.nodes():
        dist = Point(node).distance(point)
        if dist < min_dist:
            min_dist = dist
            nearest = node
    return nearest

# 起终点
origin = Point(116.35, 39.92)
destination = Point(116.45, 39.98)
start_node = find_nearest_node(G, origin)
end_node = find_nearest_node(G, destination)

# 计算最短路径
path = nx.shortest_path(G, start_node, end_node, weight="weight")
path_length = nx.shortest_path_length(G, start_node, end_node, weight="weight")
print(f"最短路径长度: {path_length:.2f}")
print(f"途经节点数: {len(path)}")

# 将路径转换为 GeoDataFrame 以便可视化
from shapely.geometry import LineString
path_line = LineString(path)
path_gdf = gpd.GeoDataFrame(
    {"type": ["最短路径"]},
    geometry=[path_line],
    crs=roads.crs
)

# 计算各种中心性指标
degree_cent = nx.degree_centrality(G)
betweenness_cent = nx.betweenness_centrality(G, weight="weight")

# 将结果映射到节点 GeoDataFrame
nodes_data = []
for node in G.nodes():
    nodes_data.append({
        "geometry": Point(node),
        "degree": degree_cent[node],
        "betweenness": betweenness_cent[node]
    })

nodes_gdf = gpd.GeoDataFrame(nodes_data, crs=roads.crs)
print(nodes_gdf.nlargest(5, "betweenness")[["betweenness", "degree"]])

import osmnx as ox

# 按城市名下载街道网络
G = ox.graph_from_place("朝阳区, 北京, 中国", network_type="drive")
print(f"节点数: {G.number_of_nodes()}")
print(f"边数: {G.number_of_edges()}")

# 将图转换为 GeoDataFrame
nodes_gdf, edges_gdf = ox.graph_to_gdfs(G)
print(f"节点列: {list(nodes_gdf.columns[:5])}")
print(f"边列: {list(edges_gdf.columns[:5])}")

# 下载建筑物轮廓
buildings = ox.features_from_place(
    "朝阳区, 北京, 中国",
    tags={"building": True}
)
print(f"建筑物数量: {len(buildings)}")

# 下载特定类型的 POI
restaurants = ox.features_from_place(
    "朝阳区, 北京, 中国",
    tags={"amenity": "restaurant"}
)
print(f"餐厅数量: {len(restaurants)}")

# 结果已经是 GeoDataFrame，可直接使用 GeoPandas 分析
print(restaurants[["name", "geometry"]].head())

# 计算最短路径
origin_node = ox.nearest_nodes(G, 116.43, 39.92)
dest_node = ox.nearest_nodes(G, 116.48, 39.95)
route = ox.shortest_path(G, origin_node, dest_node, weight="length")

# 计算路径长度
route_length = sum(
    ox.utils_graph.get_route_edge_attributes(G, route, "length")
)
print(f"路径长度: {route_length:.0f} 米")

# 计算基本路网统计
stats = ox.basic_stats(G)
print(f"路网总长度: {stats['edge_length_total']:.0f} 米")
print(f"平均街道长度: {stats['edge_length_avg']:.0f} 米")

import movingpandas as mpd
import pandas as pd
from shapely.geometry import Point
from datetime import datetime, timedelta

# 从 GeoDataFrame 创建轨迹
data = {
    "id": [1]*5 + [2]*5,
    "timestamp": [
        datetime(2024, 1, 1, 8, i*10) for i in range(5)
    ] + [
        datetime(2024, 1, 1, 8, i*10) for i in range(5)
    ],
    "geometry": [
        Point(116.3 + i*0.01, 39.9 + i*0.005) for i in range(5)
    ] + [
        Point(116.4 + i*0.01, 39.85 + i*0.005) for i in range(5)
    ]
}

gdf = gpd.GeoDataFrame(data, crs="EPSG:4326")
gdf = gdf.set_index("timestamp")

# 创建轨迹集合
traj_collection = mpd.TrajectoryCollection(gdf, traj_id_col="id")
print(f"轨迹数量: {len(traj_collection)}")
for traj in traj_collection:
    print(f"  轨迹 {traj.id}: {traj.get_length():.2f} 米")

# 轨迹分割——按时间间隔分割
splitter = mpd.ObservationGapSplitter(traj_collection)
split_trajs = splitter.split(gap=timedelta(minutes=30))
print(f"分割后轨迹数: {len(split_trajs)}")

# 轨迹速度计算
for traj in traj_collection:
    traj.add_speed(overwrite=True)
    traj_gdf = traj.to_point_gdf()
    print(f"轨迹 {traj.id} 平均速度: {traj_gdf['speed'].mean():.2f} m/s")

# 轨迹简化
for traj in traj_collection:
    generalized = mpd.DouglasPeuckerGeneralizer(traj)
    simplified = generalized.generalize(tolerance=0.001)
    print(f"简化前: {traj.size()} 点 -> 简化后: {simplified.size()} 点")

# 裁剪轨迹到指定区域
from shapely.geometry import box
study_area = box(116.3, 39.85, 116.5, 39.95)

for traj in traj_collection:
    clipped = traj.clip(study_area)
    if clipped is not None and clipped.size() > 0:
        print(f"轨迹 {traj.id} 在研究区内的点数: {clipped.size()}")

# 将轨迹转回 GeoDataFrame
all_lines = []
for traj in traj_collection:
    line_gdf = traj.to_line_gdf()
    all_lines.append(line_gdf)

result_gdf = pd.concat(all_lines)
print(f"轨迹线段总数: {len(result_gdf)}")

import geopandas as gpd
from libpysal.weights import Queen, KNN, DistanceBand

# 加载面状数据
gdf = gpd.read_file("districts.shp")

# 基于邻接关系的空间权重（Queen 方式：共享边或顶点）
w_queen = Queen.from_dataframe(gdf)
print(f"Queen 权重 - 平均邻居数: {w_queen.mean_neighbors:.2f}")

# 基于 K 近邻的空间权重
w_knn = KNN.from_dataframe(gdf, k=5)
print(f"KNN(k=5) 权重 - 每个要素恰好有 {w_knn.max_neighbors} 个邻居")

# 基于距离的空间权重
w_dist = DistanceBand.from_dataframe(gdf, threshold=5000)
print(f"距离权重 - 平均邻居数: {w_dist.mean_neighbors:.2f}")

from esda.moran import Moran, Moran_Local
import numpy as np

# 全局 Moran's I
y = gdf["income"].values
moran = Moran(y, w_queen)
print(f"Moran's I: {moran.I:.4f}")
print(f"p 值: {moran.p_sim:.4f}")
print(f"z 值: {moran.z_sim:.4f}")

if moran.p_sim < 0.05:
    if moran.I > 0:
        print("结论：存在显著的正空间自相关（聚集模式）")
    else:
        print("结论：存在显著的负空间自相关（离散模式）")
else:
    print("结论：空间分布接近随机")

# 局部 Moran's I
lisa = Moran_Local(y, w_queen)

# 将结果加入 GeoDataFrame
gdf["lisa_I"] = lisa.Is
gdf["lisa_q"] = lisa.q  # 象限：1=HH, 2=LH, 3=LL, 4=HL
gdf["lisa_p"] = lisa.p_sim

# 标记显著的聚类
significance = 0.05
gdf["lisa_cluster"] = "不显著"
sig_mask = gdf["lisa_p"] < significance

cluster_labels = {1: "高-高聚类", 2: "低-高离群", 3: "低-低聚类", 4: "高-低离群"}
for q, label in cluster_labels.items():
    mask = sig_mask & (gdf["lisa_q"] == q)
    gdf.loc[mask, "lisa_cluster"] = label

print(gdf["lisa_cluster"].value_counts())

from spreg import OLS, ML_Lag, ML_Error
import numpy as np

# 准备数据
y = gdf[["income"]].values
X = gdf[["education", "employment", "accessibility"]].values

# OLS 回归（对照组）
ols = OLS(y, X, w=w_queen, name_y="income",
          name_x=["education", "employment", "accessibility"])
print(ols.summary)

# 空间滞后模型
lag = ML_Lag(y, X, w=w_queen, name_y="income",
             name_x=["education", "employment", "accessibility"])
print(f"空间滞后系数 (rho): {lag.rho:.4f}")

# 空间误差模型
error = ML_Error(y, X, w=w_queen, name_y="income",
                 name_x=["education", "employment", "accessibility"])
print(f"空间误差系数 (lambda): {error.lam:.4f}")

import plotly.express as px
import geopandas as gpd
import json

# 加载数据
gdf = gpd.read_file("provinces.shp")

# 使用 Plotly Express 创建分级统计图
fig = px.choropleth(
    gdf,
    geojson=json.loads(gdf.to_json()),
    locations=gdf.index,
    color="population",
    color_continuous_scale="YlOrRd",
    hover_name="name",
    hover_data=["gdp", "area"],
    title="中国各省人口分布"
)
fig.update_geos(fitbounds="locations", visible=False)
fig.show()

# 点数据的交互式地图
points = gpd.read_file("cities.shp")

fig = px.scatter_geo(
    points,
    lat=points.geometry.y,
    lon=points.geometry.x,
    size="population",
    color="region",
    hover_name="city_name",
    title="中国主要城市分布",
    size_max=30
)
fig.update_layout(geo=dict(
    scope="asia",
    center=dict(lat=35, lon=105),
    projection_scale=3
))
fig.show()

import plotly.graph_objects as go
import numpy as np

# 三维地形可视化
gdf = gpd.read_file("terrain_points.shp")

fig = go.Figure(data=[go.Scatter3d(
    x=gdf.geometry.x,
    y=gdf.geometry.y,
    z=gdf["elevation"],
    mode="markers",
    marker=dict(
        size=2,
        color=gdf["elevation"],
        colorscale="Earth",
        colorbar=dict(title="海拔 (m)"),
        opacity=0.8
    ),
    text=gdf["name"],
    hovertemplate="经度: %{x:.2f}<br>纬度: %{y:.2f}<br>海拔: %{z:.0f}m"
)])

fig.update_layout(
    title="地形三维可视化",
    scene=dict(
        xaxis_title="经度",
        yaxis_title="纬度",
        zaxis_title="海拔 (m)",
        aspectmode="manual",
        aspectratio=dict(x=2, y=2, z=0.5)
    ),
    width=900,
    height=700
)
fig.show()

# 时序变化的动态可视化
gdf = gpd.read_file("yearly_data.shp")

fig = px.choropleth(
    gdf,
    geojson=json.loads(gdf.to_json()),
    locations=gdf.index,
    color="gdp",
    animation_frame="year",
    color_continuous_scale="Viridis",
    hover_name="province",
    title="各省 GDP 年度变化"
)
fig.update_geos(fitbounds="locations", visible=False)
fig.update_layout(
    updatemenus=[dict(type="buttons", showactive=False,
                      buttons=[dict(label="播放",
                                    method="animate",
                                    args=[None, {"frame": {"duration": 500}}])])]
)
fig.show()