cartopyの利用
7.
現在、Basemapの開発は終了しているので、代替となるcartopyの利用を考えても良いでしょう。ここではBasemapの代用として、似たような機能を提供しているcartopyで作図する方法を紹介します。バージョン0.17と0.18で経度線や緯度線等の作図が大きく変更され、互換性のあるコードを作成することが難しくなったため、本稿はバージョン0.18以降を中心に解説しています。バージョン0.17での挙動については7.5節に記載しました。
7.1 cartopyの基本
7.1.1地図を描く
cartopyはBasemapと似たようにmatlotlib上で地図を表示するためのパッケージです。 Basemap同様に様々な図法に対応しており、等高線や陰影、矢印、矢羽等の作図を行うこともできます。まずはcartopy_sample.pyを使い、正距円筒図法で作図を行なってみます（図７−１−１）。モジュールはimport cartopy.crs as ccrsのようにインポートします。 cartopyの場合もBasemap同様、最初にmatplotlibでプロット領域を作成するため、matplotlibもインポートします。サブプロットを生成する所では、 matplotlibのfig.add_subplotを使います。これまで同様、最初の引数が縦に並べる数、2つ目が横に並べる数、3つ目がサブプロットのうちの何番目に当たるかを表します。ここでは、新たにprojectionオプションを使い、 ccrs.PlateCarree()を与えることで、cartopyの正距円筒図法（6.2.1節参照）を適用しています。ccrs.PlateCarree()のデフォルトでは東経０度が中央に来るようになっていますが、ここでは、central_longitude=180を与えることで東経180度が中央に来るように変更しています。
ここまでの設定では図枠が表示されるだけなので、ax.メソッドのようなインスタンスメソッドを使い作図してみます。まずはax.coastlinesを使い海岸線を描きます。figやaxの作成時、海岸線作図時のそれぞれの挙動が気になる場合には、ipython3 --pylabを起動してcartopy_sample.pyの各行を順番に入力し、どのように表示されるか試しても良いでしょう。 ax.coastlines() # 海岸線を描くなお、実行時にShapelyDeprecationWarningが表示されるため、次のようなコードを追加して表示されないようにしています。このWarningは無視しておいて問題ないですが、将来のcartopy更新時に修正されるものと思います。 import warnings warnings.filterwarnings('ignore')
7.1.2経度線と緯度線を描く
Basemapの時と同様、描いた図に経度線と緯度線を加えることも可能です(図７−１−２)。作図にはcartopy_sample2.pyを使いました。まずax.gridlinesで経度線・緯度線を描きます。正距円筒図法にするため、図法を指定するcrsオプションにccrs.PlateCarree()を与えます。ここで、経度線・緯度線の設定を行うためにax.gridlinesの戻り値をglに格納しています。経度線や緯度線を引く値は、gl.xlocatorとgl.ylocatorで与えます。その際にmatplotlibのtickerに含まれるFixedLocatorを利用しており、gl.xlocator=mticker.FixedLocator(羅度線を描く位置)のように用います。こうすることで、指定した任意の位置に経度線を描くことができます。経度線を引く値np.arange(0, 360, 30)は、開始点の0から終了点の360まで30毎の値を並べることを意味しています。緯度線についても同様に行います。 import matplotlib.ticker as mticker gl = ax.gridlines(crs=ccrs.PlateCarree()) gl.xlocator = mticker.FixedLocator(np.arange(0, 360, 30)) # 経度線gl.ylocator = mticker.FixedLocator(np.arange(-90, 90, 30)) # 緯度線図７−１−２を見ると、西半球の経度線が描かれていないのが分かると思います。バージョン0.18では、経度線を引く値をnp.arange(0, 360, 30)とすると、 180より大きい値が認識されないようです。これを回避するために、np.arange(- 180, 180, 30)として−180～180まで30毎の値を並べるようにしたものが、 cartopy_sample3.pyです。図７−１−３のように、西半球でも経度線が描かれるようになりました。 gl = ax.gridlines(crs=ccrs.PlateCarree()) gl.xlocator = mticker.FixedLocator(np.arange(-180, 180, 30))   # 経度線gl.ylocator = mticker.FixedLocator(np.arange(-90, 90, 30))   # 緯度線
7.1.3経度線と緯度線の見た目を調整する
次に、経度線・緯度線の見た目を変える方法を紹介します。図7 – 1 – 4のように経度線・緯度線を細い赤破線で描くプログラムが、cartopy_sample4.pyです。線の幅はlinewidthで設定し（デフォルト値： 1 ）、細線にするため0.5にしました。線の色はcolor (図３−２−２、図３−２−３参照)、線種はlinestyleです（表３−２−１参照）。それぞれ、赤色（'r'）、破線（'--'）に設定します。不透明度alphaは0.8にしています。なおlinewidthをlw、linestyleをls、colorをcとする省略形は、ここでは使うことができません。 gl = ax.gridlines(crs=ccrs.PlateCarree(), ¥ linewidth=0.5, linestyle='--', color='r', alpha=0.8)これらのオプションを使いlinewidth=1、黒色('r')、点線(':')に設定して、幅１の黒点線で経度線・緯度線を描きます(図７−１−５)。cartopy_sample5.pyを使って作図しました。 gl = ax.gridlines(crs=ccrs.PlateCarree(), ¥ linewidth=1, linestyle=':', color='k', alpha=0.8)なおax.coastlinesにも線の色、線の幅、線種、不透明度を指定するオプションがあり、ax.coastlines(color='g')のように指定できます(緑色に変更する場合)。 7.1.4経度線と緯度線にラベルを付けるax.gridlinesには、経度線と緯度線にラベルを付けるdraw_labelsオプションがあります。draw_labels=Trueとすることでラベルが付きます(図７−１− 6 )。作図にはcartopy_sample6.pyを使いました。 gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=1, linestyle=':', color='k', alpha=0.8)同じスクリプトを使った場合でも、バージョン0.18、0.19と0.20・0.21で見た目が変わり、0.19では経度方向のラベルの向きが、0.18では緯度方向のラベルの有無が他と異なっています。なおdraw_labelsオプションの挙動は、cartopyバージョン0.17と0.18の間で大きく変わりました。バージョン0.17では、うまくラベルの付く範囲が東経０度を中心とした場合に限られており、他の範囲では7.1.5節の方法でカスタマイズする必要がありました。図７−１−６では、緯度線ラベルや経度線ラベルの一部が表示されていません。バージョン0.19では一部のラベルが斜めになっています。このように自動調整では不具合が出てしまいます。図７−１−７は、図のサイズを大きくした場合の見た目で、全てのラベルが表示され、それなりの見た目になったように思えます。バージョン0.18と0.19で同じ見た目になります。作図にはcartopy_sample7.pyを使いました。
fig = plt.figure(figsize=(8, 6)) # 図のサイズを大きくする
上と右のラベルが表示されないように、次のような設定を加えました。右のラベルはここでは表示されていないため、上のラベルだけが消えました。 gl.top_labels = False # 上のラベルを消すgl.right_labels = False # 右のラベルを消すなお同様のオプションが、下のラベル（gl.bottom_labels）や左のラベル(gl.left_labels)にも設定可能です。図の中心が東経０度の場合も試してみましょう(cartopy_sample8.py)。図７ – 1 – 8では、西経180度が左端、東経180度が右端にきました。左端の180° Wが180°のように表示され右端のラベルが表示されないため、次のようにラベルを描く範囲を150°Wから150°Eまでに変更しています。 gl.xlocator = mticker.FixedLocator(np.arange(-150, 180, 30))   # 経度線
7.1.5経度線、緯度線ラベルのカスタマイズ
ここまでは、ax.gridlinesのオプションを使った方法を紹介しましたが、カスタマイズが行いやすい方法として、matlotlibでラベルを描くax.set_xticks、 ax.set_yticks を使う方法も紹介しておきます。この方法では、経度線と緯度線はax.gridlinesで描き、ラベルはこれとは別にax.set_xticks、ax.set_yticksで描きます。まずは、経度線と緯度線を描く値とラベルを描く値を同じに設定します。共通した値を使うので、xticks、yticksという配列を作っておき、それぞれの値を格納しておきます。dlon、dlatは、それぞれ経度線、緯度線ラベルを描く間隔を表しています。緯度線ラベルを90°Sから90°Nまで30°毎に付けるため、 np.arange(-90, 90.1, dlat)としました。作図にはcartopy_sample9.pyを使いました。 # 緯度線、経度線とラベルを描く値dlon, dlat = 30, 30 xticks = np.arange(-180, 180, dlon) yticks = np.arange(-90, 90.1, dlat)まずax.gridlinesで、これまで同様に経度線と緯度線を描きます。 gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, linewidth=1, linestyle=':', color='k', alpha=0.8) gl.xlocator = mticker.FixedLocator(xticks) # 経度線を描く値gl.ylocator = mticker.FixedLocator(yticks) # 緯度線を描く値次に、ax.set_xticks、ax.set_yticksを使い、経度線と緯度線ラベルを描く値を指定します。オプションとしてcrsを与えることができ、ccrs.PlateCarree()を指定します。図７−１−８のように、180度を中心として東経側にも西経側にもラベルが描かれるようになりましたが、経度表記に直したくなります。 ax.set_xticks(xticks, crs = ccrs.PlateCarree()) # 経度線ラベルの値ax.set_yticks(yticks, crs = ccrs.PlateCarree()) # 緯度線ラベルの値cartopy_sample10.pyは、経度線・緯度線ラベルを度数表記にするためのプログラムです。経度線ラベルの表記はax.xaxis.set_major_formatterで変更します。その際に、cartopy.mpl.tickerに含まれるLongitudeFormatterを使って書式変更します。デフォルトのままでは東経０度が「０°E」ではなく「０°」と表示されてしまうため、オプションとしてzero_direction_label=Trueを与えました。図7 – 1 –10のように、経度方向のラベルが180度を中心とした数字から東経(E)、西経(W)を区別した度数表記に変わりました。緯度線ラベルについても同様に、ax.yaxis.set_major_formatterの書式を変更します。緯度の場合はLatitudeFormatterですが、赤道は「０°」表記で良いためデフォルトのまま使用しました。 from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter #目盛りの表示形式を度数表記にするax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True)) ax.yaxis.set_major_formatter(LatitudeFormatter()) 7.1.6経度線の間隔と経度線ラベルの間隔を別々にする図7 – 1 –10では、150°Wと120°Wのように一部のラベルが重なっています。ラベルの文字は少し小さいように思えますが、これ以上拡大すると隣のラベルと判別できなくなってしまうので、経度線は30度毎のまま、経度線ラベルだけ60度毎に描いてみます。経度線・緯度線を描く部分は変更ありませんが、経度線・緯度線ラベルを描くために、新たにdlon_lab、dlat_labで経度線、緯度線ラベルの間隔を指定し、xticks_lab、yticks_labという配列を作って、経度線、緯度線ラベルを描く値を格納しておきます。経度線ラベルの間隔のみ60度にしておくことで、経度線の間隔は30度のままラベルの間隔を変更できます。作図にはcartopy_samplel l.pyを用いました。 dlon_lab, dlat_lab = 60, 30 # 経度線・緯度線の間隔xticks_lab = np.arange(-180, 180, dlon_lab) # 経度線ラベルを描く値yticks_lab = np.arange(-90, 90.1, dlat_lab) # 緯度線ラベルを描く値経度線ラベル、緯度線ラベルを描く部分では、これらxticks_lab、yticks_labを使います。 ax.set_xticks(xticks_lab, crs=ccrs.PlateCarree()) # 経度線ラベル|ax.set_yticks(yticks_lab, crs=ccrs.PlateCarree()) # 緯度線ラベル目盛り線ラベルのサイズは、ax.axes.tick_paramsで設定します。オプションとしてlabelsize=12とすれば、文字サイズが12ポイントに設定されます。 |ax.axes.tick_params(labelsize=12) # 目盛り線ラベルのサイズを指定ax.axes.tick_paramsでは、他にも多くの項目を変更可能です(表7 – 1 – 1 、図7 – 1 – 12)。オプションのうち、axisとwhichは変更をx軸とy軸のどちらに適用するや主目盛り、副目盛りのどちらに適用するかを指定するものです。指定がない場合は全てに同じ変更が適用されます。目盛り線の設定は、幅をwidth (デフォルト：１)、長さをlength (デフォルト：3.5)、色をcolor（デフォルト：黒）で変更できます。目盛り線の向きは、デフォルトでは枠線の外側(direction='out')ですが、内側に変更することもでき、その場合にはdirection='in'とします。目盛り線ラベルは、目盛り線とは別に設定するようになっており、色をlabelcolor（デフォルト：黒）、ラベルを回転させる角度をlabelrotation（デフ
ax.axes.tick_params(オプション)：変更したい項目を指定ax.axes.tick_params(オプション, axis='x')：x軸で変更したい項目を指定ax.axes.tick_params(オプション, axis='y')：y軸で変更したい項目を指定ax.axes.tick_params(オプション, axis='x', which='major')：x軸の主目盛りについて、変更したい項目を指定図枠の下側、上側、左側、右側に目盛り線を付けるかどうかを決めるbottom、 top、left、rightオプションがあり、デフォルトでは下側と左側だけに目盛り線が付きます。目盛り線ラベルを描くかどうかは、labelbottom、labeltop、labelleft、 labelrightオプションで設定し、こちらもデフォルトでは下側と左側だけに付くようになっています。これらのオプションを使い、図枠の上側と右側にも目盛り線と目盛り線ラベルを付けてみます(図７−１−１３)。作図には、cartopy_sample12.pyを使いました。上側と右側に目盛り線を付けるため、top=True、right=Trueとしています。また上側と右側の目盛り線ラベルは、labeltop=True、labelright=Trueで付けます。 ax.axes.tick_params(labelsize=12, top=True, right=True, labeltop=True, labelright=True)
7.1.7陸と海を塗り分けて描く
cartopyでもBasemap同様に陸と海を塗り分けるツールが提供されています。ここでは、cartopy.featureを使って塗り分けてみます（図７−１−14)。作図にはcartopy_sample13.pyを使います。ここでは、変更される部分が分かりやすくなるように、プログラム中のax.coastlines()を無効にしているため、海岸線は描かれていません。cartopy.featureを使って作図する際には、初回のみデータのダウンロードが行われるため、ネットワークに接続されていない環境ではcartopy_sample13.pyでエラーが発生します。最初に、作図に必要なcartopy.featureをimportしてcfeatureとして利用します。ax.add_feature(cfeature.塗り潰す対象)で色を付けることができ、陸に色を付ける場合はcfeature.LAND、海に色を付ける場合がcfeature.OCEANです。なおcfeature.OCEANでは、海洋に加えて湖のうちカスピ海だけは塗り潰されています。 import cartopy.feature as cfeature ax.add_feature(cfeature.LAND) # 陸に色を付けるax.add_feature(cfeature.OCEAN) # 海に色を付ける
7.1.8海岸線と湖、川を描く
さらに、海岸線と湖、川を追加していきます（図７−１−15)。海岸線を追加するにはcfeature.COASTLINE、湖を描くにはcfeature.LAKES、川を描くにはcfeature.RIVERSを使います。黒実線の海岸線に加え、五大湖やバイカル湖、アマゾン川、ナイル川など主要な湖や河川が水色で表示されるようになりました。作図には、cartopy_sample14.pyを使いました。なお海岸線を描く場合、デフォルトではax.coastlinesよりは太い線になります。線の太さはlinewidthオプションを使うことにより変更可能で、少し細い線にするためlinewidth=0.8としています。他には不透明度を表すオプションalpha（0～１の範囲、１で不透明、０で透明）も利用可能です。 ax.add_feature(cfeature.COASTLINE, linewidth=0.8) # 海岸線を描くax.add_feature(cfeature.LAKES) # 湖を描くax.add_feature(cfeature.RIVERS) # 川を描く
7.1.9陸上に国境線を描く
陸上に国境線を描くこともでき、cfeature.BORDERSを使います。点線で国境線が描かれるようになりました(図７−１−１６)。破線にするためのオプションとしてlinestyle='--'、線の幅を細くするためlinewidth=0.8を加えています。作図には、cartopy_sample15.pyを使いました。 ax.add_feature(cfeature.BORDERS, linestyle='--', linewidth=0.8) # 国境線
7.1.10陸や海、湖を指定した色で塗り潰す
cartopy.featureはデフォルトでも綺麗に陸と海が塗り潰されているので、変更する必要はないかもしれませんが、色を指定して塗り潰すことも可能になっています。その際には、オプションとしてcolorを使います。ここでは、Basemapの時に陸と海の塗り分けを行った図６−１−９と同じ色を付けてみます（図７− l – 17)。陸を緑色にするためにcolor='g'、海と湖を水色にするためにcolor='aqua'としています。湖の周りが黒線にならない違いはありますが、図６ – 1 – 9と似た図になりました。作図には、cartopy_sample15.pyを使いました。 ax.add_feature(cfeature.LAND, color='g') # 陸を緑色で塗り潰すax.add_feature(cfeature.OCEAN, color='aqua') # 海を水色で塗り潰すax.add_feature(cfeature.LAKES, color='aqua') # 湖を水色で塗り潰すなお、川と海岸線、国境線を描く所でもcolorオプションは有効ですが、意図していない部分まで塗り潰されてしまうので、使わない方が良いでしょう。
7.2 cartopyを使った等高線、陰影の作図
ここでは、Basemapの所で行ったように、気象データの形式の１つであるNetCDFデータを読み込み、作図を行います。気象データの形式や読み込みについては、6.4節を参照してください。6.5.3節で行った時と同様に、７月の平均値を等高線、2018年７月の海面更正気圧（SLP）偏差を陰影で作図できるようにしていきます。作図には、月平均したSLPデータのNetCDFファイル(slp.mon.mean.nc)が必要なので、まだダウンロードしていない場合には、 basemap_nc2bin.pyでダウンロードしておきます。
7.2.1等高線を描く
まずはSLPデータを利用して等高線を作成します（図７−２−１）。作図にはcartopy_contour.pyを使いました。basemap_contour6.pyで図６−５−６を作図した時と同様に、NCEP再解析データのNetCDFファイルを読み込み、取り出した1次元の経度データをlon、 1次元の緯度データをlat、空間2次元、時間1次元のSLPデータをslpという変数に格納しています。次も同様に、７月のSLP平均値を計算してslp_julに格納します。時間平均を行なったので、slp_julは空間２次元のデータです。プログラム中では、作成したデータのサイズをnumpy配列名.shapeを使って書き出しており、lon: (144,)、lat: (73,)、slp: (860, 73, 144)、slp_jul: (73, 144)のように、lonとlatは1次元、slpは3次元、slp_julは2次元のデータとなっていることが分かります。このメソッドでは、配列の形状をタプルで書き出すので、1次元の場合でもカンマが付いています。 tstr = (1981 - 1948) * 12 + 6 # 開始点：1981年7月tend = (2010 - 1948 + 1) * 12 # 終了点：2010年の12月slp_jul = slp[tstr:tend:12, :, :].mean(axis=0) # 7月のSLP平年値cartopyで等高線を描くコマンドは、ax.contour(引数1,引数2,引数3)です。5.2節で紹介したmatplotlibのplt.contourと同じですが、axを定義した際にcartopyの正距円筒図法を適用したので、正距円筒図法上で等高線を描くようになっています。5.2節の場合と同様、ax.contourのl番目の引数はx軸上の位置、２番目の引数はy軸上の位置で、３番目の引数は等高線を描くためのz軸データに当たります。z軸データとしては、2次元のslp_julを用います。 x, y軸データとしては、 l次元の経度・緯度データlon、latからnp.meshgridを使って作成した２次元データx、yを使います。x、yも形状を書き出しており、x: (73, 144)、y: (73, 144)のように２次元であることが確認できます。４番目の引数として、等高線を描く値を与えることができ、そのclevsは、5.2節のときと同様に４ hPa毎に設定しています。４番目の引数を与える代わりにlevels=clevsとしても同じです。 clevs = np.arange(np.floor(slp_jul.min() - np.fmod(slp_jul.min(), 20)), ¥ np.ceil(slp_jul.max()) + 1, 4) ax.contour(x, y, slp_jul, clevs, linewidths=0.8, colors='k')等高線を描く際には、5.2節で紹介したplt.contourのオプションが利用できます（表５−２−１）。ここでは、線の幅をlinewidths=0.8（デフォルト１）、線の色をcolors='k'（黒色）に設定しました。
7.2.2等高線ラベルを描く
図６−５−６のような等高線になりましたが、等高線に気圧を示すラベルが付いていません。ラベルを付けるには、Basemapで等高線ラベルを付けた時と同様に、等高線を描いた値cs.levelsをclevelsに格納し、ax.contourの戻り値csにある、cs.clabelにclevelsを渡す方法を使います（cartopy_contour2.py）。ここでは、clevels[::5]を与えて５飛ばし(4 hPa毎に等高線を描いたので20 hPa毎)に等高線ラベルを描いています(図７−２−２)。ラベルの付く場所は、cartopyのバージョンによって変わります。ここではバージョン0.20のものを用いています、 cs = ax.contour(x, y, slp_jul, clevs, linewidths=0.8, colors='k') # 等高線clevels=cs.levels # ラベルを付ける値cs.clabel(clevels[::5], fontsize=12, fmt="%d") # 等高線ラベル
7.2.3陰影を描く
それでは、2018年７月のSLP偏差を陰影で描きます（図７−２−３）。図６− 5 – 7とほとんど同じ図になりました。作図にはcartopy_contourf.pyを用いました。6.5.3節同様に2018年7月の偏差を計算し、slp_anomに格納しておきます。陰影の範囲の決め方も同様に、正が赤、負が青で領域の範囲の中央が０になるようにしており、clevsの範囲の両端は、偏差の絶対値の最大値よりも1大きい値にしています。陰影を描くコマンドは、ax.contourf(凱数1,凱数2,凱数3)です。 matplotlibのplt.contourf と同じなので、5.3節で紹介したplt.contourfのオプションが利用できます(表５−３−２)。1番目から３番目までの引数は等高線を描くax.contourと同じで、x軸上の位置、y軸上の位置、z軸データを与えます。4番目の引数には陰影を描く値clevsを与えていますが、等高線の場合と同様、引数を与える代わりにlevels=clevsとも記述できます。陰影を描く色はcmapオプションで色テーブルの名前から設定することができ、cmap='bwr' (青～白～赤と変化）としました。色テーブルの名前は、図４−５−７にあるものが利用できます。 n = (2018 - 1948) * 12 + 6 slp_anom = slp[n, :, :] - slp_jul # 2018年７月の偏差work = max(np.abs(np.floor(slp_anom.min())), np.abs(np.ceil(slp_anom.max())) clevs = np.arange(-(work+1), work + 1) # 陰影を描く値ax.contourf(x, y, slp_anom, clevs, cmap="bwr") # 陰影を描く
7.2.4カラーバーを付ける
作成した図にカラーバーを付け、陰影の値が分かるようにしていきます。作図にはcartopy_contourf2.pyを用いました。カラーバーを付けるには、5.3.3節、6.5.3節で行なったようにplt.colorbarを使います。cartopyの場合には、これまでのようにplt.contourf、m.contourfの直後にplt.colorbarを行っても、カラーバーの作成に必要な情報がないためエラーになります。代わりにax.contourfの戻り値をimに格納しておき、plt.colorbarのl番目の引数として渡すことでカラーバーを作成します。カラーバーはデフォルトで縦方向ですが、そのままでは図枠の上下に大きくはみ出してしまうので、shrink=0.55で大きさの調整をしています。カラーバーにラベルを付ける方法は、これまでと同じでplt.colorbar の戻り値cbarに対してcbar.set_labelを使います。 im = ax.contourf(x, y, slp_anom, clevs, cmap="bwr") cbar=plt.colorbar(im, shrink=0.55) cbar.set_label('SLP anom. in Jul. 2018', fontsize=14) plt.colorbarでは大きさの調整の他にも、表５−３−１に載せた様々なオプションを使用することができます。これらのオプションを使いカラーバーを水平にして大きさやカラーバーとグラフの間の間隔を調整したものが図７−２−５です。作図にはcartopy_contourf3.pyを用いました。カラーバーは、shrink=0.8で大きさを調整し、水平方向にするために orientation='horizontal'、カラーバーとグラフの間の間隔を調整するためにpad=0.08を使いました。 cbar = plt.colorbar(im, shrink=0.8, orientation='horizontal', pad=0.08) cartopyでもmatplotlib同様に、ax.set_titleを使い図のタイトルを付けることが可能になっています（図７−２−６）。ここではfontsizeオプションを使い文字を大きくしています。作図には、cartopy_contourf4.pyを使いました。 ax.set_title("Jul. SLP & anom. in Jul. 2018", fontsize=20)
7.2.4等高線の体裁の調整
図7 − 2 − 6では、1000 hPaや1020 hPaの等高線にラベルが付いていますが、等高線の幅は全て同じになっているため、ラベルのないところでは等高線の値が読み取りにくくなっています。そこで、20 hPa毎に太線にして等高線の値がわかるように変えてみたいと思います（cartopy_contourf5.py）。等高線の間隔は4 hPaとなっており、5飛ばしにラベルを付けているため、等高線の幅も5飛ばしに太くします。等高線の幅を指定するlinewidthsオプションでは、 linewidths=0.8のような数値での指定の他に、等高線の幅をリストで与えることも可能です。そこで、linewidths=[1.2, 0.8, 0.8, 0.8, 0.8]のように太線、細線、細線、細線、細線、のリストを与えると、５飛ばしに太線が描かれるようになります(図７−２−７)。 cs = ax.contour(x, y, slp_jul, clevs, linewidths=[1.2, 0.8, 0.8, 0.8, 0.8], colors='k')等高線の幅の他に、色（colors）や線種（linestyles）もリストで与えることができ、これらの組み合わせで多様な等高線を描くことが可能です。例えば、 linewidths=[1.2, 0.8, 1.2, 0.8]、linestyles=["-", "-", "--", "--"]では、太実線、細実線、大破線、細破線の順で描かれます。 7.3 cartopyを使った矢羽、矢印の作図cartopyではBasemap同様に矢羽や矢印を描くことができます。ここでは、 7.2節で作成した７月のSLP平均値の等高線に矢羽や矢印を描いていきます。作図にはSLPデータのNetCDFファイルに加え、東西風、南北風データのNetCDFファイル(uwnd.mon.mean.nc、vwnd.mon.mean.nc)も必要なので、ダウンロードしていない場合にはbasemap_nc2bin.pyでダウンロードしておおます。
お
7.3.1矢羽を描く
矢羽を描くには、ax.barbs(引数1,引数2,引数3,引数4)を使います。 1番目の引数がx軸上の位置、 2番目の引数がy軸上の位置で、 3番目の引数が東西風データ、4番目の引数が南北風データに当たります。5.1.1節のplt.barbsと同じなので、表５−１−１のオプションが利用可能です。矢羽を描くためのサンプルプログラムは、cartopy_barbs.pyです。SLP、東西風、南北風データの３ファイルを読み込むために、file_namesとvar_namesにそれぞれのファイル名とNetCDFファイルに格納されたデータ名を識別するための変数名を定義しておきます。これらをzipに渡してループを回し、flle_name、var_nameに同じ番号の要素がセットで入るようにしています。nc.variablesで読み込んだデータは、Numpyのndarrayになっており、それらをループが回る毎にリストdに追加していきます。そのため、リストの要素d[0]がSLP、d[1]が東西風、d[2]が南北風に対応します。次に SLP、東西風、南北風それぞれについて、1981～2010年の７月平均値を計算します。例えばSLPの場合、d[0][tstr:tend:12, :, :]のように1981～2010年の７月のデータを取り出します。最初のd[0]は、リストの０番目の要素を取り出すことを意味します。取り出されたd[0]はndarrayになっているので、時間軸(axis=0)にスライス記法[tstr:tend:12, :, :]を適用することで、７月のデータを取り出すことができます。最後の.mean(axis=0)はこれまで同様、時間軸について算術平均を行うことを意味しています。 7.2.1節、7.2.2節と同様に等高線を描いてラベルを付けた後、ax.barbsにcolor='r'オプションを与えることで赤色の矢羽を描いています（図７−３−１）。また、6.5.5節のBasemapの所で解説した矢羽を間引く方法（３グリッド毎に描くために[::3, ::3]の記法を使う)、風速の小さい場所に丸印が描かれないようにする方法（sizes=dict(emptybarb=0.0)オプション）も使いました。なお、矢羽を描くために用いる東西風、南北風データの単位がm/sなので、短矢羽が5 m/s、長矢羽が10 m/s、旗矢羽が50 m/sで描かれます。 ax.barbs(x[::3, ::3], y[::3, ::3], u[::3, ::3], v[::3, ::3], sizes=dict(emptybarb=0.0), length=4, color='r')
7.3.2矢印を描く
矢印を描くには、ax.quiver( l数1,  l数2,  l数3,  l数4 )を使います。引数は矢羽の場合と同じで、x軸上の位置、y軸上の位置、東西風データ、南北風データになります。矢印を描くプログラムはcartopy_quiver.pyです。ax.quiverは5.1.4節のplt.quiverと同じなので、表５−１−２のオプションが利用可能です。ここでは、色を指定するcolorオプションを付けて赤矢印を描きました（図7 − 3 − 2 )。 ax.quiver(x[::3, ::3], y[::3, ::3], u[::3, ::3], v[::3, ::3], color='r')これだけでは矢印の大きさが分からないので、凡例を付けてみます。まずは、矢印の幅と長さを固定するため、width=.003、scale=300.を指定しました。作図にはcartopy_quive2.pyを使いました。 Q = ax.quiver(x[::3, ::3], y[::3, ::3], u[::3, ::3], v[::3, ::3], ¥ color='r', width=.003, scale=300.)矢印の凡例を描くには、ax.quiverkeyを用います。 l番目の引数のＱはax.quiverの戻り値です。２番目と３番目の引数は、図枠の経度・緯度範囲のどの場所に凡例を描くのかを表します。図枠の範囲が０～１に対応し、左下が(0.0, 0.0)、右上が(1.0, 1.0)です。4番目の引数の整数10は、10m/sの矢印を描くことを表し、５番目の引数の文字列'10'は、矢印に表示する文字列を表します。 labelposは、矢印に対してどの位置に文字列を描くのかを意味し、labelpos='N'で矢印の上側、labelpos='S'で下側に描かれます(デフォルトではlabelpos='N')。矢印と文字列との間隔はlabelsepで指定します(デフォルトはlabelsep =0.1)。デフォルトでは矢印と文字列が離れすぎているので、0.01と小さな値にしました。図７−３−３では、図枠の右上に横向きの矢印とその上に１０の文字が描かれています。 ax.quiverkey(Q, 0.97, 0.93, 10, '10', labelpos='N', color='r', labelsep=0.01)この状態では凡例の判別が難しいので、凡例の近くを塗り潰してみます（図7 – 3 – 4 )。凡例の周りが灰色の背景になったのが分かるでしょうか。ここまでは矢印を赤で描いていましたが、黒に変えました。また、矢印の頭の部分の幅を大きくして方向が判別しやすいように変えました（headwidth=3.75）。図の上に来るようにするためzorder=4としています。作図にはcartopy_quive3.pyを使いました。 Q = ax.quiver(x[::3, ::3], y[::3, ::3], u[::3, ::3], v[::3, ::3], ¥ color='k', width=.003, scale=300., headwidth=3.75, zorder=4)凡例の周りの背景を灰色に塗り潰すため、ax.add_patchで矩形領域を追加しました。矩形領域の設定は、plt.Rectangleで行いました。最初の(155, 70)は東経155度、緯度70度を左下の始点とすることを意味し、次の25と20は、経度方向に25度、緯度方向に20度の幅を持つことを意味します。そのような設定なので、東経155～180度、緯度70～90度が塗り潰されることになります。塗り潰す色はfacecolorで指定し、ここでは灰色（facecolor='gray'）としました。塗り潰した領域の端には枠線を付けないようにするため、edgecolor=Noneとしています。 rect = plt.Rectangle((155, 70), 25, 20, facecolor='gray', ¥ edgecolor=None, zorder=4) ax.add_patch(rect)図枠の内側が塗り潰されると困る場合もあると思います。そのような場合には、矢印の凡例を図枠の外に出すことも可能です（図７−３−５）。作図にはcartopy_quive4.pyを使いました。凡例の緯度方向の位置を表す３番目の引数にマイナスの値を与えると、図枠の下側に表示されます。矢印のラベルを矢印の下側に表示するため、labelpos='S'としました。 ax.quiverkey(Q, 0.97, -0.04, 10, '10', labelpos='S', color='k', labelsep=0.03)ラベルの文字の大きさを変えることも可能です（図７−３−６）。ここでは判例に用いる風速を15m/sに変え、文字の大きさを12にしました。設定は、 fontproperties={'size': 12}です。作図はcartopy_quive5.pyで行いました。なお、coordinates='axes'はデフォルト値で、省略しても同じですが、x軸上の位置、y軸上の位置を図枠の範囲で0 ～ 1として指定することを意味します。 ax.quiverkey(Q, 0.97, -0.04, 15, '15', labelpos='S', color='k', labelsep=0.03, coordinates='axes', fontproperties={'size': 12}) 7.4 cartopyで利用できる図法cartopyでもBasemapのように、様々な図法を利用して作図することができるようになっています。ここではcartopyで作図する方法についてのみ解説するので、図法については6.2節の解説を参照してください。
7.4.1正射投影図法
6.2.2節で紹介した正射投影図法（あるいは平射図法）で作図する方法です。作図を行うプログラムは、cartopy_proj.pyです。正射投影図法で作図するには、 projectionオプションにccrs.Orthographic()を与えます。図７−４−１のように地球儀風の図ができました。Basemapの所で行った図６−２−１の右上の図の場合と同様に、東経180度、北緯45度を中心に描いたので、北太平洋が手前側に来ています。正射投影図法の場合にも経度線・緯度線を描く間隔を指定でき、次のように経度線・緯度線を10度毎に付けるようにしました。ラベルを付けることはできないため、ラベルの設定はありません。さらに、ax.set_titleで図にタイトルを付けました。 ax.set_title("projection='ccrs.Orthographic'") # タイトルを付ける
7.4.2メルカトル図法
6.2.3節で紹介したメルカトル図法でも作図することができます(図７−４− 2）。作図を行うプログラムは、cartopy_proj2.pyです。メルカトル図法で作図するには、projectionオプションにccrs.Mercator()を与えます。ここでは、太平洋側を中心にするため中心経度central_longitude=180.0とし、南緯60度から北緯60度の範囲で描くため、min_latitude=-60.0、max_latitude=60.0としました。経度線は30度毎、経度線は15度毎に描きました。 ax = fig.add_subplot(1, 1, 1, ¥ projection=ccrs.Mercator(central_longitude=180.0, ¥ min_latitude=-60.0, ¥ max_latitude=60.0)) ax.set_title("projection='ccrs.Mercator'")メルカトル図法の作図では、東経０度が中心の図では、ax.gridlinesでdraw_labels=Trueとすれば、目盛り線ラベルを付けることができます。中央を東経180度とした場合には、経度帯の半分しか目盛り線ラベルがつかないため、ここでは、ラベルを付けていません。
7.4.3ランベルト正角円錐図法
6.2.4節で紹介したランベルト正角円錐図法も試してみます(図７−４−３)。作図を行うプログラムは、cartopy_proj3.pyです。ランベルト正角円錐図法で作図するには、projectionオプションにccrs.LambertConformal()を与えます。ランベルト正角円錐図法では、図の中心経度と中心緯度の設定によって図の歪みなどが変わってくるため、対象となる領域に近い点を指定した方が良いでしょう。ここでは、6.2.4節の場合と同じ、中心経度central_longitude=180.0、中心緯度central_latitude=60.0を使いました。 ax = fig.add_subplot(1, 1, 1, ¥ projection= ccrs.LambertConformal(central_longitude=180.0, ¥ central_latitude=50.0)) ax.set_title("projection='ccrs.Mercator'")ランベルト正角円錐図法の場合には、ax.gridlinesでdraw_labels=Trueとした場合、バージョン0.19では緯度・経度のラベルが表示されますが、古いバージョンではエラーとなります。ランベルト正角円錐図法では両極を同時に描けないので、地球全体を把握するための図よりは領域を限定した図に適しています。ax.set_extentを使うことで日本周辺に限定することもできます(図７−４−４)。作図はcartopy_proj4.pyで行いました。中心経度は東経135度、中心緯度は北緯35度に設定しました。 dlon、dlatを変更し、経度線、経度線は10度毎に描きました。 ax = fig.add_subplot(1, 1, 1, ¥ projection=ccrs.LambertConformal(central_longitude=135.0,¥ central_latitude=35.0)) ax.set_extent([100, 170, 10, 70]) # 領域の限定
7.4.4極投影図法
6.2.5節で紹介した極投影図法（ポーラーステレオ図法）で描くことも可能です（図７−４−５）。作図はcartopy_proj5.pyで行いました。南極を中心とする極投影図法で作図する場合、projectionオプションにccrs.SouthPolarStereo()を与えます。北極を中心とする場合には、ccrs.NorthPolarStereo()です。極投影図法では作図する領域を指定する必要があり、ここでは、南極から南緯20度までの全経度帯をax.set_extent([-180, 180, -90, -20], ccrs.PlateCarree())で指定しています。ここでccrs.PlateCarree()を省略するとうまく作図できません。 dlon, dlat = 30, 10に設定して経度線を30度毎、緯度線を10度毎に描きました。指定したのは南緯20度までですが、図7 − 4 − 5のように、南緯20度の緯度円に外接するような四角形の領域まで描かれます。 ax = fig.add_subplot(1, 1, 1, projection=ccrs.SouthPolarStereo()) ax.set_extent([-180, 180, -90, -20], ccrs.PlateCarree()) ax.set_title("projection='ccrs.SouthPolarStereo'") cartopyでは、作図範囲を南緯20度の緯度円までに限定することができ、 cartopy_proj6.pyで行っています(図７−４−６)。cartopyを呼び出した際のサブプロット(軸：axes)は、左下の座標が(0.0, 0.0)、右上の座標が(1.0, 1.0)です（図７−４−７)。サブプロットの端が図７−４−５の図枠に対応しており、この領域はax.transAxesで参照することができます。そのため、中心座標を(0.5, 0.5）とした半径0.5の円を作成すれば、この円は図枠の四角形と接することが分かります。円形の領域内だけ表示させることができれば、南緯20度の緯度円に限定できるようになります。ここでは、matplotlibに含まれているmatplotlib.pathを使い、このモジュールをmpathで参照しています。 mpath.Pathに円形の領域を渡して領域内だけ表示させることで、南緯20度の緯度円に限定しています。まずcenterに中心（0.5, 0.5）の座標を、radiusに半径0.5を入れています。次に0 ～2πの範囲で位相thetaを作成し、中心(0.0, 0.0)、半径１の円を構成する座標を計算して、vertsに格納しています。最後に、 verts * radius + centerの計算を行った結果をmpath.Pathに渡し、円形の領域を作成します(circle)。ax.set_boundaryは、図枠の領域を限定するために使っており、範囲としてサブプロットの四角形の領域に内接する円形の領域 circleを渡すことで、南緯20度の緯度円になります。
7.4.5モルワイデ図法
6.2.6節で紹介したモルワイデ図法も試すことができます(図７−４−８)。作図には、cartopy_proj7.pyを使いました。モルワイデ図法の場合、projectionオプションにccrs.Mollweide()を与えます。太平洋側を中心にするために、中心経度central_longitude=180.0としました。dlon, dlat = 30, 30と設定し、経度線、経度線を30度毎に描きました。 ax = fig.add_subplot(1, 1, 1, ¥ projection=ccrs.Mollweide(central_longitude=180.0)) ax.set_title("projection='ccrs.Mollweide'")
7.4.6 ロビンソン図法
6.2.7節で紹介したロビンソン図法も試すことができます(図７−４−９)。ロビンソン図法の場合には、projectionオプションにccrs.Robinson()を与えます。 cartopy_proj8.pyで作図しました。 ax = fig.add_subplot(1, 1, 1, ¥ projection=ccrs.Robinson(central_longitude=180.0)) ax.set_title("projection='ccrs.Robinson'")
7.4.7ランベルト正積円筒図法
6.2.8節で紹介したランベルト正積円筒図法も利用できます（図７−４−１０)。作図には、cartopy_proj9.pyを用いました。ランベルト正積円筒図法の場合には、projectionオプションにccrs.LambertCylindrical()を与えます。dlon, dlat = 30, 15と設定し、経度線を30度毎に、経度線は極域ほど幅が狭くなっていることが分かるように15度毎に描きました。 ax = fig.add_subplot(1, 1, 1, ¥ projection=ccrs.LambertCylindrical(central_longitude=180.0)) ax.set_title("projection='ccrs.LambertCylindrical'")
7.4.8ミラー図法
6.2.9節で紹介したミラー図法で描く機能もあります(図７−４−１１)。ミラー図法で描くには、projectionオプションにccrs.Miller()を与えます。赤道付近は図７−４−２と似ていますが、両極まで描くことが可能になります。作図には、 cartopy_proj10.pyを用いました。 ax = fig.add_subplot(1, 1, 1, projection=ccrs.Miller(central_longitude=180.0)) ax.set_title("projection='ccrs.Miller'")
7.5 cartopyバージョン0.17
バージョン0.17から0.18への更新で経度線や緯度線等の作図が変更され、互換性のあるコードを作成することが難しくなりました。バージョン0.18以降の便利な機能を利用するのが望ましいですが、以前のバージョンしか利用できない環境のことも考え、本節ではバージョン0.17で作図を行う場合のポイントについてまとめています。
7.5.1経度線や緯度線作図時の挙動
7.1.2節で行った経度線と緯度線の作図は、バージョン0.17でもax.gridlinesで行い、経度線や緯度線を引く値は、gl.xlocatorとgl.ylocatorで与えます。 import matplotlib.ticker as mticker gl = ax.gridlines(crs=ccrs.PlateCarree()) gl.xlocator = mticker.FixedLocator(np.arange(0, 360, 30)) # 経度線gl.ylocator = mticker.FixedLocator(np.arange(-90, 90, 30)) # 緯度線バージョン0.17の場合、cartopy_0.17_sample.pyで作図すると、図７−５− lのようになります。右端付近（東経360度付近）と上端付近（北極付近）の経度線・緯度線が描かれていないのが分かると思います。経度線を引く値np.arange(0, 360, 30)は、開始点の０から終了点の360まで30毎の値を並べることを意味していますが、終了点の値は出力されません(0.,   30., …, 330.のように最後の値が330)。バージョン0.17の場合、このような値を渡すと東経330度までしか描かれません。東経360度まで出力するには、np.arange(0, 360.1, 30)とする必要があります。 cartopy_0.17_sample2.pyでは、経度線を描く値を360.1まで、緯度線を描く値を90.1までに修正しています。バージョン0.17で動かすと、図7 − 5 − 2のように、右端付近や上端付近でも経度線・緯度線が引かれるようになります。 gl = ax.gridlines(crs=ccrs.PlateCarree()) gl.xlocator = mticker.FixedLocator(np.arange(0, 360.1, 30)) # 経度線gl.ylocator = mticker.FixedLocator(np.arange(-90, 90.1, 30)) # 緯度線
7.5.2経度線・緯度線ラベル作図時の挙動
7.1.4節で紹介した経度線・緯度線ラベルの作図では、draw_labelsオプションを使いました。cartopy_0.17_sample3.pyでは、経度線の範囲をnp.arange(- 180, 180.1, 30)で西経180度から東経180度まで30度毎に設定しました。バージョン0.17で作図したものが、図７−５−３です。 ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) # cartopy呼び出しgl.xlocator = mticker.FixedLocator(np.arange(-180, 180.1, 30)) # 経度線gl.ylocator = mticker.FixedLocator(np.arange(-90, 90.1, 30)) # 緯度線東経０度を中心とした場合にはうまくラベルが付くものの、他の範囲では失敗します。試しに東経180度を中心とした場合には、図７−５−４のように180度より東側ではラベルが付かないという不具合が起こります。ラベルを描く際に、7.1.5節、7.1.6節で紹介したax.set_xticks、ax.set_yticksを使えば、バージョン0.17でもうまくラベルをつけることができます。ここでは、経度線を30度毎、経度線ラベルを60度毎に描いた図7 – 1 – 1を0.17で描いてみます（図７−５−５）。作図に用いたcartopy_0.17_sample4.pyでは、 np.arangeに与える経度範囲を(0，360.1)に変える必要があります。なお0.18以降では、cartopy_sample11.pyのように(-180, 180)の範囲で与えていました。