Python AI 画像

cartorpyで地図を描く(地理データを可視化する)。

同様のライブラリにbasemapがありますが、2020年でサポート終了するようですのでこちらを使用しました。
データの可視化にmatplotlibも使用します。 Cartopyは、イギリス気象局(MetOffice)によりオープンソースライブラリScitoolsの一部として開発されている地図描画用のPythonライブラリです。
conda install cartopyでインストールできます。

In [103]:
# 地図を描く basemapは2020でサポート終了

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as tick
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter

地図データは、海岸線、湖川、湖、陸地、海洋ごとに、10m、50m、110mの解像度があります。
初めて利用するときに.localにダウンロードされます。

In [104]:
# 10m-coastline, 110m-coastline 地図データは./.localに保存された
land_50m=cfeature.NaturalEarthFeature('physical','land','50m',edgecolor='face',facecolor=cfeature.COLORS['land'])
#land_50m=cfeature.NaturalEarthFeature('cultual','land','50m',edgecolor='face',facecolor=cfeature.COLORS['land'])

図法には、メルカトル図法、ランベルト正角円錐図法、コンピュータでよく使われる正距円筒図法(PlateCaree)があります。

In [105]:
plt.figure()
ax = plt.axes(projection = ccrs.Mercator( # メルカトール図法です
    central_longitude = 0.0,
    min_latitude = -80.0, # メルカトール図法の場合、90は設定できない
    max_latitude = 80.0,
    globe = None,
    latitude_true_scale = 0.0))
ax.add_feature(cfeature.LAND)
#ax.add_feature(cfeature.COASTLINE)
ax.coastlines(lw = 0.5) # line width
ax.gridlines(linestyle = '-', color = 'grey')
ax.set_title('Mercator')
plt.show()
In [106]:
plt.figure(figsize = (10, 5), facecolor = 'w')
ax = plt.axes(projection = ccrs.PlateCarree(central_longitude = 180)) # 正距円筒図法、標準緯線0度

ax.stock_img() # 抑揚のある地図
    
dlon, dlat = 60, 30
xticks = np.arange(0, 360.1, dlon)
yticks = np.arange(-90, 90.1, dlat)

#グリッド線を引く
gl = ax.gridlines(crs = ccrs.PlateCarree(), draw_labels = False, linewidth = 1, alpha = 0.8)
gl.xlocator = tick.FixedLocator(xticks)    
gl.ylocator = tick.FixedLocator(yticks)

ax.set_xticks(xticks, crs = ccrs.PlateCarree())
ax.set_yticks(yticks, crs = ccrs.PlateCarree())

#目盛の表示形式を度数表記にする    
latfmt = LatitudeFormatter()
lonfmt = LongitudeFormatter(zero_direction_label = True)
ax.xaxis.set_major_formatter(lonfmt)
ax.yaxis.set_major_formatter(latfmt)
#目盛のサイズを指定
ax.axes.tick_params(labelsize = 12)

ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE, linewidth = 0.6)
#ax.coastlines(lw = 0.5)
ax.add_feature(cfeature.LAKES, alpha = 0.5)
ax.add_feature(cfeature.RIVERS)
ax.gridlines(linestyle = '-', color = 'grey')
ax.set_title('PlateCaree')
plt.show()

特定の地域(ここでは日本)だけを拡大することもできます。

In [107]:
plt.figure()
ax = plt.axes(projection = ccrs.PlateCarree()) # 正距円筒図法、標準緯線0度
ax.coastlines(resolution = '50m')

# 描画範囲
west = 120
east = 150
south = 20
north = 50
ax.set_extent([west, east, south, north], ccrs.PlateCarree())

#ax.add_feature(cfeature.LAND)
land_10m  = cfeature.NaturalEarthFeature('physical', 
                                         'land', 
                                         '10m',
                                         edgecolor = 'gray', # none
                                         facecolor = cfeature.COLORS['land'])
ax.add_feature(land_10m, linewidth = 0.3) # 10m resolution
#ax.add_feature(cfeature.OCEAN)
ocean_10m  = cfeature.NaturalEarthFeature('physical', 
                                         'ocean', 
                                         '10m',
                                         edgecolor = 'gray', # none
#                                         facecolor = 'aqua')
                                         facecolor = cfeature.COLORS['water'])
ax.add_feature(ocean_10m, linewidth = 0.3) # 10m resolution
#ax.add_feature(cfeature.COASTLINE, linewidth = 0.3)
#ax.coastlines(lw = 0.5) # line width
ax.add_feature(cfeature.LAKES, alpha = 0.5)
ax.add_feature(cfeature.RIVERS)

