第18章：与第三方库集成

本章讲解 Shapely 与 Python GIS 生态中其他重要库的集成方式，包括 GeoPandas、Fiona、matplotlib、NumPy、pyproj、rasterio 等。Shapely 作为几何计算的核心引擎，与这些库的良好配合是构建完整 GIS 工作流的关键。

18.1 与 GeoPandas 集成

18.1.1 GeoPandas 简介

GeoPandas 是基于 Pandas 的地理空间数据处理库，其核心 GeoDataFrame 的 geometry 列直接使用 Shapely 几何对象。

18.1.2 从 Shapely 到 GeoDataFrame

import geopandas as gpd
from shapely import Point, Polygon, LineString
import pandas as pd

# 从 Shapely 几何体创建 GeoDataFrame
data = {
    'name': ['北京', '上海', '广州'],
    'population': [2171, 2487, 1868],
    'geometry': [
        Point(116.46, 39.92),
        Point(121.47, 31.23),
        Point(113.26, 23.13),
    ]
}

gdf = gpd.GeoDataFrame(data, crs="EPSG:4326")
print(gdf)
print(f"\n几何列类型: {type(gdf.geometry.iloc[0])}")

18.1.3 GeoDataFrame 中的 Shapely 操作

import geopandas as gpd
from shapely import Point

# 创建城市点数据
cities = gpd.GeoDataFrame({
    'name': ['A', 'B', 'C'],
    'geometry': [Point(0, 0), Point(3, 4), Point(6, 0)]
}, crs="EPSG:4326")

# 使用 Shapely 的 buffer 操作
cities['buffer_500m'] = cities.geometry.buffer(0.005)  # 约 500m

# 使用 Shapely 的 centroid
print("质心:", cities.geometry.values[0])

# 计算两两距离
for i, city1 in cities.iterrows():
    for j, city2 in cities.iterrows():
        if i < j:
            dist = city1.geometry.distance(city2.geometry)
            print(f"{city1['name']} → {city2['name']}: {dist:.4f}")

18.1.4 空间连接

import geopandas as gpd
from shapely import Point, Polygon

# 点数据
points = gpd.GeoDataFrame({
    'id': [1, 2, 3, 4],
    'geometry': [Point(1, 1), Point(3, 3), Point(5, 5), Point(7, 2)]
})

# 面数据
polygons = gpd.GeoDataFrame({
    'zone': ['A', 'B'],
    'geometry': [
        Polygon([(0, 0), (4, 0), (4, 4), (0, 4)]),
        Polygon([(4, 0), (8, 0), (8, 4), (4, 4)]),
    ]
})

# 空间连接（点在哪个区域内）
joined = gpd.sjoin(points, polygons, how='left', predicate='within')
print(joined[['id', 'zone']])

18.1.5 批量几何操作

import geopandas as gpd
from shapely import Polygon
from shapely.ops import unary_union

# 创建多个地块
parcels = gpd.GeoDataFrame({
    'type': ['住宅', '商业', '住宅', '工业', '商业'],
    'geometry': [
        Polygon([(0, 0), (2, 0), (2, 2), (0, 2)]),
        Polygon([(2, 0), (4, 0), (4, 2), (2, 2)]),
        Polygon([(0, 2), (2, 2), (2, 4), (0, 4)]),
        Polygon([(2, 2), (4, 2), (4, 4), (2, 4)]),
        Polygon([(4, 0), (6, 0), (6, 2), (4, 2)]),
    ]
})

# 按类型分组合并
for ptype, group in parcels.groupby('type'):
    merged = unary_union(group.geometry.values)
    print(f"{ptype}: 面积={merged.area:.2f}")

# 计算总面积
print(f"总面积: {parcels.geometry.area.sum():.2f}")

18.2 与 Fiona 集成

18.2.1 使用 shape() 和 mapping()

Shapely 提供 shape() 和 mapping() 函数与 Fiona（以及任何 GeoJSON 兼容接口）交互：

from shapely.geometry import shape, mapping
from shapely import Polygon

