geopandas-geospatial

GeoPandas 地理空间分析

概述

GeoPandas 在 pandas 的基础上扩展了几何类型的空间操作，将 pandas DataFrame 与 Shapely 几何对象以及用于文件 I/O 的 Fiona 结合起来。它支持以熟悉的 pandas 风格 API 读取、写入、处理和可视化地理空间矢量数据（点、线和多边形）。

适用场景

• 读取和写入空间文件格式（Shapefile、GeoJSON、GeoPackage、Parquet）
• 对地理数据集执行空间连接（多边形内的点、最近邻）
• 执行叠加操作（交集、并集、差集、裁剪）
• 计算几何属性（面积、距离、缓冲区、质心）
• 创建分级统计地图和交互式 Web 地图
• 在不同坐标参考系之间重投影数据
• 按属性聚合空间要素（融合）
• 对于栅格数据分析，请改用 rasterio/xarray
• 对于大规模分布式地理空间处理，可考虑 Dask-GeoPandas 或 Apache Sedona

前置条件

pip install geopandas matplotlib
# Optional:
# pip install folium       — interactive maps
# pip install mapclassify  — classification schemes for choropleth
# pip install contextily   — basemaps
# pip install pyarrow      — faster I/O (2-4x speedup)
# pip install psycopg2 geoalchemy2  — PostGIS support

快速开始

import geopandas as gpd

# Read spatial data
gdf = gpd.read_file("data.geojson")
print(f"Shape: {gdf.shape}, CRS: {gdf.crs}")
print(f"Geometry types: {gdf.geometry.geom_type.unique()}")

# Reproject, compute area, save
gdf_proj = gdf.to_crs("EPSG:3857")
gdf_proj['area_m2'] = gdf_proj.geometry.area
gdf_proj.to_file("output.gpkg")

# Quick map
gdf.plot(column='population', legend=True, figsize=(10, 8))

核心 API

1. 数据 I/O

import geopandas as gpd

# Read various formats
gdf = gpd.read_file("data.shp")           # Shapefile
gdf = gpd.read_file("data.geojson")       # GeoJSON
gdf = gpd.read_file("data.gpkg")          # GeoPackage
gdf = gpd.read_file("data.gpkg", layer="roads")  # Specific layer

# Filtered reading (load only needed data)
gdf = gpd.read_file("data.gpkg", bbox=(xmin, ymin, xmax, ymax))
gdf = gpd.read_file("data.gpkg", columns=["name", "geometry"])
gdf = gpd.read_file("data.gpkg", where="population > 10000")

# Arrow acceleration (2-4x faster)
gdf = gpd.read_file("data.gpkg", use_arrow=True)

# Parquet/Feather (columnar, fast, preserves CRS)
gdf = gpd.read_parquet("data.parquet")
gdf.to_parquet("output.parquet")

# PostGIS database
from sqlalchemy import create_engine
engine = create_engine("postgresql://user:pass@host/db")
gdf = gpd.read_postgis("SELECT * FROM parcels", con=engine, geom_col='geom')
gdf.to_postgis("output_table", con=engine)

# Write
gdf.to_file("output.gpkg")             # GeoPackage (recommended)
gdf.to_file("output.shp")              # Shapefile
gdf.to_file("output.geojson", driver="GeoJSON")

2. CRS 管理

# Check current CRS
print(gdf.crs)                  # e.g., EPSG:4326
print(gdf.crs.is_geographic)    # True for lat/lon
print(gdf.crs.is_projected)     # True for meters

# Reproject (transforms coordinates)
gdf_proj = gdf.to_crs("EPSG:3857")       # Web Mercator
gdf_proj = gdf.to_crs(epsg=32633)        # UTM zone 33N

# Set CRS (only when metadata missing, does NOT transform coordinates)
gdf = gdf.set_crs("EPSG:4326")

# Estimate appropriate UTM zone
utm_crs = gdf.estimate_utm_crs()
gdf_utm = gdf.to_crs(utm_crs)

常见 EPSG 代码：