dlon, dlat = 10, 5
xticks = np.arange(west, east + 0.1, dlon)
yticks = np.arange(south, north + 0.1, dlat)    #目盛を描く緯度経度の値を設定
ax.set_xticks(xticks, crs = ccrs.PlateCarree())
ax.set_yticks(yticks, crs = ccrs.PlateCarree())
#目盛の表示形式を度数表記にする    
latfmt = LatitudeFormatter()
lonfmt = LongitudeFormatter(zero_direction_label = True)
ax.xaxis.set_major_formatter(lonfmt)
ax.yaxis.set_major_formatter(latfmt)
#目盛のサイズを指定
ax.axes.tick_params(labelsize = 12)
#グリッド線を引く
gl = ax.gridlines(crs = ccrs.PlateCarree(), draw_labels = False, linewidth = 1, color = 'grey', alpha = 0.8)
gl.xlocator = tick.FixedLocator(xticks)    
gl.ylocator = tick.FixedLocator(yticks)
#ax.gridlines(linestyle = '-', color = 'grey', linewidth = 0.5)

plt.show()

川や都道府県情報、特定点マークなども記載できます

In [108]:
fig = plt.figure(figsize=(8,4))

# physical category
land_50m  = cfeature.NaturalEarthFeature('physical', 'land', '50m', 
            edgecolor='face',                   # same color with facecolor
            facecolor=cfeature.COLORS['land'])  # use predefiend color of cartopy
ocean_50m = cfeature.NaturalEarthFeature('physical', 'ocean', '50m', 
            edgecolor='face',                   # same color with facecolor
            facecolor=cfeature.COLORS['water']) # use predefiend color of cartopy
lakes_50m = cfeature.NaturalEarthFeature('physical', 'lakes', '50m', 
            edgecolor='face',                   # same color with facecolor
            facecolor=cfeature.COLORS['water']) # use predefiend color of cartopy
river_50m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '50m',
            edgecolor=cfeature.COLORS['water'], # use predefiend color of cartopy
            facecolor='none')                   # no filled color
    
# cultural category
countries_50m  = cfeature.NaturalEarthFeature('cultural', 'admin_0_countries', '50m', 
                 edgecolor='gray',
                 facecolor='none')  # no filled color
states_50m  = cfeature.NaturalEarthFeature('cultural', 'admin_1_states_provinces_lines', '50m', 
              edgecolor='gray',
              facecolor='none')  # no filled color
states_10m  = cfeature.NaturalEarthFeature('cultural', 'admin_1_states_provinces_lines', '10m', 
              edgecolor='gray',
              facecolor='none')  # no filled color
    
ax = plt.axes(projection=ccrs.PlateCarree()) # 正距円筒図法、標準緯線0度
#ax.add_feature(cfeature.LAND)
ax.set_extent([130, 145, 30, 40], ccrs.PlateCarree())
ax.coastlines(resolution = '10m', color = 'gray')
ocean_10m  = cfeature.NaturalEarthFeature('physical', 
                                         'ocean', 
                                         '10m',
                                         edgecolor = 'gray',
                                         facecolor=cfeature.COLORS['water'])
ax.add_feature(ocean_10m, linewidth = 0.3) # 10m resolution
#ax.add_feature(cfeature.OCEAN)
ax.add_feature(states_10m, color = 'y', linewidth = 0.6)
#ax.add_feature(river_50m)
ax.gridlines(linestyle = '--', color = 'grey', linewidth = 0.5)
ax.set_title('Japan')

lat = 36
lon = 140
ax.plot(lon, lat, 'm.', markersize = 10, markerfacecolor = 'r', markeredgecolor = 'b')
ax.text(lon+0.1, lat+0.11, 'tsukuba', fontsize = 12)

plt.show()

半球。地球の無限遠方に光源をおいて平面に投影する正射方位図法。

In [109]:
plt.figure(figsize = (10, 5), facecolor = 'w')
ax = plt.axes(projection = ccrs.Orthographic(central_longitude = 136, central_latitude = 35)) # 正射方位図法

ax.stock_img() # 抑揚のある地図
    
#グリッド線を引く
gl = ax.gridlines(crs = ccrs.Orthographic(), draw_labels = False, linewidth = 1, alpha = 0.8)

ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE, linewidth = 0.4)
ax.add_feature(cfeature.LAKES, alpha = 0.5)
ax.add_feature(cfeature.RIVERS)
ax.coastlines(lw = 0.5)
ax.gridlines(linestyle = '-', color = 'grey')
ax.set_title('Orthographic')
plt.show()
In [110]:
import datetime
from cartopy.feature.nightshade import Nightshade

世界地図上に昼夜を描く

In [111]:
fig = plt.figure(figsize = (10, 5))
ax = fig.add_subplot(1, 1, 1, projection = ccrs.PlateCarree()) # (正距円筒図法 central_longitude=0.0, globe=None)
#ax = fig.add_subplot(1, 1, 1, projection = ccrs.Orthographic()) # (正投影図法 central_long/lat=0.0, globe=None)
ax.gridlines(linestyle = '-', color = 'grey') # 経緯度線表示

date = datetime.datetime(2019, 8, 19, 3) # 表示する日時を指定

ax.set_title('Night time shading for {}'.format(date)) # タイトルに日時も表示
ax.stock_img() # 抑揚のある地図
ax.add_feature(Nightshade(date, alpha = 0.2))

plt.show()
In [ ]: