Pythonで世界地図-2

Pythonで世界地図
https://qiita.com/ty21ky/items/c7a93b83cbbfa8a82a17
の続き

#####地球を見る位置（東経、北緯、高度）を指定して地図を作成する。

```python:map.py
#!/usr/bin/python
# coding: UTF-8

from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
import sys

def get_input(prompt):
    if sys.hexversion > 0x03000000:
        return input(prompt)
    else:
        return raw_input(prompt)

# 東経・北緯・地球を見る高度を入力
lon_0 = float(get_input('東経を入力。西経の場合はマイナスを付ける (東経):'))
lat_0 = float(get_input('北緯を入力。南緯の場合はマイナスを付ける (北緯):'))
h = float(get_input('カメラの高度 (m) (h):'))

if h < 150000:  #経線、緯線のピッチ
    pitch = 1
elif h < 1000000:
    pitch = 5
else:
    pitch = 10

# 地図表示
m = Basemap(projection='nsper',lon_0=lon_0,lat_0=lat_0,satellite_height=h,resolution='i')
m.drawcoastlines()
m.drawmapboundary(fill_color='aqua')
m.fillcontinents(color='coral',lake_color='aqua')
m.drawcountries()
m.drawstates()

m.drawmeridians(np.arange(0, 360, pitch)) #経線　表示　pitch度毎
m.drawparallels(np.arange(-90, 90, pitch)) #緯線　表示　pitch度毎
m.drawmapboundary(fill_color='aqua')
plt.title('east longitude=%s, north latitude=%s \n Line spacing=%s,Height=%g m' %(lon_0,lat_0,pitch,h),fontsize=15) #タイトル

plt.show()
```

気象衛星「ひまわり５号」から見た地球

![気象衛星ひまわり５号から見た地球.png](https://qiita-image-store.s3.amazonaws.com/0/238800/66ca3d65-1697-7064-12b6-24194ce9f5a1.png)

10000mの上空から見た地球

![10000m上空の飛行機から見た地球.png](https://qiita-image-store.s3.amazonaws.com/0/238800/7c077610-3f5f-1e8e-6129-71e8e40fbf0f.png)

#####大圏コース

```python:GreatCircleRoute.py
#!/usr/bin/python
# coding: UTF-8

from mpl_toolkits.basemap import Basemap, pyproj
import matplotlib.pyplot as plt
import numpy as np
import sys

def get_input(prompt):
    if sys.hexversion > 0x03000000:
        return input(prompt)
    else:
        return raw_input(prompt)

#出発地と目的地の緯度、経度を入力
lon1 = float(get_input('出発地の東経を入力。西経の場合はマイナスを付ける (東経):'))
lat1 = float(get_input('出発地の北緯を入力。南緯の場合はマイナスを付ける (北緯):'))
lon2 = float(get_input('目的地の東経を入力。西経の場合はマイナスを付ける (東経):'))
lat2 = float(get_input('目的地の北緯を入力。南緯の場合はマイナスを付ける (北緯):'))
 
# 地図表示
map = Basemap(projection='robin', lat_0=0, lon_0=150,resolution='l', area_thresh=1000.0)

map.drawgreatcircle(lon1,lat1,lon2,lat2,linewidth=2,color='b')  #大圏コース

# ポイント間の距離計算
gc = pyproj.Geod(a=map.rmajor,b=map.rminor)
az1, az2, dist12 = gc.inv(lon1,lat1,lon2,lat2)
 
map.drawcoastlines()
map.drawcountries()
map.fillcontinents(color='coral')
map.drawmapboundary()
 
map.drawmeridians(np.arange(0, 360, 10))
map.drawparallels(np.arange(-90, 90, 10))

plt.title('Great Circle from lon1=%s, lat1=%s \nto lon2=%s, lat2=%s \nDistance between two points=%s km' %(lon1,lat1,lon2,lat2,dist12/1000),fontsize=15)
 
plt.show()
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/238800/453a6875-82bd-84b0-a55c-47c1365e674b.png)

#####スケールの表示
Google Mapと比べたら、非常に誤差が大きいような気がします。
2点間距離を測定する方が確実で、これはあくまでも目安程度。

```python:scale.py
#!/usr/bin/python
# coding: UTF-8

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

# 緯度経度で範囲を指定する
north = 46.
south = 30.
east = 147.
west = 128.

# 地図の表示
m = Basemap(projection='merc',llcrnrlat=south,urcrnrlat=north, resolution='l', area_thresh=0.0,llcrnrlon=west,urcrnrlon=east)
# 海岸線を引く
m.drawcoastlines()

# 陸地を茶色に, 湖を水色に
m.fillcontinents(color='coral', lake_color='aqua')

# 海を水色に
m.drawlsmask(ocean_color='aqua')

m.drawmeridians(np.arange(0, 360, 10)) #緯線　表示　例.10度毎
m.drawparallels(np.arange(-90, 90, 10)) #経線　表示　例.10度毎

m.drawmapscale(140,32, 140,32, 1000, barstyle='fancy', units='km') #scale bar

plt.show()
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/238800/7fa08fdc-303a-81c1-55b5-176911cb5a9f.png)

#####昼と夜の領域

```python:daynight.py
#!/usr/bin/python
# coding: UTF-8

import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from pytz import timezone
from datetime import datetime

# 昼と夜の領域

map = Basemap(projection='robin',lon_0=150)
# 海岸線
map.drawcoastlines()
map.drawmeridians(np.arange(0, 360, 10)) #緯線　表示　例.10度毎
map.drawparallels(np.arange(-90, 90, 10)) #経線　表示　例.10度毎
# 色の指定
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')

# 夜間領域を陰影付けし、地図が表示されるようにアルファ透明度を設定します。
# UTCで現在の時刻を使用します。
date = datetime.utcnow()
CS=map.nightshade(date)

# JSTで現在の時刻を使用します(変換の仕方がよくわからない)
date1 = timezone('Asia/Tokyo').localize(datetime.now())
plt.title('Day/Night Map for %s (JST)' % date1.strftime("%d %b %Y %H:%M:%S"))
plt.show()
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/238800/e5f39b0c-cf1e-f591-3e22-d11a53d661cc.png)