tripcolor

cartopyでunstructured gridの描画を行うには,matplotlibのtripcolorを使う.Basemapではpcolorにtri=Trueオプションを付けることで描画が可能だったが,そういう簡単な方法は用意されていないようなので,matplotlibでの扱いと同じようにする.特に欠損値があると,matplotlibの場合と同様にmaskを明示的に設定することが必要だ(下の例).なお下の例では,枠がおかしなところに表示されているが,これはどうもcartopyのバグだろうと思う,このスクリプトの途中にあって今コメントアウトされているplt.gca().outline_patch.set_visible(False)を有効にすれば枠は描かれない.

import xarray as xr

import numpy.ma as ma

import numpy as np

import matplotlib.pyplot as plt

import matplotlib.tri as tri

import cartopy.crs as ccrs

fn='/data/CMIP5/lnk_nc_xyzt/ACCESS1-0/zos.hist.Omon2D/zos_Omon_ACCESS1-0_historical_r1i1p1_185001-200512.nc'

xrds=xr.open_dataset(fn)

xrdt=xrds['zos']

lats_np=np.array(xrds.coords['lat'])

lons_np=np.array(xrds.coords['lon'])

npda=np.array(xrdt.isel(time=slice(0,1)))

tsz,ysz,xsz=npda.shape

npda=npda.reshape(ysz,xsz)

plt.ion()

ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180.),frameon=False)

# plt.gca().outline_patch.set_visible(False) # This can remove the frame, which appears wrong place.

# set triangles

triang = tri.Triangulation(lons_np.ravel(),lats_np.ravel())

# set mask for triangle

npda=npda.ravel()

mask=np.zeros(triang.triangles[:,0].shape)

mask[np.isnan(npda[triang.triangles[:,0]])]=1

mask[np.isnan(npda[triang.triangles[:,1]])]=1

mask[np.isnan(npda[triang.triangles[:,2]])]=1

triang.set_mask(mask)

# set value for triangles

z=npda[triang.triangles]

z=(z[:,0]+z[:,1]+z[:,2])*0.33

# draw triangle with color

cs = ax.tripcolor(triang,z)

ax.axis([0, 360, -90, 90])

plt.title('cartopy tripcolor')

plt.savefig('fig_unstructured_cartopy_withnan')