Python + Jupyter で地図上に円を描く

Python + Jupyter で地図上に絵を描くには folium というライブラリを使います。例えばこんなふうにすると稚内から 1000 km の地点に円を描く事が出来ます。

import folium

WAKKANAI = [45.4171, 141.6761] # 稚内 (緯度、経度)

m = folium.Map(location=WAKKANAI, zoom_start=4)
folium.Circle(location=WAKKANAI, radius=1_000_000).add_to(m)
m

ただ、これは見て楽しむ以外あまり役に立たないでしょう。folium は円の描画を Python ではなくブラウザ側で行うので、例えばせっかく求めた座標を Python を使って「稚内と大阪の両方から 1000 km 離れた点を求める」ような問題に使えないからです。

Python で位置情報を扱うためには shapely というライブラリを使います。shapely には円を描画する機能が無いので、「稚内の位置からのオフセットを求める」事で代用してみます。

shapely には地図上の距離を指定する機能も無いので、1000 km を角度に換算します。folium での座標の並びが「緯度・経度」なのに対して、shapely では「経度・緯度」なので逆にしないといけないのがちょっと面倒です。

import math
from shapely.geometry import Point, mapping

EQUATORIAL_RADIUS = 6_378_137 # 地球の半径
distance_degree = 1_000_000 / (EQUATORIAL_RADIUS * 2 * math.pi) * 360 # 1000 km を角度に変換。

point_wakkanai = Point(reversed(WAKKANAI))
circle_wakkanai = point_wakkanai.buffer(distance_degree) # 経度,緯度の順にするため reverse

m = folium.Map(location=WAKKANAI, zoom_start=4)
folium.GeoJson(circle_wakkanai).add_to(m)
m

なんか変ですね。。。円ではなく多角形なのは仕方ないとして、タマゴ型になってしまいました。

folium は座標を緯度経度で示しますが、shapely は直交座標とみなして図形を操作します。緯度経度というのは球面を表すための物なので、shapely で求めた座標を地図上にそのまま描いてしまうと歪むのです。緯度が高い場所の経度1度は、赤道上の経度1度よりも短くなります。

Python には、座標変換を行うための pyproj というライブラリがあります。pyproj で座標変換を行うためには、変換前と変換後の座標系が必要です。pyproj では座標系を "proj-string" https://proj.org/usage/quickstart.html という物で示します。

変換前の座標系はややこしいですが、中心からの距離と角度が正しくなる正距方位図法を使いましょう。稚内を中心とする正距方位図法を表す "proj-string" は次のようになります。

aeqd_wakkanai = "+proj=aeqd +lat_0=45.4171 +lon_0=141.6761 +x_0=0 +y_0=0"

変換後の座標系は、電子地図や GPS データ使われる WGS84 を表す proj-stinrg は次のようになります。

wgs84 = "+proj=longlat +datum=WGS84"

これを使って正距方位図法に描いた円を WGS84 に変換すると、WGS84 の地図上に円が描けます。

import pyproj
from shapely.ops import transform

def circle(point, radius):

    aeqd = f'+proj=aeqd +lat_0={point.y} +lon_0={point.x} +x_0=0 +y_0=0' # point 中心の正距方位図法
    wgs84 = "+proj=longlat +datum=WGS84" # WGS84
    project = pyproj.Transformer.from_crs(aeqd, wgs84).transform # 正距方位図法からWGS84への変換。

    circle_aeqd = point.buffer(radius)
    return transform(project, circle_aeqd)

circle_wakkanai = circle(point_wakkanai, 1_000_000)

m = folium.Map(location=WAKKANAI, zoom_start=4)
folium.GeoJson(circle_wakkanai).add_to(m)
m

一応、稚内と大阪の両方から 1000km 離れた場所を探してみます。

OSAKA = [34.6869, 135.5259]

point_osaka = Point(reversed(OSAKA))
circle_osaka = circle(point_osaka, 1_000_000)

intersection = circle_wakkanai.boundary.intersection(circle_osaka.boundary) # 二つの円の交点

m = folium.Map(location=WAKKANAI, zoom_start=4)
folium.GeoJson(circle_wakkanai).add_to(m)
folium.GeoJson(circle_osaka).add_to(m)
folium.GeoJson(intersection).add_to(m)
m

検算として距離を測ります。

詳細は私もよく知りませんが、pyproj の機能で距離も測れます。

geod = pyproj.Geod(ellps='WGS84')
print("大阪からの距離: ", 
    geod.inv(point_osaka.x, point_osaka.y, intersection[0].x, intersection[0].y)[2])
print("稚内からの距離: ", 
    geod.inv(point_wakkanai.x, point_wakkanai.y, intersection[0].x, intersection[0].y)[2])

大阪からの距離:  999640.1189235467
稚内からの距離:  998819.0687679343

ちょっと違いますが、こんなもんでしょうか? 誤差がどこで出るのか分かったらまた書きます。

参考

• folium: https://leafletjs.com/
• shapely: https://shapely.readthedocs.io/en/stable/
• pyproj: https://pyproj4.github.io/pyproj/stable/index.html
• proj-string: https://proj.org/usage/quickstart.html
• Creating buffer circle x kilometers from point using Python?
• 正距方位図法: https://ja.wikipedia.org/wiki/%E6%AD%A3%E8%B7%9D%E6%96%B9%E4%BD%8D%E5%9B%B3%E6%B3%95