# GeoJSON 字典 → Shapely 几何
geojson_dict = {
    "type": "Polygon",
    "coordinates": [[(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]]
}

geom = shape(geojson_dict)
print(f"类型: {geom.geom_type}")
print(f"面积: {geom.area}")

# Shapely 几何 → GeoJSON 字典
poly = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
geojson = mapping(poly)
print(f"GeoJSON: {geojson}")

18.2.2 读取 Shapefile

import fiona
from shapely.geometry import shape

# 读取 Shapefile
# with fiona.open("buildings.shp") as src:
#     for feature in src:
#         geom = shape(feature['geometry'])
#         props = feature['properties']
#         print(f"名称: {props.get('name')}")
#         print(f"面积: {geom.area:.2f}")
#         print(f"类型: {geom.geom_type}")

# 示例：模拟 Fiona 输出
features = [
    {
        'geometry': {'type': 'Polygon', 'coordinates': [[(0,0),(1,0),(1,1),(0,1),(0,0)]]},
        'properties': {'name': '建筑A', 'floors': 5}
    },
    {
        'geometry': {'type': 'Polygon', 'coordinates': [[(2,0),(3,0),(3,1),(2,1),(2,0)]]},
        'properties': {'name': '建筑B', 'floors': 3}
    }
]

for feature in features:
    geom = shape(feature['geometry'])
    print(f"{feature['properties']['name']}: 面积={geom.area:.2f}")

18.2.3 写入 GeoJSON 文件

from shapely.geometry import mapping
from shapely import Point, Polygon
import json

# 创建 GeoJSON FeatureCollection
features = []

geometries = [
    ("公园A", Point(116.4, 39.9).buffer(0.01)),
    ("公园B", Point(116.5, 39.8).buffer(0.015)),
]

for name, geom in geometries:
    feature = {
        "type": "Feature",
        "geometry": mapping(geom),
        "properties": {"name": name, "area": geom.area}
    }
    features.append(feature)

geojson = {
    "type": "FeatureCollection",
    "features": features
}

# 输出 GeoJSON
print(json.dumps(geojson, indent=2, ensure_ascii=False)[:500])

18.3 与 matplotlib 可视化

18.3.1 基本绘图

import matplotlib.pyplot as plt
from shapely import Point, Polygon, LineString
from shapely.plotting import plot_polygon, plot_line, plot_points

# 创建几何体
polygon = Polygon([(0, 0), (4, 0), (4, 3), (0, 3)])
line = LineString([(0, 0), (2, 3), (4, 0)])
point = Point(2, 1.5)

fig, axes = plt.subplots(1, 3, figsize=(12, 4))

# 绘制多边形
plot_polygon(polygon, ax=axes[0], add_points=False, color='lightblue', alpha=0.5)
axes[0].set_title('多边形')

# 绘制线
plot_line(line, ax=axes[1], add_points=True, color='red')
axes[1].set_title('线')

# 绘制点
axes[2].plot(point.x, point.y, 'ko', markersize=10)
axes[2].set_title('点')

for ax in axes:
    ax.set_aspect('equal')

plt.tight_layout()
plt.savefig('/tmp/shapely_plot.png', dpi=100)
plt.close()

18.3.2 使用 Matplotlib Patches

import matplotlib.pyplot as plt
from matplotlib.patches import Polygon as MplPolygon
from matplotlib.collections import PatchCollection
from shapely import Polygon, MultiPolygon
import numpy as np

def shapely_to_mpl_patch(polygon, **kwargs):
    """将 Shapely Polygon 转换为 Matplotlib Patch"""
    exterior = np.array(polygon.exterior.coords)
    patch = MplPolygon(exterior, **kwargs)
    return patch

# 创建多个多边形
polys = [
    Polygon([(0, 0), (2, 0), (2, 2), (0, 2)]),
    Polygon([(3, 0), (5, 0), (5, 2), (3, 2)]),
    Polygon([(1, 3), (4, 3), (4, 5), (1, 5)]),
]

fig, ax = plt.subplots(figsize=(8, 6))

colors = ['#ff6b6b', '#4ecdc4', '#45b7d1']
for poly, color in zip(polys, colors):
    patch = shapely_to_mpl_patch(poly, facecolor=color, edgecolor='black', alpha=0.7)
    ax.add_patch(patch)

ax.set_xlim(-1, 6)
ax.set_ylim(-1, 6)
ax.set_aspect('equal')
ax.set_title('Shapely + Matplotlib')
plt.savefig('/tmp/shapely_patches.png', dpi=100)
plt.close()

18.3.3 缓冲区可视化

import matplotlib.pyplot as plt
from shapely import Point, LineString
from shapely.plotting import plot_polygon, plot_line

fig, ax = plt.subplots(figsize=(8, 8))

# 点缓冲区
center = Point(5, 5)
for dist in [1, 2, 3]:
    buf = center.buffer(dist)
    plot_polygon(buf, ax=ax, add_points=False, 
                 color=f'C{dist}', alpha=0.2)
ax.plot(center.x, center.y, 'ko', markersize=5)

# 线缓冲区
line = LineString([(0, 0), (3, 2), (6, 1)])
buf = line.buffer(0.5)
plot_polygon(buf, ax=ax, add_points=False, color='yellow', alpha=0.3)
plot_line(line, ax=ax, add_points=False, color='red')

ax.set_aspect('equal')
ax.set_title('缓冲区可视化')
plt.savefig('/tmp/shapely_buffer_viz.png', dpi=100)
plt.close()

18.4 与 NumPy 集成

18.4.1 坐标数组操作

import numpy as np
from shapely import get_coordinates, set_coordinates, points
from shapely import Point, LineString, Polygon

# 从几何体提取坐标到 NumPy 数组
polygon = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
coords = get_coordinates(polygon)
print(f"坐标数组形状: {coords.shape}")
print(coords)

# 修改坐标生成新几何体
new_coords = coords * 2  # 所有坐标乘以 2
new_polygon = set_coordinates(polygon, new_coords)
print(f"原面积: {polygon.area}, 新面积: {new_polygon.area}")

18.4.2 向量化几何创建

import numpy as np
from shapely import points, linestrings, polygons, box

# 批量创建点
coords = np.random.uniform(0, 100, (1000, 2))
pts = points(coords)
print(f"创建 {len(pts)} 个点")

# 批量创建矩形
boxes = [box(i, 0, i+1, 1) for i in range(10)]
print(f"创建 {len(boxes)} 个矩形")

# 向量化属性计算
import shapely
areas = shapely.area(boxes)
print(f"面积: {areas}")

18.4.3 NumPy 数组中的几何体

import numpy as np
import shapely
from shapely import Point

# Shapely 2.0+ 支持几何体的 NumPy 数组
geom_array = np.array([
    Point(0, 0), Point(1, 1), Point(2, 2), Point(3, 3)
])

# 向量化操作
buffers = shapely.buffer(geom_array, 0.5)
print(f"缓冲区面积: {shapely.area(buffers)}")

# 向量化谓词
target = Point(1.5, 1.5).buffer(2)
results = shapely.contains(target, geom_array)
print(f"包含结果: {results}")

18.5 与 pyproj 坐标转换

18.5.1 基本坐标转换

from shapely import Point, Polygon
from shapely.ops import transform
from pyproj import Transformer

# 创建坐标转换器（WGS84 → Web Mercator）
transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)

