気象×Python　〜3次元描画〜

気象データは地上だけじゃなくて、何層にも等圧面データがあるんで3次元的に描けるか試してみました。とりあえず自分の備忘録として、ゆるく書きます。

3次元の降水データはこんなん[(JAXAより引用してます。)](http://www.sapc.jaxa.jp/use/rd/3d_precipitation/)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/354874/f8a2dd93-9e14-6678-f160-498f27b6ec27.png)
#1.　使用データ
▶[京都大学生存圏研究所（RISH: Research Institute for Sustainable Humanosphere）](http://database.rish.kyoto-u.ac.jp/arch/jmadata/gpv-netcdf.html)から取得したメソ数値予報モデルGPV（MSM）
▶NetCDF形式
16層からなる（1000,975,950,925,900,850,800,700,600,500,400,300,250,200,150,100）
ちなみに2014/12/17は雪がめっちゃ降った日。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/354874/bd911342-d3c5-17b9-e340-ca26a0ddc12c.png)
#2.　コード

```temp_3d.py
# -*- coding: utf-8 -*-
import matplotlib
matplotlib.rcParams['backend'] = 'TkAgg'
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import netCDF4
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.mplot3d import Axes3D

nc_p = netCDF4.Dataset('20141217_p.nc', 'r')
temp = nc_p['temp'][0][:][:][:] -273.15
p = nc_p['p'][:]
lat_p = nc_p['lat'][:]
lon_p = nc_p['lon'][:]

fig = plt.figure(facecolor='white',dpi=300)
ax = fig.gca(projection='3d')
m = Basemap(projection="cyl", resolution="i", llcrnrlat=30,urcrnrlat=45, llcrnrlon=130, urcrnrlon=145,fix_aspect=False, ax=ax)
ax.add_collection3d(m.drawcoastlines(linewidth=1))
ax.add_collection3d(m.drawcountries(linewidth=1))
ax.view_init(azim=290, elev=35)

meridians = np.arange(130, 145 + 5, 5)
parallels = np.arange(30, 45 + 5, 5)
pressure = np.arange(100, 1000 + 100 , 100)
ax.set_xticks(meridians)
ax.set_xticklabels(meridians,fontsize=5)
ax.set_yticks(parallels)
ax.set_yticklabels(parallels,fontsize=5)
ax.set_zticks(pressure)
ax.set_zticklabels(pressure, fontsize=5)
ax.set_zlabel('hPa', labelpad=0.5,fontsize=5)
ax.invert_zaxis()

X, Y =  np.meshgrid(lon_p , lat_p)
level = np.arange(-45,20,5)
#試しに2層で。
#1000,975,950,925,900,850,800,700,600,500,400,300,250,200,150,100
m.contourf(X,Y, temp[5,:,:], offset=850, cmap=cm.jet , alpha=0.7) #850hPa
img = m.contourf(X, Y, temp[9,:,:], offset=500, cmap=cm.jet , alpha=0.7) #500hPa
# (left,bottom,right,top)
cax = fig.add_axes([0.22, 0.05, 0.65, 0.045]) 
fig.colorbar(img,cax, orientation="horizontal")
```
日本地図の位置が事案
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/354874/4e27bf14-b8ce-81ff-2d74-af3a3be8b997.png)
ちょいと修正。

```new_temp_3d.py
#途中から
meridians = np.arange(130, 145 + 5, 5)
parallels = np.arange(30, 45 + 5, 5)
pressure = np.arange(100, 1000 + 100 , 100)
ax.set_xticks(meridians)
ax.set_xticklabels(meridians,fontsize=5)
ax.set_yticks(parallels)
ax.set_yticklabels(parallels,fontsize=5)
ax.set_zticks(pressure)
ax.set_zticklabels(pressure[::-1], fontsize=5)
ax.set_zlabel('hPa', labelpad=0.5,fontsize=5)

X, Y =  np.meshgrid(lon_p , lat_p)
level = np.arange(-45,20,5)
m.contourf(X,Y, temp[5,:,:], offset=1100-850, cmap=cm.jet , alpha=0.7) #850hPa
img = m.contourf(X, Y, temp[9,:,:], offset=1100-500, cmap=cm.jet , alpha=0.7) #500hPa
cax = fig.add_axes([0.22, 0.05, 0.65, 0.045]) 
fig.colorbar(img,cax, orientation="horizontal")
```
正しい位置には戻ったが、正直微妙。
もっとそれっぽい可視化方法ありそうだけど、、誰か教えてくださいまし。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/354874/631d3838-d3d2-f882-1ff4-efef6347a779.png)
以上！！！