pyproj

PyProj - 坐标投影与大地测量

PyProj 对任何需要高精度的空间分析都至关重要。它可以在数千种不同的参考系统之间转换坐标，并使用最精确的公式（Vincenty、Karney）计算地球椭球表面上的距离。

适用场景

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

参考文档

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

核心原则

CRS（坐标参考系统）

用于定义数字（坐标）如何映射到现实世界。可以通过 EPSG 代码（EPSG:4326）、PROJ 字符串或 WKT（熟知文本）定义。

Transformer

用于将大量点从一个 CRS 转换到另一个 CRS 的优化对象。批量操作时，始终使用 Transformer，而不是逐次调用。

Geod

用于执行“椭球计算”的类。在不投影到平面地图的情况下计算地球上的距离时，请使用它。

快速参考

安装

pip install pyproj

标准导入

import pyproj
from pyproj import CRS, Transformer, Geod

基本模式 - 坐标转换

from pyproj import Transformer

# 1. Define the transformation (WGS84 to Web Mercator)
# always_xy=True ensures (lon, lat) order instead of PROJ default (lat, lon)
transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)

# 2. Transform coordinates
lon, lat = -74.006, 40.7128 # NYC
x, y = transformer.transform(lon, lat)

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

关键规则

✅ 应该做

• 始终使用 always_xy=True - 默认情况下，PROJ 6+ 使用 CRS 定义的坐标顺序（通常是纬度、经度）。设置 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

# ❌ BAD: Creating a transformer inside a loop (Extremely slow!)
for lon, lat in coordinates:
    t = Transformer.from_crs("EPSG:4326", "EPSG:32633") # Initialization overhead
    x, y = t.transform(lon, lat)

# ✅ GOOD: Initialize once, transform in bulk
t = Transformer.from_crs("EPSG:4326", "EPSG:32633", always_xy=True)
lons, lats = zip(*coordinates)
xs, ys = t.transform(lons, lats)

# ❌ BAD: Calculating distance on a flat projection over long distances
# x1, y1 = t.transform(lon1, lat1)
# x2, y2 = t.transform(lon2, lat2)
# dist = np.sqrt((x1-x2)**2 + (y1-y2)**2) # ❌ Distorted result!

# ✅ GOOD: Use Geod for geodesic distance
g = Geod(ellps='WGS84')
az12, az21, dist = g.inv(lon1, lat1, lon2, lat2) # ✅ Accurate meters

使用 CRS（pyproj.CRS）

检查与比较

from pyproj import CRS

crs = CRS.from_epsg(4326)

# Accessing properties
print(crs.name)            # "WGS 84"
print(crs.area_of_use)     # Bounding box of applicability
print(crs.axis_info)       # Order of axes

# Check units
units = crs.axis_info[0].unit_name # 'degree' or 'metre'

# Comparing CRS
is_same = crs.equals(CRS.from_user_input("EPSG:4326"))

坐标转换（pyproj.Transformer）

处理大型数据集

import numpy as np
from pyproj import Transformer

# Random points in NYC area
lons = -74.0 + np.random.rand(10000)
lats = 40.7 + np.random.rand(10000)

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

# Transform entire arrays at once (Vectorized)
xs, ys = transformer.transform(lons, lats)

大地测量计算（pyproj.Geod）

距离、面积和路径

from pyproj import Geod

g = Geod(ellps='WGS84')

# 1. Inverse Transformation: Get distance and azimuth between points
# NYC to London
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. Forward Transformation: Find point given start, bearing, and distance
# Start at NYC, go 1000km East (bearing 90)
lon_new, lat_new, back_az = g.fwd(lon1, lat1, 90, 1000000)

# 3. Intermediate points (Great Circle path)
path = g.npts(lon1, lat1, lon2, lat2, npts=10) # 10 points along the path

# 4. Area of a polygon on the ellipsoid
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²")

实用工作流

1. 自动确定 UTM 分区

