🧭

Pythonで方位角を16方位に変換する

2023/09/08に公開

Python

pyproj

tech

イントロダクション

本記事では、Pythonで2地点間の方位角を16方位に変換する方法について記載します。

16方位の図（出典：日本地図センター）

こちらの記事のJavaScriptでのコードを参考にさせていただきました。

動作環境

macOS Ventura 13.5.1
Python 3.11.4

コードと説明

方位角が0〜360°で表される時は、下記のように算出できます。

from decimal import Decimal, ROUND_HALF_UP

def azimuth_to_direction(azimuth):
    # 16方位の方角を定義する
    directions = ['北', '北北東','北東', '東北東',
                  '東', '東南東', '南東', '南南東',
                  '南', '南南西', '南西', '西南西',
                  '西', '西北西', '北西', '北北西',
                  '北']
 
    # Decimalモジュールで小数第一位を四捨五入
    directions_index = Decimal(azimuth / 22.5).quantize(Decimal('1.'), rounding=ROUND_HALF_UP)

    return directions[int(directions_index)]

処理の手順は冒頭の記事のとおりですが、Pythonでは一般的な四捨五入をするのにround()関数ではなくDecimalモジュールを使用します。

より厳密に計算するのであれば、22.5で割る前にDecimal型に変換しておきます。
さらに、方位角がfloat型のときは、decimal型に変換するときに近似値となってしまうことがあるため、一度str型を挟みます。

value を float にする場合、二進浮動小数点数値が損失なく正確に等価な Decimal に変換されます。この変換はしばしば 53 桁以上の精度を要求します。例えば、 Decimal(float('1.1')) は Decimal('1.100000000000000088817841970012523233890533447265625') に変換されます。

# str型を挟んでDecimal型に変換
azimuth = Decimal(str(azimuth))

directions_index = (azimuth / Decimal('22.5')).quantize(Decimal('1.'), rounding=ROUND_HALF_UP)

return directions[int(directions_index)]

実際の例

pyprojのGeodクラスを用いて新千歳空港から見た那覇空港の方位角を求め^[1]、それを16方位に変換してみます。Geod.inv()メソッドの返り値に含まれる方位角は-180～180°表記のようなので、azimuth_to_direction()関数を修正し、0〜360°表記に変換してから処理するようにしました。

from decimal import Decimal, ROUND_HALF_UP

from pyproj import Geod

def azimuth_to_direction(azimuth):
    # 方位角を-180～180°から、0〜360°表記にする
    if azimuth < 0:
        azimuth += 360

    directions = ['北', '北北東','北東', '東北東',
                  '東', '東南東', '南東', '南南東',
                  '南', '南南西', '南西', '西南西',
                  '西', '西北西', '北西', '北北西',
                  '北']

    azimuth = Decimal(str(azimuth))
    directions_index = (azimuth / Decimal('22.5')).quantize(Decimal('1.'), rounding=ROUND_HALF_UP)

    return directions[int(directions_index)]

# GRS80楕円体を使用
g = Geod(ellps='GRS80')

# 新千歳空港の緯度経度（Wikipediaより）
cts_lat, cts_lon = 42.775, 141.692222
# 那覇空港の緯度経度（Wikipediaより）
oka_lat, oka_lon = 26.193333, 127.639722
 
f_azimuth, b_azimuth, distance = g.inv(cts_lon, cts_lat, oka_lon, oka_lat)

print(f_azimuth, azimuth_to_direction(f_azimuth))

コードを実行してみると、方位角は-140.67921357699439°、16方位では「南西」となりました。南西は-146.25〜-123.75°なので、正しく判定されていることが分かります。

脚注

pyproj 3.6.0で検証 ↩︎

イントロダクション

動作環境

コードと説明

実際の例

Discussion