geopandas

GeoPandas - 地理空间数据分析

GeoPandas 让你可以使用熟悉的 Pandas API 执行空间连接、几何操作和坐标转换。它将“几何”视为 DataFrame 中的另一列，但这列数据能够计算面积、距离和交集。

适用场景

• 读取和写入空间格式（Shapefile、GeoJSON、GeoPackage、KML）。
• 执行空间连接（例如，“哪些点位于这个多边形内？”）。
• 进行坐标系转换（将经纬度转换为米制坐标）。
• 几何分析（计算缓冲区、质心、凸包）。
• 专题制图（分级统计地图）。
• 计算空间关系（包含、重叠、接触、位于内部）。
• 处理 OpenStreetMap 数据或卫星衍生的矢量数据。

参考文档

官方文档：https://geopandas.org/
交互式教程：https://geopandas.org/en/stable/gallery/index.html
搜索模式：gpd.read_file、gdf.to_crs、gpd.sjoin、gdf.buffer、gdf.explore

核心原则

GeoDataFrame

GeoDataFrame 是至少包含一个 GeoSeries 列（通常名为 geometry）的 pandas.DataFrame。每一行代表一个要素（点、线或多边形）。

坐标参考系（CRS）

没有 CRS 的数据只是网格上的数字。要进行实际的地理计算（例如以 km² 为单位计算面积），必须定义 CRS（例如 WGS84 - EPSG:4326 或 UTM）。

谓词和集合运算

空间分析依赖二元谓词（相交、位于内部、包含）和集合论运算（并集、交集、差集）。

快速参考

安装

pip install geopandas pyarrow pyproj fiona shapely

常用导入

import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Point, LineString, Polygon

基本模式 - 加载和绘图

import geopandas as gpd

# Load built-in dataset
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

# Filter and Project
europe = world[world.continent == 'Europe']
europe = europe.to_crs(epsg=3035) # Equal Area projection for Europe

# Plot
europe.plot(column='pop_est', legend=True, cmap='viridis')

关键规则

✅ 应该做

• 始终检查 CRS - 执行任何空间操作前，先验证 gdf.crs。
• 投影后再测量 - 计算面积或距离前，使用 UTM 等投影 CRS（米/英尺）。
• 使用空间索引 - 对于大型数据集，使用 gdf.sindex，或确保使用 sjoin 来加快查询。
• 验证几何 - 使用 gdf.is_valid 查找有问题的多边形（自相交）。
• 为可视化简化几何 - 使用 gdf.simplify() 加快复杂边界的绘图。
• 使用 .explore() - 在 Jupyter 中快速创建交互式地图（使用 Leaflet/Folium）。

❌ 不应该做

• 不要以度为单位测量面积 - 绝不要在经纬度 CRS（EPSG:4326）上计算 .area。结果将以“平方度”为单位（毫无意义）。
• 不要用循环迭代 - 避免逐行循环；使用向量化空间操作。
• 不要忽略拓扑关系 - 注意“接触”和“相交”是不同的（边界与内部）。
• 不要忘记设置活动几何 - 如果 GeoDataFrame 有多个几何列，请使用 gdf.set_geometry() 指定要使用的几何列。

反模式（绝不要这样做）

# ❌ BAD: Manual distance calculation on Lat/Lon (ignores Earth's curvature)
def dist(p1, p2):
    return np.sqrt((p1.x - p2.x)**2 + (p1.y - p2.y)**2)

# ✅ GOOD: Reproject and use vectorized distance
gdf = gdf.to_crs(epsg=3857) # Web Mercator (meters)
distances = gdf.distance(other_point)

# ❌ BAD: Manually checking points in polygons
for i, poly in countries.iterrows():
    for j, pt in cities.iterrows():
        if poly.geometry.contains(pt.geometry):
            print("Found")

# ✅ GOOD: Spatial Join (Optimized with spatial index)
cities_with_country = gpd.sjoin(cities, countries, predicate='within')

几何创建和操作

根据坐标创建

# From a Pandas DataFrame with Lat/Lon
df = pd.DataFrame({'City': ['NY', 'London'], 'Lat': [40.7, 51.5], 'Lon': [-74.0, -0.1]})
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.Lon, df.Lat))
gdf.set_crs(epsg=4326, inplace=True)

几何操作

# Centroids and Envelopes
gdf['centroid'] = gdf.centroid
gdf['envelope'] = gdf.envelope # Bounding box

# Buffering (Creating a zone around features)
# Warning: Do this in a projected CRS (meters)
stations_buffered = metro_stations.to_crs(epsg=32633).buffer(500) # 500 meters

