cartopyで地図を描くときに軸ラベルとグリッドを整える

cartopyでグリッド線や軸目盛りの体裁を整えるための自作関数を作成します

はじめに

matplotlibで地図を描く際に使用されてきたBasemapは開発が終了しcartopyへと移行することが決定しているようです。基本的にはcartopyの方がコードが簡潔に書けるのですが、グリッド線の調整を行うのが少し面倒だったので、本記事ではcartopyの公式ドキュメントのギャラリーを参考にして調整を行ってくれる関数を作成し、簡単に描画できるようにします。

環境

描画はJupyterNotebook上で行っています。
matplotlib :3.0.0, 2.2.2
cartopy :0.16.0
numpy :1.15.1, 1.14.3

ひとまず軸まわりを調整して地図を描いてみる

以後コードはmatplotlibのオブジェクト指向インターフェイスで記述していきます。

地図にグリッド線を引く

test.ipynb

#In[1]
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LatitudeFormatter,LongitudeFormatter
import cartopy.feature as cfeature
import numpy as np
import matplotlib.ticker as mticker

#geoaxesを生成
fig=plt.figure(figsize=(10,5),facecolor='w')
ax=fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=180))

ax.set_global() #全球を表示する
ax.coastlines() #海岸線を描く

#グリッドと軸目盛を描く緯度経度を設定するための配列
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 = mticker.FixedLocator(xticks)    
gl.ylocator = mticker.FixedLocator(yticks)

draw_labelsはFalseとします。Trueにすると一応は軸目盛を描いてくれるのですが、標準では度数表記になってない、度数表記にしたとき経度の中心を180°にしたときEWが重なってしまうなど、調整がややこしいのであきらめました。

軸目盛りを度数表記にする

上述したように、ax.gridlinesからも軸目盛りが書けるのですが、軸目盛りはticksから調整することにします。

#上の続き
#目盛を描く緯度経度の値を設定
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)

見栄えを確認してみる

OUT[1]

軸まわりを調整する処理を関数化する

関数にする際は基本的に引数にax:cartopy.mpl.geoaxesをとるようにすればよいだけです。ただ外部からある程度見た目を調整できるようにしておきたいので、ラインの色や太さなどを調整するための引数も設定します。
またグリッドと軸ラベル両方に対応している図法はメルカトル図法と正距円筒図法の2つだけで、それ以外ではグリッドしか対応していません。汎用性を高めるために渡されたaxの投影法を判定する処理を加えます。

carto_kit.py

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

def set_geogrid(ax,dlon=60,dlat=30
               ,linewidth=0.5,labelsize=15,linestyle='-'
               ,color='grey' ,alpha=0.8 ):
    """
    parameter
    -------------
    ax        :cartopy.mpl.geoaxes
    dlon      :float  grid interval of longitude
    dlat      :float  grid interval of latitude
    linewidth,fontsize,labelsize,alpha :float
    color     :string

    return 
    -------------
    ax
    """
    ax.coastlines()
    gl = ax.gridlines(crs=ccrs.PlateCarree(),color=color
                      , draw_labels=False,linestyle=linestyle
                      ,linewidth=linewidth, alpha=alpha)
    xticks=np.arange(0,360.1,dlon)
#    xtciks=np.arange(-180,180.1,dlon) #中心を経度0°にするとき
    yticks=np.arange(-90,90.1,dlat)
    gl.xlocator = mticker.FixedLocator(xticks)    
    gl.ylocator = mticker.FixedLocator(yticks)
    if (type(ax.projection)==type(ccrs.PlateCarree())):
        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=labelsize)

    return ax

メルカトル図法に対してはyticksが90°をとれないので対応できていません。

ついでに

陸・海・川・湖を塗り分けて描くadd_featureもまとめて関数にしてしまいます。(ほぼほぼギャラリーのコピペですが…)

carto_kit.py

#上のファイルに追記
def set_feature(ax):
    '''
    parameter
    -----------
    ax    :cartopy.mpl.geoaxes
    return
    ----------
    ax    :as above
    '''
    ax.add_feature(cfeature.LAND)
    ax.add_feature(cfeature.OCEAN)
    ax.add_feature(cfeature.COASTLINE)