def get_utm_crs(lon, lat):
    """Returns the correct UTM CRS for a given lon/lat point."""
    utm_zone = int((lon + 180) / 6) + 1
    hemisphere = 'north' if lat >= 0 else 'south'
    # UTM EPSG ranges: 32601-32660 (North), 32701-32760 (South)
    base = 32600 if hemisphere == 'north' else 32700
    return CRS.from_epsg(base + utm_zone)

# Usage:
# my_crs = get_utm_crs(-74.0, 40.7) # Returns EPSG:32618 (UTM 18N)

2. 精确的点对点距离矩阵

import numpy as np
from pyproj import Geod

def calculate_distance_matrix(lons, lats):
    """Computes a symmetric distance matrix using the ellipsoid."""
    g = Geod(ellps='WGS84')
    n = len(lons)
    matrix = np.zeros((n, n))
    
    for i in range(n):
        for j in range(i + 1, n):
            _, _, d = g.inv(lons[i], lats[i], lons[j], lats[j])
            matrix[i, j] = matrix[j, i] = d
    return matrix

3. 将本地测量值转换为 GPS 坐标

def local_to_gps(start_lon, start_lat, dx_meters, dy_meters):
    """Moves a GPS point by a local offset in meters."""
    g = Geod(ellps='WGS84')
    # Move in X (East-West)
    lon_tmp, lat_tmp, _ = g.fwd(start_lon, start_lat, 90 if dx_meters > 0 else 270, abs(dx_meters))
    # Move in Y (North-South)
    lon_final, lat_final, _ = g.fwd(lon_tmp, lat_tmp, 0 if dy_meters > 0 else 180, abs(dy_meters))
    return lon_final, lat_final

性能优化

对迭代器使用 itransform

如果数据来自生成器，并且不想将其全部加载到内存中，请使用 itransform。

# Generator of (x, y) tuples
coords_gen = ((lon, lat) for lon, lat in my_source)

for x, y in transformer.itransform(coords_gen):
    process(x, y)

并行转换

Transformer 是线程安全的。可以将它们与 concurrent.futures 配合使用。

from concurrent.futures import ThreadPoolExecutor

with ThreadPoolExecutor() as executor:
    results = list(executor.map(lambda p: transformer.transform(*p), list_of_points))

常见问题与解决方案

“无穷大”问题

将 CRS 有效区域之外的点进行投影（例如，在墨卡托投影中投影纬度为 90° 的点）会得到 inf 或 nan。

# ✅ Solution: Check bounds before transforming
if not crs.area_of_use.contains(lon, lat):
    print("Point outside projection bounds!")

纬度/经度与经度/纬度

在历史上，GIS 用户会说“纬度/经度”，但从数学上讲，X=经度，Y=纬度。

# ❌ Problem: transformer.transform(40.7, -74.0) -> Wrong result!
# ✅ Solution: Use always_xy=True and pass (Lon, Lat)
transformer = Transformer.from_crs(4326, 3857, always_xy=True)
transformer.transform(-74.0, 40.7)

大地测量公式的精度

Geod 默认精度很高，但在复杂情况下若要达到亚毫米级精度：

# Use the Karney algorithm (built-in for modern pyproj)
g = Geod(ellps='WGS84') # This is already highly accurate

最佳实践

1. 创建 Transformer 时始终使用 always_xy=True，以确保 (lon, lat) 顺序一致
2. 创建 Transformer 后重复使用，以进行批量操作
3. 使用 Geod 在椭球体上计算距离，而不是投影到平面坐标
4. 优先使用 EPSG 代码，而不是 Proj4 字符串，以获得更高的稳健性并支持现代基准
5. 通过传入数组/列表而不是逐点循环来进行矢量化转换
6. 转换坐标前，检查 CRS 的有效性和适用范围
7. 在不同参考系统之间转换时，注意基准转换
8. 使用 itransform 对生成器中的大型数据集进行内存高效处理

PyProj 是坐标精度方面的权威工具。通过抽象地球复杂的球面和椭球三角计算，它让科学家能够专注于数据，同时确保空间计算在地理上有效。