pyproj

PyProj 是 PROJ（制图投影与坐标转换库）的 Python 接口，用于处理不同坐标参考系统（CRS）之间的转换以及椭球面上的大地测量数学计算。

用途

PyProj 适用于任何需要高精度空间分析的场景，主要用途包括：

将 GPS 坐标（WGS84）转换为本地投影系统（如 UTM、Albers、State Plane）。
在地球椭球面上计算精确的距离、方位角和面积（大地测量数学）。
定义自定义坐标参考系统（CRS）。
处理不同垂直基准面之间的转换。
处理“大圆”路径和大地线。
识别特定位置的最佳 UTM 分区。

主要能力

PyProj 提供三个核心对象：

CRS（坐标参考系统）：定义数字（坐标）如何映射到现实世界。可通过 EPSG 代码（如 EPSG:4326）、PROJ 字符串或 WKT（Well-Known Text）定义。可访问名称、适用范围、轴信息、单位等属性，并支持比较。
Transformer：针对将大量点从一个 CRS 转换到另一个 CRS 而优化的对象。应始终使用 Transformer 进行批量操作，而非一次性调用。
Geod：用于执行“椭球数学”的类。用于在不投影到平面地图的情况下计算地球上的距离。支持逆变换（获取两点间距离和方位角）、正变换（给定起点、方位角和距离求点）、中间点（大圆路径）以及椭球面上多边形的面积和周长。

使用方式

安装：

pip install pyproj

标准导入：

import pyproj
from pyproj import CRS, Transformer, Geod

基本转换模式：

from pyproj import Transformer

# 1. 定义转换（WGS84 到 Web Mercator）
# always_xy=True 确保 (lon, lat) 顺序，而非 PROJ 默认的 (lat, lon)
transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)

# 2. 转换坐标
lon, lat = -74.006, 40.7128 # 纽约
x, y = transformer.transform(lon, lat)

print(f"Projected Meters: X={x:.2f}, Y={y:.2f}")

大地测量计算示例：

from pyproj import Geod

g = Geod(ellps='WGS84')

# 1. 逆变换：获取两点间距离和方位角
# 纽约到伦敦
lon1, lat1 = -74.006, 40.7128
lon2, lat2 = -0.1278, 51.5074
az12, az21, dist = g.inv(lon1, lat1, lon2, lat2)
print(f"Distance: {dist/1000:.2f} km")

# 2. 正变换：给定起点、方位角和距离求点
# 从纽约出发，向东 1000 公里（方位角 90）
lon_new, lat_new, back_az = g.fwd(lon1, lat1, 90, 1000000)

# 3. 中间点（大圆路径）
path = g.npts(lon1, lat1, lon2, lat2, npts=10) # 路径上的 10 个点

# 4. 椭球面上多边形的面积
lons = [-10, 10, 10, -10]
lats = [-10, -10, 10, 10]
area, perimeter = g.polygon_area_perimeter(lons, lats)
print(f"Area: {abs(area)/1e6:.2f} km²")

关键规则与最佳实践

应该做的：

使用 always_xy=True：默认情况下，PROJ 6+ 使用 CRS 定义的顺序（通常是 lat, lon）。设置 always_xy=True 强制使用一致的经度、纬度（X, Y）顺序。
复用 Transformer：创建 Transformer 开销较大。创建一次，用于数据集中的所有点。
使用 Geod 计算距离：如果需要两个 GPS 点之间的距离，不要将它们投影到地图上；使用 Geod.inv() 计算椭球距离。
检查 CRS 有效性：使用 CRS.from_user_input(code).is_valid 验证投影字符串。
向量化：将 NumPy 数组或坐标列表传递给 transformer.transform()。这比循环遍历点快得多。

不应该做的：

为短距离投影：对于全局分析，不要投影到平面地图来计算距离；地图投影会引入变形（例如墨卡托使两极看起来巨大）。
硬编码 Proj4 字符串：优先使用 EPSG 代码（如 "EPSG:4326"），因为它们更稳健并包含现代基准面转换。
忽略基准面：记住在不同基准面之间移动（如 NAD27 到 WGS84）需要特定的转换参数。

反模式（切勿）：

from pyproj import Transformer, Geod

# ❌ 错误：在循环内创建 transformer（极慢！）
for lon, lat in coordinates:
    t = Transformer.from_crs("EPSG:4326", "EPSG:32633") # 初始化开销
    x, y = t.transform(lon, lat)

# ✅ 正确：初始化一次，批量转换
t = Transformer.from_crs("EPSG:4326", "EPSG:32633", always_xy=True)
lons, lats = zip(*coordinates)
xs, ys = t.transform(lons, lats)

# ❌ 错误：在平面投影上计算长距离
# x1, y1 = t.transform(lon1, lat1)
# x2, y2 = t.transform(lon2, lat2)
# dist = np.sqrt((x1-x2)**2 + (y1-y2)**2) # ❌ 结果变形！

# ✅ 正确：使用 Geod 计算大地线距离
g = Geod(ellps='WGS84')
az12, az21, dist = g.inv(lon1, lat1, lon2, lat2) # ✅ 精确米数

适用场景

空间数据科学、地理信息系统（GIS）分析。
需要高精度坐标转换的科研与工程应用。
处理大规模坐标数据集（支持向量化与迭代器转换）。
需要在地球椭球面上进行精确距离、面积、方位角计算的场景。

限制与注意事项

将点投影到 CRS 有效区域之外（例如在墨卡托中投影 90° 纬度）会产生 inf 或 nan。转换前应检查边界。
历史上 GIS 用户说“Lat/Lon”，但数学上 X=Lon, Y=Lat。使用 always_xy=True 并传入 (Lon, Lat) 以避免错误结果。
在不同参考系统之间转换时，需注意基准面转换。
对于来自生成器的大型数据集，可使用 itransform 进行内存高效处理。
Transformer 是线程安全的，可与 concurrent.futures 一起使用。

参考文档

官方文档：https://pyproj4.github.io/pyproj/
PROJ（引擎）：https://proj.org/
EPSG 注册表：https://epsg.io/（查找 CRS 代码必备）
搜索模式：pyproj.Transformer, pyproj.CRS, pyproj.Geod

技能介绍

用途

主要能力

使用方式

关键规则与最佳实践

适用场景

限制与注意事项

参考文档