#    ax.add_feature(cfeature.BORDERS, linestyle=':')
    ax.add_feature(cfeature.LAKES, alpha=0.5)
    ax.add_feature(cfeature.RIVERS)
    return ax

関数を使って地図をプロットしてみる

私は自作の関数を使用するときは、その都度ファイルが保存されているディレクトリを一時的にモジュールサーチパスに追加してからimportするようにしています。(特にWindows環境で環境変数をいじる勇気がないので...)
ここではcarto_kit.pyを~/Pylib下に保存したことにします。

carto_kit.py(再掲)

~/Pylib/carto_kit.py

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

def set_geogrid(ax,dlon=60,dlat=30
               ,linewidth=0.5,labelsize=15,linestyle='-'
               ,color='grey' ,alpha=0.8 ):
    """
    parameter
    -------------
    ax        :cartopy.mpl.geoaxes
    dlon      :float  grid interval of longitude
    dlat      :float  grid interval of latitude
    linewidth,fontsize,labelsize,alpha :float
    color     :string

    return 
    -------------
    ax
    """
    ax.coastlines()
    gl = ax.gridlines(crs=ccrs.PlateCarree(),color=color
                      , draw_labels=False, linestyle=linestyle
                      , linewidth=linewidth, alpha=alpha)
    xticks=np.arange(0,360.1,dlon)
    yticks=np.arange(-90,90.1,dlat)
    gl.xlocator = mticker.FixedLocator(xticks)    
    gl.ylocator = mticker.FixedLocator(yticks)
    if (type(ax.projection)==type(ccrs.PlateCarree())):
        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=labelsize)
    return ax

def set_feature(ax):
    '''
    parameter
    -----------
    ax    :cartopy.mpl.geoaxes
    return
    ----------
    ax    :as above
    '''
    ax.add_feature(cfeature.LAND)
    ax.add_feature(cfeature.OCEAN)
    ax.add_feature(cfeature.COASTLINE)
#    ax.add_feature(cfeature.BORDERS, linestyle=':')
    ax.add_feature(cfeature.LAKES, alpha=0.5)
    ax.add_feature(cfeature.RIVERS)
    return ax

plot_map.ipynb

#In[1]
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import sys
sys.path.append("~/Pylib") #~/Pylibをpythonのモジュールサーチパスに追加
import carto_kit as ckit

fig=plt.figure(figsize=(5,2.5),facecolor='w')
ax=fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=180))
ax.set_global()
ckit.set_geogrid(ax,linewidth=1,color='k',linestyle='--',labelsize=10)
ckit.set_feature(ax)

OUT[1]

余談(特定の領域だけを描画したときの挙動)

ax.set_extent([lon1,lon2,lat1,lat2],crs=ccrs.PlateCarree())でlon1~lon2,lat1~lat2の範囲の地図を描画することができます。ところがax.set_extentのあとにticksの指定をしてしまうと領域の定義が更新されてしまい、グローバルに表示されてしまうようです。

map_plot.ipynb

#In[2]
fig=plt.figure(figsize=(10,5),facecolor='w')
ax=fig.add_subplot(1,2,1,projection=ccrs.PlateCarree(central_longitude=180))
ax.set_extent([0,180,0,90],crs=ccrs.PlateCarre())
ckit.set_geogrid(ax)
ckit.set_feature(ax)
ax.set_title('ax',fontsize=18)

ax2=fig.add_subplot(1,2,2,projection=ccrs.PlateCarree(central_longitude=180))
ckit.set_geogrid(ax2)
ckit.set_feature(ax2)
ax2.set_extent([0,180,0,90],crs=ccrs.PlateCarree()) #軸目盛の調整を終えてから描画する領域を定義
ax2.set_title('ax2',fontsize=18)

OUT[2]

またax.set_extent([lon1,lon2,lat1,lat2],crs=ccrs.PlateCarree())のlon1とlon2は異なる値でないと地図が棒になります。特にlon1=0,lon2=360とした場合、地図上では同じ値とみなされるようです。lon2=360.1などとすることにより棒は回避できます。

追記：マルチプロットに対応する