# Unifying overlapping polygons
total_area = gdf.union_all()

空间查询

空间连接（sjoin）

# Find which district each school belongs to
schools_in_districts = gpd.sjoin(schools, districts, how="inner", predicate="within")

# sjoin types:
# 'intersects' (default), 'contains', 'within', 'touches', 'overlaps', 'crosses'

叠加（集合运算）

# Intersection of two layers (e.g., protected area vs. forest)
forest_in_park = gpd.overlay(forests, parks, how='intersection')

# 'union', 'intersection', 'difference', 'symmetric_difference'

坐标参考系（CRS）

重投影

# Checking the current CRS
print(gdf.crs)

# Reprojecting to a specific EPSG code
# EPSG:4326 -> WGS84 (Degrees, used for GPS)
# EPSG:3857 -> Web Mercator (Meters, used for web maps)
gdf_meters = gdf.to_crs(epsg=3857)

# Match CRS of another layer
gdf_2 = gdf_2.to_crs(gdf_1.crs)

可视化

静态和交互式地图

# Layered plotting
fig, ax = plt.subplots(figsize=(10, 10))
base = countries.plot(ax=ax, color='white', edgecolor='black')
cities.plot(ax=base, marker='o', color='red', markersize=5)

# Interactive exploration (requires folium)
cities.explore(column='population', cmap='magma', m=None)

实用工作流

1. 邻近分析（点是否位于缓冲区内）

def find_entities_near_road(roads, entities, distance_m=1000):
    """Find all entities within 1km of any road."""
    # 1. Project to a metric CRS (e.g., UTM)
    roads_m = roads.to_crs(epsg=3857)
    entities_m = entities.to_crs(epsg=3857)
    
    # 2. Create buffer
    road_buffer = roads_m.buffer(distance_m)
    
    # 3. Create a GeoDataFrame from buffer to use in sjoin
    buffer_gdf = gpd.GeoDataFrame(geometry=road_buffer, crs=roads_m.crs)
    
    # 4. Spatial Join
    nearby = gpd.sjoin(entities_m, buffer_gdf, predicate='within')
    return nearby

2. 将数据裁剪到边界范围

def clip_data(data, boundary):
    """Clip a large vector dataset to a specific boundary polygon."""
    return gpd.clip(data, boundary)

# Usage: city_parks = clip_data(national_parks, city_limits)

3. 计算面积覆盖百分比

def calculate_land_use_pct(region, land_use_layer):
    """Calculate what % of 'region' is covered by each land use type."""
    # Ensure CRS matches and is projected
    land_use_layer = land_use_layer.to_crs(region.crs)
    
    # Intersect region with land use
    intersections = gpd.overlay(land_use_layer, region, how='intersection')
    
    # Calculate area
    intersections['area'] = intersections.area
    total_area = region.area.sum()
    
    return intersections.groupby('class')['area'].sum() / total_area * 100

性能优化

使用空间索引（sindex）

# Check if a point is within any polygon in a large GDF efficiently
spatial_index = countries.sindex

# Find possible matches using bounding boxes first
possible_matches_index = list(spatial_index.intersection(target_point.bounds))
possible_matches = countries.iloc[possible_matches_index]

# Precise check only on candidates
precise_match = possible_matches[possible_matches.intersects(target_point)]

读取大型文件（Parquet）

# GeoJSON is slow to read/write. GeoParquet is significantly faster and smaller.
gdf.to_parquet("large_data.parquet")
gdf = gpd.read_parquet("large_data.parquet")

常见问题和解决方案

CRS 不匹配

# ❌ Problem: sjoin returns 0 results even if data looks overlapping
# ✅ Solution: Align CRS
if cities.crs != districts.crs:
    cities = cities.to_crs(districts.crs)

无效几何（自相交）

# ❌ Problem: Overlay or Area calculation fails/gives weird results
# ✅ Solution: Fix with buffer(0) or check validity
invalid = gdf[~gdf.is_valid]
gdf['geometry'] = gdf['geometry'].buffer(0) # Common trick to fix minor topology errors

缓冲区导致内存耗尽

# ❌ Problem: Buffering millions of points with high resolution
# ✅ Solution: Use low resolution or simplify
gdf.buffer(100, resolution=4) # Default is 16, 4 is often enough for analysis

GeoPandas 弥合了传统 GIS 软件与 Python 数据科学工具栈之间的鸿沟，让你只需编写一行 Pandas 代码，就能轻松进行空间分析。