はじめに

pyprojの使い方を検索したところ、ある２地点の緯度経度から距離を計算する方法について書かれたweb記事等はたくさん見つかるのですが、ある緯度経度の出発点からある方向にある距離だけ進んだ先の到着点の緯度経度を計算する方法について書かれたweb記事は、ほとんど見つかりませんでした。

結局、pyprojの公式ドキュメントを参照して該当するメソッドを見つけたのですが、せっかくなんでQiitaに書いておこう（自分のためにもなるし）と思ったのがきっかけです。

やりたかったこと

『ラグランジュ的に大気中の粒子の動きをシミュレーションしたい』というのがモチベーションです。

例えば火山灰のような微小粒子が大気中を風に運ばれるようなことを考えます。ここで風速と時間ステップ（5分とか30分とか1時間とか）があれば、進む距離が計算できます。また風向があれば、その反対方向が進む方向となります。

※風向は風が吹いてくる方向で表すので、微小粒子が運ばれる方向＝風が吹いていく方向は風向の反対方向となります。

そうすると、始点の緯度経度と進む方向（方位角）と距離を与えれば、終点の緯度経度を計算してくれる関数が必要です。もし地球を球面座標で近似できるなら自力で実装でもよいのですが、地球楕円体を真面目に扱うならpyprojに何か良いメソッドはないものか、と思って調べてみたら案の定ありました。

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記事が出てきますので、ここでの説明は省略します。

ret = grs80.inv(lon0, lat0, lon1, lat1)
print(f"進む方向：{ret[0]:.2f}度、進む距離：{ret[2]:.2f}m")

# [出力]
# 進む方向：45.00度、進む距離：18000.00m

方位角と距離を再現することができました。

さてこれで使い方としては以上なのですが、気象データ、特に数値予報GPVデータを使う場合、風のデータはUV成分（東西成分・南北成分）になっています。この場合、まずUV成分を風向風速へ変換するという作業が発生します。

そこで、このやり方もついでに記述しておきます。

ここでは私の独自Pythonモジュール『wxparams』を使います。wxparamsについてはこちらのQiita記事をご参照下さい。

最初にサンプルコードを掲載します。その下に解説を書きますので、ご覧下さい。

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メソッドを使うことで、到着点の緯度経度を求めることができます。この例では西北西の風により東南東に進むことになるので、緯度の値は小さく、経度の値は大きくなるはずです。

実際に出力を見ると、想定通りの値で到着点の緯度経度が求められていることがわかります。