Pythonで出発点と進む距離・方向から到着点の緯度経度を計算する方法
はじめに
以前『火山灰追跡モデルPUFF』についてnoteを書きました。原論文を見ながらゼロからPythonで実装するという、なかなかハードで良い経験でした。
このとき火山灰の動きを計算する方法として、Pythonのpyprojライブラリを活用しました。pyprojは地球が楕円体であることを考慮して、地理的な位置情報を使った計算ができるライブラリです。
pyprojの使い方をweb検索したところ、ある２地点の緯度経度から距離を計算する方法について書かれたweb記事等はたくさん見つかりました。
例えば大阪管区気象台と東京管区気象台の距離は、以下のコードで計算できます。
import pyproj
grs80 = pyproj.Geod(ellps='GRS80')
# 東京管区気象台
tokyo_lat = 35.69
tokyo_lon = 139.75
# 大阪管区気象台
osaka_lat = 34.68
osaka_lon = 135.52
grs80.inv(osaka_lon, osaka_lat, tokyo_lon, tokyo_lat)# 大阪から見た方位角、東京から見た方位角、２地点間の距離(m)
(72.5716380992533, -104.99012020591334, 401201.948838876)
つまり、大阪管区気象台から東京管区気象台は、約73度の方向に約401kmの距離のところにある、ということがわかります。
（東京管区気象台から見ると、約255度の方向に大阪管区気象台がある）
こんな感じで位置情報を扱うときに非常に便利なライブラリですが、今回私が必要としているのは、ある緯度経度の出発点からある方向にある距離だけ進んだ先の到着点の緯度経度を計算する方法であり、これについて書かれたweb記事はほとんど見つかりませんでした。
結局、pyprojの公式ドキュメントを参照して該当するメソッドを見つけました。そこで、どれだけニーズがあるのかわかりませんが、web記事無いなら自分で書いてしまおうということで、備忘録にもなるのでnoteにまとめたいと思います。
やりたかったこと
今回『ラグランジュ的に大気中の粒子の動きをシミュレーションしたい』というのがモチベーションでした。
PUFFモデルでは、火山灰粒子が大気中を風で運ばれることを考えます。ここで風速と時間ステップ（5分とか10分とか）があれば、進む距離が計算できます。また風向があれば、その反対方向が進む方向(※1)となります。
(※1)
風向は風が吹いてくる方向で表すので、火山灰粒子が運ばれる方向＝風が吹いていく方向は風向の反対方向となります。
そうすると、始点の緯度経度と進む方向（方位角）と距離を与えれば、終点の緯度経度を計算してくれるメソッドが必要になります。
地球を球面座標(※2)で近似してしまえば自力で実装することはそれほど難しくないのですが、地球楕円体を真面目に扱ってみたいと思い、それならpyprojに何か良いメソッドがあるはず、と思って調べてみたら案の定ありました。
ではその秘蔵の（？）やり方を次章でご紹介いたします。
(※2)
球面座標でも実装し、実際に計算を行ってみました。その結果、pyprojで実装した場合と計算結果に大きな違いはありませんでした。PUFFで扱う時間・空間スケールの範囲であれば、球面座標でも十分そうでした。
pyprojを使って計算する
ではやり方です。先述のコードですでにご披露していますが、まずは下準備から。
ここでは測地系としてGRS80を指定することにします。
import pyproj
grs80 = pyproj.Geod(ellps='GRS80')解説用サンプルデータとして、出発点の緯度経度を北緯30度・東経130度とします。また風向225度（南西風）の、風速30m/sとします。時間ステップは10分とします。
すると進む方向（方位角）と距離は以下のように求まります。
# 緯度経度
lat0 = 30.0
lon0 = 130.0
# 風向風速
wind_speed = 30.0
wind_dir = 225.0
# 時間ステップ（秒）
time_step = 10 * 60
# 方位角と距離
azimuth = wind_dir - 180
distance = wind_speed * time_step
print(f"出発点：{lat0:.2f}N, {lon0:.2f}E")
print(f"進む方向：{azimuth:.2f}度、進む距離：{distance:.2f}m")出発点：30.00N, 130.00E
進む方向：45.00度、進む距離：18000.00m
では本題の、出発点と方位角・距離を与えて到着点の緯度経度を計算する方法ですが、fwdメソッドを使えば簡単に求まります。
fwdメソッドには引数として出発点の経度・緯度、方位角・距離を与えます。戻り値は到着点の経度・緯度・逆方位角（到着点から出発点を見た時の方位角）です。
南西風なので北東方向に進みますから、緯度経度ともに出発点より数値は大きくなるはずです。
lon1, lat1, back_azimuth = grs80.fwd(lon0, lat0, azimuth, distance)
print(f"到着点：{lat1:.2f}N, {lon1:.2f}E")到着点：30.11N, 130.13E
確かに、出発点より到着点の緯度経度の方が大きな数値で求まりました。
なお今回は各変数とも値1つだけ与えていますが、numpy配列で複数の値を与えれば、複数の計算を一度に行うこともできます。
さてこれで良さそうですが、念のためinvメソッドで検算してみます。invメソッドについては冒頭で書いたとおり、検索すればたくさんweb記事が出てきますので、ここでは細かい説明は省略します。
fwd_azimuth, back_azimuth, distance = grs80.inv(lon0, lat0, lon1, lat1)
print(f"進む方向：{fwd_azimuth:.2f}度、進む距離：{distance:.2f}m")進む方向：45.00度、進む距離：18000.00m
方位角と距離を再現することができました。
数値予報GPVの風を使う
さてこれでpyprojの使い方としては以上なのですが、気象データ、特に数値予報GPVデータを使う場合、風のデータはUV成分（東西成分・南北成分）になっています。この場合、まずUV成分を風向風速へ変換するという作業が発生します。
そこで、このやり方もついでに記述しておきます。
ここでは私の独自Pythonモジュール『wxparams』を使います。wxparamsについては以下のnoteやQiita記事をご参照下さい。
最初にサンプルコードを掲載します。その下に解説を書きますので、ご覧下さい。
なお、数値予報GPVから風のデータを取得するところからやると大変なので、ここでは簡単にサンプルデータを与えてしまいます。
import wxparams as wx
import numpy as np
# 風のUV成分
u_wind = np.array([10.0])
v_wind = np.array([-5.0])
# UV成分から風向風速を計算する関数
distance, azimuth = wx.UV_to_SpdDir(-1.0 * u_wind * time_step, -1.0 * v_wind * time_step)
print(f"進む方向：{azimuth[0]:.2f}度、進む距離：{distance[0]:.2f}m")
# 到着点の緯度経度を求める
lon1, lat1, back_azimuth = grs80.fwd(lon0, lat0, azimuth[0], distance[0])
print(f"到着点：{lat1:.2f}N, {lon1:.2f}E")進む方向：116.57度、進む距離：6708.20m
到着点：29.97N, 130.06E
風のUV成分として、U成分（東西風）を10m/s、V成分（南北風）を-5m/sとします。U成分は西風が、V成分は南風がプラスの値となります。このUV成分の場合、西北西の風を表していることになります。
ここでwxparamsでは引数としてnumpy配列を想定していますので、値1つの場合でもnumpy配列の形にして与えます。戻り値もnumpy配列になります。
次にUV成分から風向風速を計算するUV_to_SpdDir関数を使っています。ここで引数としてUV成分に時間ステップを掛け算することで、戻り値は風速ではなく進む距離になります。またUV成分それぞれ符号を逆にすることで、戻り値が風向ではなく進む方向（方位角）にすることができます。
こうして得られた距離と方位角を使い、上述のfwdメソッドを使うことで、到着点の緯度経度を求めることができます。この例では西北西の風により東南東に進むことになるので、緯度の値は小さく、経度の値は大きくなるはずです。
実際に出力を見ると、想定通りの値で到着点の緯度経度が求められていることがわかります。
以上、pyprojとwxparamsを組み合わせてGPVデータを活用する具体例のご紹介でした。最後までお読みいただき、ありがとうございました。
※ 本noteは以前Qiitaに投稿した記事を、note用に加筆修正したものです。
この記事が気に入ったらサポートをしてみませんか？