代码	名称	用途
4326	WGS 84	GPS 坐标、Web 数据
3857	Web Mercator	Web 地图（Google/OSM 瓦片）
326xx	UTM 分区（N）	面积/距离计算
5070	Albers 等积投影（美国）	保持面积的美国地图

3. 几何操作

# Buffer (expand/erode geometry by distance)
buffered = gdf.geometry.buffer(100)      # 100 units (meters if projected)
eroded = gdf.geometry.buffer(-50)        # Negative = erosion

# Simplify (reduce complexity)
simplified = gdf.geometry.simplify(tolerance=10, preserve_topology=True)

# Centroid, convex hull, envelope
centroids = gdf.geometry.centroid
hulls = gdf.geometry.convex_hull
bounds = gdf.geometry.envelope

# Union all geometries
unified = gdf.geometry.union_all()

# Affine transformations
rotated = gdf.geometry.rotate(angle=45, origin='center')
scaled = gdf.geometry.scale(xfact=2.0, yfact=2.0)
translated = gdf.geometry.translate(xoff=100, yoff=50)

# Geometric properties
areas = gdf.geometry.area           # Use projected CRS for accuracy
lengths = gdf.geometry.length       # Perimeter for polygons
is_valid = gdf.geometry.is_valid    # Validate geometry
total = gdf.geometry.total_bounds   # [minx, miny, maxx, maxy]

4. 空间分析

# Spatial join (combine datasets by spatial relationship)
joined = gpd.sjoin(points_gdf, polygons_gdf, predicate='intersects')
joined = gpd.sjoin(gdf1, gdf2, predicate='within')
joined = gpd.sjoin(gdf1, gdf2, predicate='contains', how='left')

# Nearest neighbor join
nearest = gpd.sjoin_nearest(gdf1, gdf2, max_distance=1000, distance_col='dist')

# Overlay operations (set-theoretic)
intersection = gpd.overlay(gdf1, gdf2, how='intersection')
union = gpd.overlay(gdf1, gdf2, how='union')
difference = gpd.overlay(gdf1, gdf2, how='difference')
sym_diff = gpd.overlay(gdf1, gdf2, how='symmetric_difference')

# Dissolve (aggregate by attribute)
dissolved = gdf.dissolve(by='region', aggfunc='sum')
dissolved = gdf.dissolve(by='region', aggfunc={'population': 'sum', 'area': 'mean'})

# Clip to boundary
clipped = gpd.clip(gdf, boundary_gdf)

# Distance calculations (use projected CRS)
distances = gdf.geometry.distance(single_point)

# Spatial predicates
within_mask = gdf1.geometry.within(gdf2.geometry)
intersects_mask = gdf1.geometry.intersects(gdf2.geometry)

5. 可视化

import matplotlib.pyplot as plt

# Basic plot
gdf.plot(figsize=(10, 8))

# Choropleth map
gdf.plot(column='population', cmap='YlOrRd', legend=True, figsize=(12, 8))

# Classification schemes (requires mapclassify)
gdf.plot(column='population', scheme='quantiles', k=5, legend=True)
gdf.plot(column='population', scheme='fisher_jenks', k=5, legend=True)

# Multi-layer map
fig, ax = plt.subplots(figsize=(12, 10))
polygons_gdf.plot(ax=ax, color='lightblue', edgecolor='black')
points_gdf.plot(ax=ax, color='red', markersize=10)
roads_gdf.plot(ax=ax, color='gray', linewidth=0.5)
ax.set_title('Multi-layer Map')
ax.set_axis_off()

# Interactive map (requires folium)
m = gdf.explore(column='population', cmap='YlOrRd', legend=True,
                tooltip=['name', 'population'])
m.save('map.html')

# Multi-layer interactive
m = gdf1.explore(color='blue', name='Layer 1')
gdf2.explore(m=m, color='red', name='Layer 2')
import folium
folium.LayerControl().add_to(m)

