Xarrayで断面図を作成するときに海底地形を描く方法について記載します。
海底地形のデータをDatasetに作成
ny,nx=grd_s.vgrid.h[1:-1,1:-1].shape
ary0=np.zeros((2,ny,nx))
ary0[0,:,:]=-6000.0
ary0[1,:,:]=-grd_s.vgrid.h[1:-1,1:-1]
ary1=np.zeros((2,ny,nx))
ary1[0,:,:]=1.0
ary1[1,:,:]=0.0
dsA=xr.Dataset()
vlist=["z_rho","z_w","z_v","z_u"]
for vname in vlist:
dsA[vname]=ds0a[vname]
dsA["btm_pos"]=xr.DataArray(ary0,dims=["bh","yh","xh"])
dsA["btm_fil"]=xr.DataArray(ary1,dims=["bh","yh","xh"])
vlist=vlist+["btm_pos"]
dsA=dsA.set_coords(vlist)
#dsA.to_zarr("{0}/test".format(datdir),mode="a")- grd_sというオブジェクトのhという変数に海底の深度情報がはいっているものとします。
- ary0として、全地点の最大水深のデータ（実際の最大水深より深めにとっておいて描画で切る方が良い）と、水深データの配列を作ります。
- ary1として、塗り分けをするための0,1配列を作ります。
- ary0を”btm_pos”、ary1を”btm_fil”という変数名でDatasetに格納します。
- “btm_pos”の方を非次元座標として保存します。
- 必要があればZarrに保存します。
from matplotlib.colors import ListedColormap, BoundaryNorm
thres=1.e10
colors=["white","grey"]
cmap=ListedColormap(colors)
norm=BoundaryNorm([-thres,0.0,thres],cmap.N)fig,ax=plt.subplots(1,1)
ds0d["btm_fil"].isel(xh=600).plot(x="lat_rho",y="btm_pos",cmap=cmap,norm=norm,add_colorbar=False)
ds0d["SiO3_vadv"].isel(xh=600,time=0).plot(x="lat_rho",y="z_rho",robust=True)
ax.set_ylim([-5500,0.0])
ax.set_xlim([10,20])
解説
- 海底を塗るためのカラーマップを作成します。上記の例では海底をgray、その他を
- 先に海底地形を描画します。”btm_fil”の変数を描画し、鉛直座標として”btm_pos”を指定します。
- その後、変数を描画します。