# 转换点
point_wgs84 = Point(116.46, 39.92)  # 北京（经纬度）

# 使用 transform 函数
point_mercator = transform(transformer.transform, point_wgs84)
print(f"WGS84: {point_wgs84}")
print(f"Web Mercator: {point_mercator}")

18.5.2 多边形坐标转换

from shapely import Polygon
from shapely.ops import transform
from pyproj import Transformer

transformer = Transformer.from_crs("EPSG:4326", "EPSG:32650", always_xy=True)

# WGS84 多边形
polygon_wgs84 = Polygon([
    (116.0, 39.5), (117.0, 39.5),
    (117.0, 40.5), (116.0, 40.5)
])

# 转换到 UTM Zone 50N
polygon_utm = transform(transformer.transform, polygon_wgs84)
print(f"WGS84 面积（度²）: {polygon_wgs84.area:.6f}")
print(f"UTM 面积（m²）: {polygon_utm.area:.0f}")
print(f"UTM 面积（km²）: {polygon_utm.area / 1e6:.2f}")

18.5.3 使用 transform_coordseq（Shapely 2.1+）

from shapely import transform as shapely_transform
from shapely import Point, LineString
from pyproj import Transformer
import numpy as np

transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)

def proj_transform(coords):
    """坐标转换函数，接收 (N, 2) 数组"""
    x, y = transformer.transform(coords[:, 0], coords[:, 1])
    return np.column_stack([x, y])

line_wgs84 = LineString([(116, 39), (117, 40), (118, 39)])

# 注意：Shapely 2.1+ 的 shapely.transform 函数
# line_mercator = shapely_transform(line_wgs84, proj_transform)

18.6 与 rasterio 集成

18.6.1 几何体作为掩膜

# 使用 Shapely 几何体裁剪栅格数据
from shapely import Polygon, box

# 创建掩膜几何体
mask_polygon = Polygon([
    (116.0, 39.5), (117.0, 39.5),
    (117.0, 40.5), (116.0, 40.5)
])

