pyproj

• pyproj Documentation
• GitHub - pyproj4/pyproj: Python interface to PROJ

地理的な位置を示す情報を扱うときに、座標系・測地系変換を行ったり、2点間の距離・方位角を計算したりできる。

• 座標系・測地系の変換
• 緯度経度からの2点間の距離・方位角計算

• 座標系とは？ | 空間情報クラブ
• 測地系とは？ | 空間情報クラブ

ざっくり言うと、座標系とは、ある位置を示す数字の組（たとえば $(35.16, 134.27)$ とか）があったとして、 「位置を特定するには、数字をどう使ったらいいか」というルールのことだ。

地球上の1点を示す時、よく座標系として「緯度経度」が用いられる。 「本初子午線からの東西方向の角度」「赤道からの南北方向の角度」2つの角度で地球上の位置を表し、どこであろうとおよそ正確に表せる。

でも、緯度経度だと、ある2点間の距離を知りたい時などは難しくて不便だ。 要求する精度にもよるが、日本国内くらいなら数学の座標のようにX軸Y軸で表せる「平面直角座標」を使っても実用上問題なく、距離なども計算しやすい。 このように座標系によって利点・欠点あるため、普段から複数の座標系が使い分けられている。

また、座標系とは別に測地系というものもある。
大まかには「測量の前提となるルール」のことで、いろんな場所の緯度経度を測量によって特定する際、どこを基準として、どのように算出したかを指す。 一般には、座標系のルールの中に「どんな測地系を使ったか」も含まれる。

測地系には、「地球をどんな形と見做すか」も基本的に定義に含まれる。 地球表面はごつごつして起伏があるし、さらには地殻変動で常に変化しているため、実際の地球を100%正確に表現するのは不可能だ。 なので、「大体こんな感じ」という楕円体に補正を加えて近似したものを使う。 これが異なれば、同じ緯度経度であっても指し示す位置が異なってくるため、使用した測地系も座標系には重要な要素の1つとなる。

日本では2002年の測量法改正により、測地系が日本測地系（Bessel楕円体）から世界測地系（GRS80楕円体）に変更され、誤差の大きいところでは500m以上のズレが生じることとなった。現在は世界測地系への移行が進んでいるが、未だに日本測地系での情報もあるため、注意の上、必要なら変換せねばならない。

座標系間を変換する計算式はあるが、かなり複雑。（例:平面直角座標への換算 計算式）

こんな計算を実装してられないので、PROJというオープンソースのライブラリが存在する。これのPython版モジュールが、pyprojである。複雑な計算を楽に、正確にできるので、使わない手は無い。

もう1つ、座標系の関係で出てくる用語として、CRS（Coordinate Reference System、座標参照系）がある。特にGISを使うとよく出てくる。

だいたい「座標系」と同じ説明、同じ使われ方をされているように見える。

座標系は数学など一般に使われる用語でもあるため、区別したい場合は CRS と言えば地理的な座標系であることが明確になるくらい？

• Transformer — pyproj documentation

用意するのは、変換元座標系、変換先座標系、変換対象の座標。

まず pyproj.Transformer という座標変換用クラスに変換元・変換先座標系を与えてインスタンス化し、 その後 Transformer.transform() に変換対象の座標を与えるのが基本的な流れとなる。

座標系の与え方は様々あるが、以下ではEPSGコードの文字列で与えている。

変換対象の座標は、同じ長さのリストやnumpy.arrayでもよい。複数座標を一度に変換できる。

import pyproj

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

transformer.transform(135.5, 39.5)
# => (15083791.002488568, 4793547.459104806)

transformer.transform([135.5, 135.6, 135.7], [39.5, 39.6, 39.7])
# => ([15083791.002488568, 15094922.951567894, 15106054.900647221],
#     [ 4793547.459104806,  4807984.493190501,  4822442.387442776])

pyproj ver.1.xx の記述

import pyproj

wgs84 = pyproj.Proj(init='EPSG:4326') # WGS84 緯度経度
rect6 = pyproj.Proj(init='EPSG:2448') # 平面直角座標6系

lat, lng = 34.705185, 135.498468 # 梅田駅

# WGS84緯度経度から、平面直角座標6系のXY座標に変換
x, y = pyproj.transform(wgs84, rect6, lng, lat)

print(x, y)
# => -45943.1565325506 -143527.06696254978

• Geod — pyproj documentation

楕円体を表すpyproj.Geodインスタンスを作ると、Geod.inv(from_Lng, from_Lat, to_Lng, to_Lat)で「方位角」「逆方向の方位角」「距離」が一度に算出される。方位角が両方向なのは、球面上では正方向の方位角の180°逆向きが逆方向の方位角とはならないため。

方位角は、真北を0とした時計回りのdegreeで表される。引数にradians=Trueを渡すとradianで返される。

与える緯度経度は、長さの等しいリストまたはnumpy.arrayでもよい。複数の座標対をまとめて高速に計算してくれる。

import pyproj

