【Science】水平線までどのくらいの距離？Part2ーマップで表示

こんにちは！
前回に引き続き、Cartopyでマップに表示したいと思います

[前回の記事]

水平線までどのくらいの距離なのか？

runtascience.hatenablog.com

地図
コード

今回はCartopyを使ってどこまで見られるのか表示してみました！

都道府県などの境を表示するためにGADMのシェープファイルを利用させていただきました。

[ シェープファイル利用(Japan⇒Shapefile)]

gadm.org

地図

大気スケール

まずは大気スケールです。

日本列島の端から端まではだいたい2000 kmです

したがって、半径が1000 kmだとだいたい日本列島全体が見れますね

飛行機の高度の10 kmくらいの時はいちばん小さい円です。

それでも関東地方全体は軽く眺められますね。

f:id:RuntaScience:20200620103556p:plain

スケール	H		D
スケール	名前	m	km
大気	熱圏	10000	1133.23
	中間層	8000	1012.8
	成層圏	5000	799.75
	対流圏	1000	357.099

山スケール

次に山スケールです。

今回は富士山とエベレストのみにします。(他はそこまで高さが変わらないので)

赤の点は富士山頂です。

富士山とエベレストで見られる範囲を比べると100 km程しか変わらないのでそこまで違いが見られませんね。それにしても、山に登るだけでここまで見れるのはすごいですよね。

f:id:RuntaScience:20200620103549p:plain

建物スケール

次に建物スケールです。

赤い点は東京スカイツリーです。(634 m)

世界1位のブルジュハリファとスカイツリーを比較してもあまり見える距離は変わらないんですね。意外です。

しかし、東京タワーと比べると差が大きくなりますね。

都内の大きな高層マンションでは水平距離で24区見れそうですね、障害物がなければですが(笑

f:id:RuntaScience:20200620103539p:plain

スケール	H		D
スケール	名前	m	km
建物	ブルジュハリファ	828	102.7
	東京スカイツリー	634	89.9
	東京タワー	333	65.1
	福岡タワー	234	54.6
	40階	150	43.7
	20階	60	27.7
	10階	30	19.6
	5階	12	12.4
	3階	9	10.7
	2階	6	8.7

人スケール

最後に人スケールです。

平均身長あたりだと150ー170 cmなので見られる距離は4キロほどです。

東京ディズニーリゾートのある舞浜からです。

実際東京湾に立つと、羽田空港から千葉県が見られますし、千葉県から東京や川崎が見られます。

それは、水平線ではなく建物の高さがあるためや地上の高低差ならびに大気の屈折があるためと考えられます。

f:id:RuntaScience:20200620103532p:plain

スケール	H		D
	名前	m	km

身長	10cm	0.1	1.1
	30cm	0.3	2
	50cm	0.5	2.5
	70cm	0.7	3
	90cm	0.9	3.4
	110cm	1.1	0.7
	130cm	1.3	4.1
	150cm	1.5	4.4
	170cm	1.7	4.7
	190cm	1.9	4.9
	210cm	2.1	5.2
	230cm	2.3	5.4
	250cm	2.5	5.6
	270cm	2.7	5.9

コード

まずは描写に必要なモジュールをインポートします

import matplotlib.pyplot as plt       
import cartopy.crs as ccrs               
import cartopy.feature as cfeature 
#shape file
import cartopy.io.shapereader as shpreader

次に、地球状のどこかかの点から円を描くので

半径(r km)と中心の緯度経度から、円の緯度経度座標に変換します

# 円
def circle(r,S,T):
    #r:radius(km), (S,T):centerDegree
    s,t = S*111, T*111
    points = 1000  
    theta = np.linspace(0, 2 * np.pi, points)
    x = r * np.sin(theta) + s
    y = r * np.cos(theta) + t
    
    #x,y=>Degree
    X, Y = x/111, y/111
    return (X,Y)

Keys

「km⇄緯度経度」は赤道半径を利用した値で変換しました

km = Degree ✕ 111

それでは描写してみます

def main():
    fig = plt.figure(figsize=(5,5))
    plt.rcParams["font.size"] = 18
    # TokyoStation
    ts_lon = 139.767125
    ts_lat = 35.681236
    
    ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
    ax.set_extent([ts_lon-1, ts_lon+1, ts_lat-1, ts_lat+1], crs=ccrs.PlateCarree())
    ax.set_title("Title")

    ax.add_feature(cfeature.BORDERS, linestyle=":")
    #都道府県境と市町村境
    fname1 = "data/gadm36_JPN_shp/gadm36_JPN_1.shp" 
    fname2 = "data/gadm36_JPN_shp/gadm36_JPN_2.shp"
    adm1_shapes = list(shpreader.Reader(fname1).geometries())
    adm2_shapes = list(shpreader.Reader(fname2).geometries())
    
    #データのプロット
    ax.add_geometries(adm1_shapes, ccrs.PlateCarree(),
                  edgecolor="k",facecolor="w",alpha=0.5) 
    ax.add_geometries(adm2_shapes, ccrs.PlateCarree(),
                  edgecolor="gray", facecolor="w",  alpha=0.5) 
    
    #TokyoStation
    
    ax.scatter(ts_lon,ts_lat,s=20,c="r",alpha=1, zorder=10)
    c = circle(100, ts_lon,ts_lat)
    ax.plot(c[0],c[1],c="b",alpha=1, zorder=10)
    
    plt.show()

if __name__ == "__main__":
    main()

参考コード

qiita.com

geopandas.readthedocs.io

それでは🌏

runtascience.hatenablog.com