# Basemap (requires contextily)
import contextily as ctx
gdf_wm = gdf.to_crs(epsg=3857)
ax = gdf_wm.plot(alpha=0.5, figsize=(10, 10))
ctx.add_basemap(ax)

核心概念

数据结构

• GeoSeries：包含 Shapely 几何对象并提供空间方法（面积、距离、缓冲区等）的 Pandas Series
• GeoDataFrame：包含一个或多个几何列的 Pandas DataFrame。其中一列是供空间方法使用的“活动几何列”

from shapely.geometry import Point

# Create from coordinates
gdf = gpd.GeoDataFrame(
    {'name': ['A', 'B'], 'value': [10, 20]},
    geometry=[Point(0, 0), Point(1, 1)],
    crs="EPSG:4326"
)

# Multiple geometry columns
gdf['centroid'] = gdf.geometry.centroid
gdf = gdf.set_geometry('centroid')  # Switch active geometry

空间操作的 CRS 规则

• 在执行任何空间操作前，务必检查 CRS：print(gdf.crs)
• 在执行空间连接、叠加或距离计算前，确保 CRS 一致
• 使用投影 CRS（米）计算面积/距离——地理 CRS（度）会导致结果错误
• set_crs() 只添加元数据；to_crs() 会转换坐标

空间索引

GeoPandas 会为 sjoin、overlay 和其他空间操作自动创建空间索引（R-tree）。手动查询方法：

sindex = gdf.sindex
possible_idx = list(sindex.intersection((xmin, ymin, xmax, ymax)))

常见工作流

工作流 1：加载 → 转换 → 分析 → 导出

import geopandas as gpd

# Load
gdf = gpd.read_file("parcels.shp")
print(f"CRS: {gdf.crs}, Rows: {len(gdf)}")

# Transform to projected CRS for measurements
gdf = gdf.to_crs(gdf.estimate_utm_crs())

# Analyze
gdf['area_ha'] = gdf.geometry.area / 10000  # hectares
gdf['perimeter_m'] = gdf.geometry.length
print(f"Total area: {gdf['area_ha'].sum():.1f} ha")

# Export
gdf.to_file("parcels_analyzed.gpkg")

工作流 2：空间连接与聚合

# Count points per polygon
points_in_poly = gpd.sjoin(points_gdf, polygons_gdf, predicate='within')
counts = points_in_poly.groupby('index_right').agg(
    point_count=('geometry', 'size'),
    total_value=('value', 'sum')
)
result = polygons_gdf.merge(counts, left_index=True, right_index=True, how='left')
result['point_count'] = result['point_count'].fillna(0)
print(f"Polygons with points: {(result['point_count'] > 0).sum()}/{len(result)}")

工作流 3：多源数据整合

# Read from different sources, ensure matching CRS
roads = gpd.read_file("roads.shp")
buildings = gpd.read_file("buildings.geojson")
parcels = gpd.read_postgis("SELECT * FROM parcels", con=engine, geom_col='geom')

target_crs = roads.crs
buildings = buildings.to_crs(target_crs)
parcels = parcels.to_crs(target_crs)

# Find buildings within 50m of roads
buildings_near_roads = gpd.sjoin_nearest(
    buildings, roads, max_distance=50, distance_col='road_dist'
)
print(f"Buildings near roads: {len(buildings_near_roads)}/{len(buildings)}")

关键参数

参数	函数	默认值	作用
predicate	sjoin	'intersects'	空间关系：intersects、within、contains、touches、crosses
how	sjoin、overlay	'inner'	连接类型：inner、left、right
max_distance	sjoin_nearest	None	搜索半径限制（可提升性能）
k	sjoin_nearest	1	要查找的最近邻数量
tolerance	simplify	必填	Douglas-Peucker 容差（以 CRS 单位计）
preserve_topology	simplify	True	防止自相交
resolution	buffer	16	缓冲曲线的线段数
aggfunc	dissolve	'first'	非几何列的聚合函数
use_arrow	read_file	False	启用 Arrow 加速（快 2-4 倍）
scheme	plot	None	分类方式：quantiles、equal_interval、fisher_jenks
k	plot（使用 scheme 时）	5	分类区间数

