正距円筒図法 (PlateCaree)で線や点をプロットするときはplt.plot()などを使えば、プロットすることができます。普通のグラフと同じ感覚で描けるので楽ですね。しかし、正射方位図法(Orthographic)を使う場合は、もう一工夫する必要があります。
今回は「世界の都市総合ランキング」(Global Power City Index, GPCI)
の2019年度版から、
http://www.mori-m-foundation.or.jp/english/ius2/gpci2/
交通・アクセスの分野でtop5にランクインした市を線でつなげてみました。出来上がりはこんな感じです。
スポンサーリンク
Pythonコード
以下のstackoverflowをベースにしています。ここに組み合わせ(順番を考慮しない順列)を使って地球に線を引いていきます。
#!/usr/bin/env python3
#必要なものを宣言
import matplotlib.pyplot as plt
from matplotlib.patheffects import withStroke
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np
import shapely.geometry as sgeom
import itertools
#保存する場所を指定
DIR1 = '/home/User/city/figs'
#解像度の指定
set_dpi = 100
#対象となる市のリスト
city_list = ['Paris', 'London', 'New York', 'Shanghai', 'Frankfurt']
#対象の市の緯度経度
point = np.array([[2.21, 48.51], [0.07, 51.30],[-74.00, 40.42], [121.29, 31.10], [8.40, 50.06]])
nsite = len(point)
#組み合わせをするため、一度listに戻す
line_list = point.tolist()
#5つの市から2つの市を選ぶときの組み合わせをlist表示
c_list = list(itertools.combinations(line_list, 2))
#組み合わせのlistを２列のarrayに変換
line = np.array((c_list), np.float64).reshape(-1, 2)
#sgeom.LineStringの引数はlistしかとらないので、
#再度listに戻す
line2 = line.tolist()
lines = sgeom.LineString(line2)
#ここから作図
plt.figure(figsize=(12,9))
#ここでclassを指定しておかないと、projected_lineが使えない
class HighResPC(ccrs.PlateCarree):
pass
#ccrs.Orthographicに対応できるShapely geometryをつくる
projected_line = HighResPC().project_geometry(lines, ccrs.Geodetic())
#projectionと中心とする緯度経度の指定
ax = plt.axes(projection=ccrs.Orthographic(central_latitude=82, central_longitude=80))
#陸地、海、海岸線
ax.add_feature(cfeature.LAND,facecolor='#999999',
linewidth=0.1,edgecolor='#999999',zorder=1)
ax.add_feature(cfeature.OCEAN,zorder=0,edgecolor='none', alpha=0.8)
ax.coastlines(lw=0.5)
#それぞれの市を結ぶ線を追加
ax.add_geometries(
[projected_line], HighResPC(),
edgecolor='#EFEA3A', facecolor='none', alpha=0.5, lw=0.8)
for isite in range(nsite):
#市に丸い点を追加
ax.scatter(point[isite,0], point[isite,1], 50, marker='o',
c='#EFEA3A', transform=ccrs.Geodetic(), zorder=3)
#市の名称のedgecolor（縁取りの色)を指定
effects=[withStroke(linewidth=1.0,
foreground='k')]
#位置の微調整
if isite == 0: #パリ
x_diff = -4
y_diff = -14
elif isite == 4: #フランクフルト
x_diff = 4
y_diff = -8
else: #その他
x_diff = 4
y_diff = 4
#市の名称をプロット
ax.text(point[isite,0]+x_diff, point[isite,1]+y_diff,city_list[isite],
transform=ccrs.Geodetic(), fontsize=18, horizontalalignment='left',
color='#EFEA3A', path_effects=effects)
#海底の形状を表示
ax.stock_img()
#これがないと球体にならないので必ずいれる！
ax.set_global()
#自動調節
plt.tight_layout()
#pngで保存
fileo1 = DIR1 +'access_globalcity'+'.png'
plt.savefig(fileo1,dpi=set_dpi, transparent=True)
組み合わせを使う理由
コードを見た方の中には、なぜ組み合わせを使っているのだろうか？リストや配列をなんどもいじる必要があるのだろうか？と疑問に思うかもしれません。（もっと簡単にする方法があれば、教えてほしいですが・・・！）
ということで、組み合わせとax.add_geometries(・・・)の関係を以下の画像にまとめました。簡単にいえば直線を引く順番をきめるために組み合わせを使っているということです。
スポンサーリンク
使用したリスト・配列を解説
リスト化したり配列に戻したりする工程を２回行っているので、コードだけを見ると少し混乱してしまうかもしれませんね。
[projected_line]のもとになるline2が出来上がる工程まで、使用したリスト・配列を順を追って解説します。
points (np.array)
もとの配列です。プリントするとdtypeがfloat64で５×２の配列ができます。
points = [[ 2.2100e+00 4.8510e+01]
[ 7.0000e-02 5.1300e+01]
[-7.4000e+01 4.0420e+01]
[ 1.2129e+02 3.1100e+01]
[ 8.4000e+00 5.0060e+01]]
c_list (list)
組み合わせは５つの市(n=5)から２つ(r=2)を選ぶのでnCrを計算すると１０通りになります。
リストのサイズもちゃんと１０になってますね。
c_list = [([2.21, 48.51], [0.07, 51.3]), ([2.21, 48.51], [-74.0, 40.42]), ([2.21, 48.51], [121.29, 31.1]),
([2.21, 48.51], [8.4, 50.06]), ([0.07, 51.3], [-74.0, 40.42]), ([0.07, 51.3], [121.29, 31.1]),
([0.07, 51.3],[8.4, 50.06]), ([-74.0, 40.42], [121.29, 31.1]), ([-74.0, 40.42], [8.4, 50.06]),
([121.29, 31.1], [8.4, 50.06])]
len(c_list) = 10
line (np.array)
c_listの括弧をはずしたいので、２列にリシェイプしたものがlineです。
１０個あった括弧それぞれには２つの市の緯度経度が入っているので、１０×２で２０行の配列ができています。
line = [[ 2.2100e+00 4.8510e+01]
[ 7.0000e-02 5.1300e+01]
[ 2.2100e+00 4.8510e+01]
[-7.4000e+01 4.0420e+01]
[ 2.2100e+00 4.8510e+01]
[ 1.2129e+02 3.1100e+01]
[ 2.2100e+00 4.8510e+01]
[ 8.4000e+00 5.0060e+01]
[ 7.0000e-02 5.1300e+01]
[-7.4000e+01 4.0420e+01]
[ 7.0000e-02 5.1300e+01]
[ 1.2129e+02 3.1100e+01]
[ 7.0000e-02 5.1300e+01]
[ 8.4000e+00 5.0060e+01]
[-7.4000e+01 4.0420e+01]
[ 1.2129e+02 3.1100e+01]
[-7.4000e+01 4.0420e+01]
[ 8.4000e+00 5.0060e+01]
[ 1.2129e+02 3.1100e+01]
[ 8.4000e+00 5.0060e+01]]
line.shape = (20, 2)
line2 (list)
lineをリスト化したものがline2です。c_listのような括弧が外れていることがわかります。
これでsgeom.LineStringをすれば、地球に線を張り付ける準備ができます！
line2 = [[2.21, 48.51], [0.07, 51.3], [2.21, 48.51], [-74.0, 40.42], [2.21, 48.51], [121.29, 31.1], [2.21, 48.51],
[8.4, 50.06], [0.07, 51.3], [-74.0, 40.42], [0.07, 51.3], [121.29, 31.1], [0.07, 51.3], [8.4, 50.06],
[-74.0, 40.42], [121.29, 31.1], [-74.0, 40.42], [8.4, 50.06], [121.29, 31.1], [8.4, 50.06]]