【Pythonで地図描画】Basemapで実装されたコードをCartopyに書き換える

はじめに

可視化の際、地図を描画して見た目を良くするために、すでに開発が終了したBasemapを使い続けています。流石にそろそろ、後継のCartopyに乗り換えたいということで、実装してみました。本記事は、元のBasemapを用いたコードをどのように書き換えたのかを順に紹介します。本編に入る前に、Cartopyをインストールしておきます。（各自の環境に合わせてインストールして下さい。）

conda install -c conda-forge cartopy

最終的に生成される画像はこのような感じです。（緯度・経度線は３度ごとに引いています。）

元のBasemapを利用したコード

以下に元のコードを示します。predictには64x64の地域（メッシュ）毎の値（降水量など）が入ります。

Basemap

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

def draw_map(predict, lonRange=np.array([初期値,初期値]), latRange=np.array([初期値,初期値]), cmap="seismic", cbar_label="intensity", lonNum=64, latNum=64, norm=[0.0,1.0]):
    # 地域の描画
    m = Basemap(projection='merc',
                resolution='h',
                llcrnrlon=lonRange[0],
                llcrnrlat=latRange[0],
                urcrnrlon=lonRange[1],
                urcrnrlat=latRange[1])
    # 海岸線、陸地線などの色設定
    m.drawcoastlines(color='gray')
    m.drawcountries(color='gray')
    m.drawmapboundary(fill_color='#eeeeee')
    # 5度ごとに緯度線を描く
    m.drawparallels(np.arange(latRange[0], latRange[1], 10), labels = [1, 0, 0, 0], fontsize=10, color='white')
    # 5度ごとに経度線を描く
    m.drawmeridians(np.arange(lonRange[0], lonRange[1], 10), labels = [0, 0, 0, 1], fontsize=10, color='white')
    # 座標の描画
    X = np.linspace(lonRange[0], lonRange[1], lonNum)
    Y = np.linspace(latRange[1], latRange[0], latNum)
    lon, lat = np.meshgrid(X, Y)
    if len(predict):
        m.pcolormesh(lon, lat, predict, latlon=True, cmap=cmap, norm=norm)
    cbar = plt.colorbar()
    cbar.set_label(cbar_label)
    plt.show()

順に書き換える

地域の描画

まずは、地域の描画から始めます。Basemapでは以下のように設定していました。

Basemap

# 地域の描画
m = Basemap(projection='merc',
            resolution='h',
            llcrnrlon=lonRange[0],
            llcrnrlat=latRange[0],
            urcrnrlon=lonRange[1],
            urcrnrlat=latRange[1])

上記のコードは、左右上下の終端値を設定していました。
今回は、インタフェースを用いて下記のように書き換えました。

Cartopy

fig = plt.figure()
proj = ccrs.PlateCarree() # 正距円筒図法
ax = fig.add_subplot(1,1,1,projection=proj)
# 地域の描画
ax.set_extent((lonRange[0], lonRange[1], latRange[0], latRange[1]), proj)

海岸線の描画

元々のBasemapを用いたコードは海岸線の他に陸地線なども描画していましたが、特にしないため、今回は省略しました。

Basemap

# 海岸線、陸地線などの色設定
m.drawcoastlines(color='gray') # <-今回は海岸線のみ描画する
m.drawcountries(color='gray')
m.drawmapboundary(fill_color='#eeeeee')

Cartopy

# 海岸線を描画
ax.coastlines()

緯度・経度線を描く

ここは、上手い実装方法が無いかな〜と模索中です。。。

Basemap

# 5度ごとに緯度線を描く
m.drawparallels(np.arange(latRange[0], latRange[1], 10), labels = [1, 0, 0, 0], fontsize=10, color='white')
# 5度ごとに経度線を描く
m.drawmeridians(np.arange(lonRange[0], lonRange[1], 10), labels = [0, 0, 0, 1], fontsize=10, color='white')

Cartopy

# 5度ごとに緯度・経度線を描く
gl=ax.gridlines(crs=proj, draw_labels=True, linewidth=1, alpha=0.8)
gl.xlocator = tck.FixedLocator(np.arange(lonRange[0], lonRange[1], 5)) # 5度ごとに経度線を描く
gl.ylocator = tck.FixedLocator(np.arange(latRange[0], latRange[1], 5)) # 5度ごとに緯度線を描く
gl.top_labels = False # 上部のラベルを非表示
gl.right_labels = False # 右側のラベルを非表示

うーん、Basemapの方が見やすい。

座標の描画

座標の描画については、ほとんど変更はありません。
一部、pcolormeshに差異があります。

Basemap

# 座標の描画
X = np.linspace(lonRange[0], lonRange[1], lonNum)
Y = np.linspace(latRange[1], latRange[0], latNum)
lon, lat = np.meshgrid(X, Y)
if len(predict):
    m.pcolormesh(lon, lat, predict, latlon=True, cmap=cmap, norm=norm)

Cartopy

# 座標の描画
X = np.linspace(lonRange[0], lonRange[1], lonNum)
Y = np.linspace(latRange[1], latRange[0], latNum)
lon, lat = np.meshgrid(X, Y)
if len(predict):
    plt.pcolormesh(lon, lat, predict, cmap=cmap, vmin=norm[0], vmax=norm[1])

。。。上手くNormalizationが設定できなかったので、とりあえず上記のようにしています。
あと、ax.pcolormeshにすると、この後のカラーバーを設定できませんでした。

最終的なコード

最終的には下記のようなコードになりました。思いの外、すぐに書き換えられました。いくつか要検討な部分がありますが、とりあえずこれにて作業終了ということにします。このように実装した方が良いなどのアドバイスがありましたら、コメント頂けると幸いです！

Cartopy

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as tck
import cartopy.crs as ccrs

def draw_map(predict, lonRange=np.array([初期値,初期値]), latRange=np.array([初期値,初期値]), cmap="seismic", cbar_label="intensity", lonNum=64, latNum=64, norm=[0.0,1.0]):
    fig = plt.figure()
    proj = ccrs.PlateCarree() # 正距円筒図法
    ax = fig.add_subplot(1,1,1,projection=proj)
    # 地域の描画
    ax.set_extent((lonRange[0], lonRange[1], latRange[0], latRange[1]), proj)
    # 海岸線を描画
    ax.coastlines()
    # 5度ごとに緯度・経度線を描く
    gl=ax.gridlines(crs=proj, draw_labels=True, linewidth=1, alpha=0.8)
    gl.xlocator = tck.FixedLocator(np.arange(lonRange[0], lonRange[1], 5)) # 5度ごとに経度線を描く
    gl.ylocator = tck.FixedLocator(np.arange(latRange[0], latRange[1], 5)) # 5度ごとに緯度線を描く
    gl.top_labels = False # 上部のラベルを非表示
    gl.right_labels = False # 右側のラベルを非表示

    # 座標の描画
    X = np.linspace(lonRange[0], lonRange[1], lonNum)
    Y = np.linspace(latRange[1], latRange[0], latNum)
    lon, lat = np.meshgrid(X, Y)
    if len(predict):
        plt.pcolormesh(lon, lat, predict, cmap=cmap, vmin=norm[0], vmax=norm[1])
    cbar = plt.colorbar()
    cbar.set_label(cbar_label)

    plt.show()

参考記事

Basemapの実装時に参考にしたブログ：pythonでの地図表示
Cartopyインストール時に参考にしたサイト：Cartopyの公式ページ