最佳实践

1. 在空间操作前始终检查 CRS——CRS 不匹配会导致结果错误或为空
2. 使用投影 CRS 计算面积和距离——地理 CRS（度）不适用于测量
3. 空间连接前匹配 CRS——在调用 gpd.sjoin(gdf1, gdf2) 前运行 gdf2 = gdf2.to_crs(gdf1.crs)
4. 在复杂操作前使用 .is_valid 验证几何——无效几何会导致未显式提示的错误
5. 优先使用 GeoPackage，而不是 Shapefile——没有 10 个字符的列名限制，支持多个图层，性能更好
6. 读取时进行筛选——对于大型文件，使用 bbox、columns、where 仅加载所需数据
7. 在 sjoin_nearest 中设置 max_distance——对大型数据集进行不设上限的最近邻搜索会很慢
8. 修改几何时使用 .copy()——避免对原始 GeoDataFrame 产生意外副作用

常用配方

配方：统计多边形内的点

适用场景：按区域（县、流域）汇总点观测数据（采样点、观测值）。

import geopandas as gpd

regions = gpd.read_file("regions.geojson")
points = gpd.read_file("observations.geojson").to_crs(regions.crs)

joined = gpd.sjoin(points, regions, how="inner", predicate="within")
counts = joined.groupby("index_right").size().rename("point_count")
regions_with_counts = regions.join(counts).fillna(0)
regions_with_counts["point_count"] = regions_with_counts["point_count"].astype(int)
print(regions_with_counts[["name", "point_count"]].sort_values("point_count", ascending=False).head())

配方：为点创建缓冲区并融合重叠区域

适用场景：根据点位置创建服务区或集水区。

import geopandas as gpd

facilities = gpd.read_file("facilities.geojson")
utm_crs = facilities.estimate_utm_crs()   # Project to meters
facilities_m = facilities.to_crs(utm_crs)

# 1 km buffer
buffers = facilities_m.copy()
buffers["geometry"] = facilities_m.geometry.buffer(1000)

# Dissolve overlapping buffers into single polygon
service_area = buffers.dissolve()
service_area_wgs84 = service_area.to_crs("EPSG:4326")
service_area_wgs84.to_file("service_area.geojson", driver="GeoJSON")
print(f"Service area: {service_area_wgs84.geometry.area.sum():.0f} sq degrees")

故障排除

问题	原因	解决方案
空间连接结果为空	GeoDataFrame 之间的 CRS 不匹配	确保 CRS 匹配：gdf2 = gdf2.to_crs(gdf1.crs)
面积／距离值错误	使用了地理 CRS（度）	重新投影到投影 CRS：gdf.to_crs(gdf.estimate_utm_crs())
读取时出现 DriverError	缺少驱动程序或文件损坏	检查文件是否存在；尝试显式指定 driver="GeoJSON"
合并后几何列丢失	调用了 df.merge(gdf)，而不是 gdf.merge(df)	始终在 GeoDataFrame 上调用 merge
overlay 时出现 TopologicalError	几何无效	使用 gdf.geometry = gdf.geometry.buffer(0) 修复
sjoin_nearest 速度慢	大型数据集未设置距离限制	设置 max_distance 参数
Folium 地图空白	几何不在 EPSG:4326 中	重新投影到 WGS84：gdf.to_crs(epsg=4326).explore()
Shapefile 列名被截断	Shapefile 有 10 个字符的限制	改用 GeoPackage：gdf.to_file("out.gpkg")

参考资料

• GeoPandas 文档：https://geopandas.org/
• GeoPandas GitHub：https://github.com/geopandas/geopandas
• Shapely 文档：https://shapely.readthedocs.io/
• EPSG 注册表：https://epsg.io/

相关技能

• matplotlib-scientific-plotting — 高级地图样式设置和图像导出
• polars-dataframes — 空间操作前后的高性能表格分析
• folium — 超越 geopandas.explore() 的高级交互式 Web 制图