grs80 = pyproj.Geod(ellps='GRS80')  # GRS80楕円体

lat1, lng1 = 34.705185, 135.498468  # 梅田駅
lat2, lng2 = 35.170897, 136.881537  # 名古屋駅

azimuth, bkw_azimuth, distance = grs80.inv(lng1, lat1, lng2, lat2)

print(azimuth, bkw_azimuth, distance)
# => 67.36528283503415 -111.84261446690202 136506.23110143896

• EPSGコードとは？ | 空間情報クラブ
• Spatial Reference List -- Spatial Reference

座標系の意味合いは既述の通りだが、この組み合わせは歴史的経緯や精度の問題で、種類がいっぱいあるため、ややこしい。

それを一発で言い表せるようにまとめられたのが、EPSGコード（European Petroleum Survey Group）。

CRSに個別のコードを紐付け、番号だけで「何の測地系を使い」「緯度経度か、平面座標か」「平面座標なら原点はどこか、どのように投影したか」などが一意になる。

pyprojでは、Projオブジェクトを作る際、これを渡すのがやりやすい。

bessel = pyproj.Proj(init='EPSG:4301') # 日本測地系 緯度経度
grs80 = pyproj.Proj(init='EPSG:4612')  # 世界測地系 緯度経度
wgs84 = pyproj.Proj(init='EPSG:4326')  # WGS84 緯度経度

日本でよく使われる代表的なところは、下記のリンクにまとめられている。

• GISのための測地成果、測地系、楕円体、投影座標系、EPSGコードのまとめ - 自然環境保全のための周辺技術

pandasのデータフレームに緯度経度が入っている場合、pyprojはpandas.Seriesには残念ながら対応していないため、Series.valuesでNumpyの配列に変換してから渡す必要がある。

import pandas as pd
import pyproj

df = pd.read_csv('input.csv')

bessel = pyproj.Proj('+init=EPSG:4301')
wgs84 = pyproj.Proj('+init=EPSG:4326')
transformed = pyproj.transform(bessel, wgs84, row['LNG'].values, row['LAT'].values)

結果はNumpy配列が2つ入ったタプルで返る。（変換に高度も含む場合は3つ）

タプルなので、これまた直接はpandas.DataFrameに代入は出来ない。一度分解して、1カラムずつ与えてやる。

# ..略

# ×
df[['x', 'y']] = pyproj.transform(bessel, wgs84, row['lng'].values, row['lat'].values)

# ○
x, y = pyproj.transform(bessel, wgs84, row['lng'].values, row['lat'].values)
df['x'] = x
df['y'] = y

• 参考
  • python - Coordinates from UTM to Latitude and Longitude in pandas - Stack Overflow

Geod.inv()で複数の距離・方位角計算をするときに、渡すコレクションの違いによる速度差。（検証コードは後述）

inv()は、以下の配列であれば渡したのと同じ配列に変換して返してくれるが、内部ではArrayに一度変換しているようだ。

• NumPy (numpy.ndarray)
• List
• Tuple
• Array (array.array)
• Listを実行前にArrayに変換

• nは配列サイズ、mは繰返し回数
• Arrayが一番速い
• NumPy配列が微差で次いで速い
• ListやTupleは若干遅い
  • 実行前にそのためだけにArrayに変換してもお釣りが来る
• そこまでこだわる価値があるかというとそうでもない
  • 5%くらいなので、大量ループが必要な場合は一考

検証コード

import numpy as np
import pyproj
import timeit
import array

n = 100000
m = 100
inv = pyproj.Geod(ellps='GRS80').inv
cmd = 'inv(from_lng, from_lat, to_lng, to_lat)'

from_lng = np.random.rand(n) * 360 - 180
from_lat = np.random.rand(n) * 180 - 90
to_lng = np.random.rand(n) * 360 - 180
to_lat = np.random.rand(n) * 180 - 90

print(timeit.timeit(cmd, globals=globals(), number=m))

from_lng = from_lng.tolist()
from_lat = from_lat.tolist()
to_lng = to_lng.tolist()
to_lat = to_lat.tolist()

print(timeit.timeit(cmd, globals=globals(), number=m))

from_lng = tuple(from_lng)
from_lat = tuple(from_lat)
to_lng = tuple(to_lng)
to_lat = tuple(to_lat)

print(timeit.timeit(cmd, globals=globals(), number=m))

from_lng = array.array('d', from_lng)
from_lat = array.array('d', from_lat)
to_lng = array.array('d', to_lng)
to_lat = array.array('d', to_lat)

print(timeit.timeit(cmd, globals=globals(), number=m))

from_lng = from_lng.tolist()
from_lat = from_lat.tolist()
to_lng = to_lng.tolist()
to_lat = to_lat.tolist()

cmd = 'inv(array.array("d", from_lng), array.array("d", from_lat),' \
      '    array.array("d", to_lng), array.array("d", to_lat))'

print(timeit.timeit(cmd, globals=globals(), number=m))