# 与 rasterio 配合使用
# import rasterio
# from rasterio.mask import mask as rio_mask
#
# with rasterio.open("dem.tif") as src:
#     out_image, out_transform = rio_mask(
#         src,
#         [mapping(mask_polygon)],
#         crop=True
#     )
#     print(f"裁剪后形状: {out_image.shape}")

# 创建边界框用于窗口读取
bbox = box(116.0, 39.5, 117.0, 40.5)
print(f"边界框: {bbox.bounds}")

18.7 与 scipy 空间分析

18.7.1 KDTree vs STRtree

import numpy as np
from shapely import STRtree, points

# 使用 Shapely STRtree
np.random.seed(42)
coords = np.random.uniform(0, 100, (1000, 2))
geom_points = points(coords)

tree = STRtree(geom_points)

# 查询最近邻
from shapely import Point
query_point = Point(50, 50)
nearest_idx = tree.nearest(query_point)
print(f"STRtree 最近邻索引: {nearest_idx}")

# scipy KDTree 对比
# from scipy.spatial import KDTree
# kdtree = KDTree(coords)
# dist, idx = kdtree.query([50, 50])
# print(f"KDTree 最近邻: 索引={idx}, 距离={dist:.4f}")

18.8 GeoJSON 与 Web 集成

18.8.1 生成 Leaflet/Folium 兼容数据

from shapely import Point, Polygon
from shapely.geometry import mapping
import json

# 创建几何体
features = []

# 添加点要素
for name, lon, lat in [("北京", 116.46, 39.92), ("上海", 121.47, 31.23)]:
    features.append({
        "type": "Feature",
        "geometry": mapping(Point(lon, lat)),
        "properties": {"name": name}
    })

# 添加面要素
area = Polygon([(116, 39), (117, 39), (117, 40), (116, 40)])
features.append({
    "type": "Feature",
    "geometry": mapping(area),
    "properties": {"name": "区域A"}
})

geojson = {"type": "FeatureCollection", "features": features}
print(json.dumps(geojson, indent=2, ensure_ascii=False)[:500])

# 可直接用于 Folium
# import folium
# m = folium.Map(location=[35, 115], zoom_start=5)
# folium.GeoJson(geojson).add_to(m)
# m.save("map.html")

18.9 与 PostGIS 数据库交互

18.9.1 通过 psycopg2 读写几何体

from shapely import Point, Polygon
from shapely import wkb, wkt

# 写入 PostGIS
polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
wkb_hex = polygon.wkb_hex

# SQL 示例
sql_insert = f"""
INSERT INTO parcels (name, geom)
VALUES ('地块A', ST_GeomFromWKB(decode('{wkb_hex}', 'hex'), 4326))
"""
print("插入 SQL:")
print(sql_insert[:100] + "...")

# 从 PostGIS 读取
# cursor.execute("SELECT ST_AsBinary(geom) FROM parcels")
# for row in cursor:
#     geom = wkb.loads(row[0])
#     print(f"类型: {geom.geom_type}, 面积: {geom.area}")

18.9.2 通过 SQLAlchemy + GeoAlchemy2

# SQLAlchemy + GeoAlchemy2 集成示例
# from sqlalchemy import create_engine
# from geoalchemy2 import Geometry
# from shapely import Point
#
# engine = create_engine("postgresql://user:pass@localhost/gisdb")
#
# # 使用 WKT 格式
# point = Point(116.46, 39.92)
# sql = f"INSERT INTO cities (name, geom) VALUES ('北京', ST_GeomFromText('{point.wkt}', 4326))"

18.10 集成最佳实践

18.10.1 推荐的库组合

18.10.2 数据流转模式

数据读取                    核心处理                    输出
┌──────────┐           ┌────────────┐           ┌──────────┐
│ Fiona    │──shape()─→│  Shapely   │──mapping()─→│ GeoJSON  │
│ GeoPandas│           │  几何操作   │            │ Fiona    │
│ PostGIS  │──wkb──→   │  空间分析   │──wkb──→    │ PostGIS  │
│ WKT/WKB  │           │  缓冲/裁剪  │            │ matplotlib│
└──────────┘           └────────────┘           └──────────┘

18.10.3 注意事项

1. CRS 一致性：Shapely 不处理坐标参考系统，确保在操作前统一 CRS
2. 单位一致性：面积、距离计算结果取决于坐标单位（度 vs 米）
3. 几何有效性：在库间传递几何体时验证有效性
4. 版本兼容：确保 Shapely 2.x 与其他库的版本兼容
5. 内存管理：大数据集使用流式处理，避免一次性加载所有几何体