[[pyproj]]

pyproj

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

主にできること

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

座標系・測地系の変換

座標系・測地系とは

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

地球上の1点を示す時、よく「緯度経度」が用いられる(経緯度)。2つの角度で地球上の位置を表し、どこであろうとおよそ正確に表せる。でも、ある2点間の距離を知りたい時などは、計算が難しくて不便だ。要求する精度にもよるが、日本国内くらいならXYで表せる「平面直角座標」を使っても実用上問題ない。このように座標系によって利点・欠点あるため、普段から我々は複数の座標系を使い分けている。

また、座標系とは別に測地系というものもある。座標から具体的な位置を特定する際、実際の地球を100%正確に表現するのは不可能だ。地球表面はごつごつして起伏があるし、さらには地殻変動で常に変化している。なので、「大体こんな感じ」という楕円体に補正を加えて近似したものを使う。この近似の仕方が、測地系である。

同じ緯度経度の数値でも、測地系(地球をどんな楕円体と見なすか)によって、指し示す位置は変わってくる。

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

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

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

座標系・測地系の変換方法

pyproj.Projクラスは、座標系・測地系を表現するクラスである。

変換は、pyproj.transform()に変換元と変換先のProjインスタンスと、X座標、Y座標を渡すだけ。XY座標は同じ長さのリストやnumpy.arrayでもよい。

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

東西方向、南北方向の順に与えるので、緯度経度の場合は「経度、緯度」の順になる。

ちなみに日本の平面直角座標は、奇妙なことに南北方向をX軸、東西方向をY軸で表すと定義されている(わかりやすい平面直角座標系|国土地理院)。一方、pyprojで返ってくるのは一貫して「東西、南北」の順である。

2点の方位角・距離の計算

楕円体を表す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

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

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

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

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

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

Tips

pandasのDataFrameを変換

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

コレクションの違いによる速度

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

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

  • NumPy (numpy.ndarray)
  • List
  • Tuple
  • Array (array.array)
  • Listを実行前にArrayに変換
NumPy List Tuple Array 変換
n=100000, m=100 16.456 18.198 18.104 16.195 17.187
n=100, m=100000 16.680 17.722 17.761 15.986 16.775
  • nは配列サイズ、mは繰返し回数
  • Arrayが一番速い
  • NumPy配列が微差で次いで速い
  • ListやTupleは若干遅い
    • 実行前にそのためだけにArrayに変換してもお釣りが来る
  • そこまでこだわる価値があるかというとそうでもない
    • 5%くらいなので、大量ループが必要な場合は一考

検証コード

programming/python/packages/pyproj.txt · 最終更新: 2018/08/10 by ikatakos
CC Attribution 4.0 International
Driven by DokuWiki Recent changes RSS feed Valid CSS Valid XHTML 1.0