地図をマルチプロットする際に、緯度経度のラベルを書かないようにするオプションを加えます。
以下set_geogridの引数にleftとbottomをbool値として加え、関数内部でmajor_formatterにplt.NullFormatterを渡すかどうかの判定をすることによって、軸ラベルの非表示を実現させています。

carto_kit.py(改訂)

# coding: utf-8

import cartopy.crs as ccrs
from cartopy.mpl.ticker import LatitudeFormatter,LongitudeFormatter
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.feature as cfeature
import numpy as np
import matplotlib.ticker as mticker
import matplotlib.path as mpath
import matplotlib.pyplot as plt

def set_geogrid(ax,resolution='110m'
               ,dlon=60,dlat=30
               ,manual_ticks=False,xticks=None,yticks=None
               ,bottom=True,left=True,right=False,top=False
               ,linewidth=0.5,fontsize=15,labelsize=15
               ,color='grey' ,alpha=0.8,linestyle='-' ):
    """
    parameter
    -------------
    ax        :cartopy.mpl.geoaxes
    dlon      :float  grid interval of longitude
    dlat      :float  grid interval of latitude
    linewidth,fontsize,labelsize,alpha :float
    color     :string
    resolution :string '10m','50m','110m'
    bottom    :boolean draw  xaxis ticklabel
    lfet      :boolean draw  yaxis ticklabel

    return 
    -------------
    ax
    """
#    以下4行で軸ラベルの位置を変更できると思うのですが正常に動作しないのでコメントアウトしています。  
#    plt.rcParams['ytick.left']=plt.rcParams['ytick.labelleft']=left
#    plt.rcParams['ytick.right']=plt.rcParams['ytick.labelright']=right
#    plt.rcParams['xtick.top']=plt.rcParams['xtick.labeltop']=top
#    plt.rcParams['xtick.bottom']=plt.rcParams['xtick.labelbottom']=bottom

    ax.coastlines(resolution=resolution)
    gl = ax.gridlines(crs=ccrs.PlateCarree()
                      , draw_labels=False,
                      linewidth=linewidth, alpha=alpha
                      , color=color,linestyle=linestyle)
    if manual_ticks == False: 
        xticks=np.arange(0,360.1,dlon)
        yticks=np.arange(-90,90.1,dlat)
    gl.xlocator = mticker.FixedLocator(xticks)    
    gl.ylocator = mticker.FixedLocator(yticks)

    if (type(ax.projection)==type(ccrs.PlateCarree())): 
        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)
        if (bottom==False):
            ax.xaxis.set_major_formatter(plt.NullFormatter())
        if (left==False):
            ax.yaxis.set_major_formatter(plt.NullFormatter())
        ax.axes.tick_params(labelsize=labelsize)
    return ax

#set_featureについては変更がないので省略

　プロット例

multi_plot.ipynb

#In[1]

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import sys 
sys.path.append('~/Pylib/')
import carto_kit as ckit

row,col=2,3

fig=plt.figure(figsize=(18,6),facecolor='w')
ax=[0]*(row*col)
ax_num=0
for i in range(row):
    if(i == 1):
        bottom=True
    else:
        bottom=False
    for j in range(col):
        if(j == 0):
            left=True
        else:
            left=False
        ax[ax_num]=fig.add_subplot(row,col,ax_num+1,projection=ccrs.PlateCarree(central_longitude=180))
        ax[ax_num]=ckit.set_geogrid(ax[ax_num],bottom=bottom,left=left,labelsize=12)
        ax[ax_num].set_title('({},{},{})'.format(row,col,ax_num+1))
     #   ax[ax_num].xaxis.set_major_locator(plt.NullLocator()) #目盛り自体を消したいときはNullLocatorが使える
        ckit.set_feature(ax[ax_num])
        ax_num += 1

OUT[1]

軸ラベルの位置は左と下に固定されており、それをon/offするという動作になっています。
本当は軸ラベルの位置もコントロールできるといいのですが、rcParamsを書き換えるやり方では思った通りに動いてくれません。

まとめ

cartopyはプロットの際に座標を緯度、経度で渡せる点でBasemapよりも明快になっていると思います。
ただ軸の調整に関してはBasemapの方が分かりやすかったようにも思います。

何かおかしいところがありましたらご指摘いただけると助かります。

参考URL

Cartopy Gallery