import geopandas as gpd

# 读取 Shapefile
gdf = gpd.read_file("cities.shp")

# 读取 GeoJSON
gdf = gpd.read_file("boundaries.geojson")

# 读取 GeoPackage
gdf = gpd.read_file("data.gpkg", layer="buildings")

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

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

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

# 距离计算
distances = gdf1.distance(gdf2)

# 面积计算
areas = gdf.area

# 质心计算
centroids = gdf.centroid

# 查看 CRS
print(gdf.crs)  # EPSG:4326

# 投影转换
gdf_projected = gdf.to_crs(epsg=3857)

# 自定义 CRS
from pyproj import CRS
custom_crs = CRS.from_proj4("+proj=aea +lat_1=25 +lat_2=47 +lon_0=105")
gdf_custom = gdf.to_crs(custom_crs)

import matplotlib.pyplot as plt

# 基础绘图
gdf.plot(column='population', cmap='YlOrRd', legend=True)
plt.title("人口分布图")
plt.show()

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

# 空间索引自动创建
sindex = gdf.sindex

# 查询某个范围内的要素
bbox = (116.0, 39.0, 117.0, 40.0)
candidates = list(sindex.query(bbox))

import geopandas as gpd

# 读取城市用地规划数据
land_use = gpd.read_file("land_use.shp")

# 计算各类用地面积
land_use['面积_平方公里'] = land_use.to_crs(epsg=32650).area / 1e6

# 按用地类型统计
summary = land_use.groupby('用地类型')['面积_平方公里'].sum()
print(summary)

# 读取监测站点和行政区划数据
stations = gpd.read_file("monitoring_stations.geojson")
districts = gpd.read_file("districts.geojson")

# 空间连接：将站点关联到所在区域
joined = gpd.sjoin(stations, districts, predicate='within')

# 按区域计算平均污染指数
avg_pollution = joined.groupby('区域名称')['PM25'].mean()

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

# 计算道路长度
roads['长度_km'] = roads.to_crs(epsg=32650).length / 1000

# 按道路等级统计里程
mileage = roads.groupby('道路等级')['长度_km'].sum()

# 创建道路缓冲区（模拟噪音影响范围）
noise_zones = roads.to_crs(epsg=32650).buffer(50)  # 50米缓冲区

# 读取房产和配套设施数据
properties = gpd.read_file("properties.geojson")
schools = gpd.read_file("schools.geojson")

# 计算每个房产到最近学校的距离
properties_proj = properties.to_crs(epsg=32650)
schools_proj = schools.to_crs(epsg=32650)

# 计算最近学校距离
from shapely.ops import nearest_points
properties_proj['到最近学校距离'] = properties_proj.geometry.apply(
    lambda x: schools_proj.distance(x).min()
)

import geopandas as gpd
import matplotlib.pyplot as plt

# ========================================
# 步骤 1：读取数据
# ========================================
# 使用 GeoPandas 内置的世界各国数据集
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

# 查看数据概览
print("数据形状:", world.shape)
print("\n列名:", world.columns.tolist())
print("\n前5行数据:")
print(world.head())
print("\n坐标参考系:", world.crs)

# ========================================
# 步骤 2：数据探索
# ========================================
# 查看数据类型
print("\n数据类型:")
print(world.dtypes)

# 基本统计
print("\n人口统计:")
print(world['pop_est'].describe())

# 按大洲分组统计
continent_stats = world.groupby('continent').agg({
    'pop_est': ['sum', 'mean', 'count'],
    'gdp_md_est': 'sum'
}).round(0)
print("\n大洲统计:")
print(continent_stats)

# ========================================
# 步骤 3：空间分析
# ========================================
# 计算各国面积（使用等面积投影）
world_aea = world.to_crs('+proj=aea +lat_1=20 +lat_2=60 +lon_0=0')
world['area_km2'] = world_aea.area / 1e6  # 转换为平方公里

# 计算人口密度
world['pop_density'] = world['pop_est'] / world['area_km2']

# 找出人口密度最高的 10 个国家
top10_density = world.nlargest(10, 'pop_density')[['name', 'pop_density', 'pop_est']]
print("\n人口密度 TOP 10:")
print(top10_density)

# ========================================
# 步骤 4：空间操作
# ========================================
# 获取各国质心
world['centroid'] = world.centroid

# 计算中国的邻国
china = world[world['name'] == 'China']
neighbors = world[world.touches(china.unary_union)]
print(f"\n中国的邻国 ({len(neighbors)} 个):")
print(neighbors['name'].tolist())

# ========================================
# 步骤 5：可视化
# ========================================
fig, axes = plt.subplots(2, 2, figsize=(20, 16))

# 子图1：世界政区图
world.plot(ax=axes[0, 0], color='lightblue', edgecolor='gray')
axes[0, 0].set_title('世界政区图', fontsize=14)
axes[0, 0].set_axis_off()

# 子图2：人口分布图
world.plot(column='pop_est', ax=axes[0, 1], legend=True,
           legend_kwds={'label': '人口数量', 'shrink': 0.6},
           cmap='YlOrRd', missing_kwds={'color': 'lightgrey'})
axes[0, 1].set_title('世界人口分布', fontsize=14)
axes[0, 1].set_axis_off()

# 子图3：人均GDP
world.plot(column='gdp_md_est', ax=axes[1, 0], legend=True,
           legend_kwds={'label': 'GDP (百万美元)', 'shrink': 0.6},
           cmap='Greens', missing_kwds={'color': 'lightgrey'})
axes[1, 0].set_title('世界 GDP 分布', fontsize=14)
axes[1, 0].set_axis_off()

# 子图4：人口密度
world.plot(column='pop_density', ax=axes[1, 1], legend=True,
           legend_kwds={'label': '人口密度 (人/km²)', 'shrink': 0.6},
           cmap='OrRd', scheme='quantiles',
           missing_kwds={'color': 'lightgrey'})
axes[1, 1].set_title('世界人口密度分布', fontsize=14)
axes[1, 1].set_axis_off()

plt.tight_layout()
plt.savefig('world_analysis.png', dpi=150, bbox_inches='tight')
plt.show()

# ========================================
# 步骤 6：数据导出
# ========================================
# 导出为 GeoJSON
world.to_file("world_with_density.geojson", driver="GeoJSON")

# 导出为 GeoPackage
world.to_file("world_analysis.gpkg", layer="countries", driver="GPKG")

print("\n分析完成！数据已导出。")