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（Well-Known Text）で定義できます。

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 を使用する - 2 つの 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 を使用します
2. Transformer は一度だけ作成し、一括処理で再利用します
3. 楕円体上の距離計算には、平面座標への投影ではなく Geod を使用します
4. 堅牢性と最新の測地系への対応を高めるため、Proj4 文字列より EPSG コードを優先します
5. 個々の点をループで処理するのではなく、配列やリストを渡して変換をベクトル化します
6. 座標を変換する前に、CRS の有効性と適用範囲を確認します
7. 異なる基準系間で変換する場合は、測地系変換に注意します
8. ジェネレーターから大規模なデータセットをメモリ効率よく処理するには、itransform を使用します

PyProj は、座標精度を確保するための決定版ツールです。地球の複雑な球面三角法や楕円体三角法を抽象化することで、空間計算の地理的な妥当性を保ちながら、科学者がデータに専念できるようにします。