RuntaScience diary

気象系データを扱う学生 旅が好きです

Welcome to my blog

About

Cartopy-Pythonを用いた、Merra-2の利用

こんにちは!

今回は、NASA(アメリカ航空宇宙局, National Aeronautics and Space Administration)のGES DISCで無料で入手可能なnetCDFのデータをCartopyを用いて表示させたいと思います。

netCDFの扱い方

netCDFとは何だ??という方は、以下を参照ください。
netCDF データ格納形式の基礎—ヘルプ | ArcGIS for Desktop

データのダウンロード

GES DISCからのデータのダウンロードの仕方は、ホームページに描いてあるので割愛します。

disc.gsfc.nasa.gov

pythonでデータの中身を見る

初めに、必要なものをインポートします

from netCDF4 import Dataset #netCDFを扱うため
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt #グラフ描写のため
import numpy as np #数値計算用

私は、1ヶ月あたりの気温の月平均値のみを使いました。2019年12月のデータです。

nc = Dataset('MERRA2_400.tavgM_2d_slv_Nx.201912.nc4.nc4', mode='r')
nc

表示させるとたくさん情報が出てきます

<class 'netCDF4._netCDF4.Dataset'>
・・・
SouthernmostLatitude: -90.0
NorthernmostLatitude: 90.0
WesternmostLongitude: -180.0
EasternmostLongitude: 179.375
LatitudeResolution: 0.5
LongitudeResolution: 0.625
DataResolution: 0.5 x 0.625
・・・
https://goldsmr4.gesdisc.eosdis.nasa.gov:443/opendap/MERRA2_MONTHLY/M2TMNXSLV.5.12.4/2019/MERRA2_400.tavgM_2d_slv_Nx.201912.nc4.nc4?T2M,time,lat,lon   
dimensions(sizes): time(1), lat(361), lon(576)   
variables(dimensions): float32 T2M(time,lat,lon), float64 lat(lat), float64 lon(lon), int32 time(time)  

内容をしっかり確認します
T2M:2 m地点の気温。(時間,緯度,経度)で格納されている。
lat :緯度。(緯度)で格納されている。
lon:経度。(経度)で格納されている。

それでは、データを変数に入れていきます。

lon = nc.variables["lon"][:].data
lat = nc.variables["lat"][:].data
t2m0= nc.variables["T2M"][:,:,:].data

.dataを設定すると、maskedarrayというのがなくなり、普通のarrayになります。
それぞれの形も確認しておきましょう

print(t2m0.shape)
print(lat.shape)
print(lon.shape)

今回は、月平均値を使ったのでtimeの次元が1でしたが、日平均や月平均値をとる場合は以下の操作を施します。
t2m.shapeをすると、(時間,緯度,経度)の配列数がでるので、それぞれが(0,1,2)番目であることから、
axis(軸)=0(時間軸)の平均をとります

t2m = t2m0.mean(axis=0)

単位の変換(K⇒℃)

t2mはケルビンなので℃に変換します

temp = t2m -273.15

グラフの作成

それではグラフを作成していきます

Cartopyを用いた基礎的な内容は、以下を見てください

runtascience.hatenablog.com

#緯度経度グリッド作成のため
import matplotlib.ticker as mticker
#緯度経度線を度数表示にするため
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

グラフ-global

def main():
    fig = plt.figure(figsize=(10,5))
    plt.rcParams["font.size"] = 18

    ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
    ax.set_title("Temperature Dec. 2019")
    ax.coastlines(resolution="110m") 
    ax.add_feature(cfeature.BORDERS, linestyle=":")
    

    #緯度・経度線のグリッドの設定(参考*2)
    gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                  linewidth=2, color='gray', alpha=0.5, linestyle=':')
    gl.xlabels_top = False
    gl.ylabels_right = False
    dlon, dlat = 30,30
    xticks = np.arange(-180, 180.1, dlon)
    yticks = np.arange(-90, 90.1, dlat)
    gl.xlocator = mticker.FixedLocator(xticks)
    gl.ylocator = mticker.FixedLocator(yticks)
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    gl.xlabel_style = {'size': 12}
    gl.ylabel_style = {'size': 12}

    
    #データ描写(参考*1)
    clevs = np.arange(-50,50,5)
    im = ax.contourf(lon, lat, temp,clevs,cmap="jet")

    #カラーバーの設定
    cb = plt.colorbar(im, orientation="vertical", pad=0.02, aspect=25, shrink=0.85)
    cb.set_label("Temperature ($^\circ$C)",rotation=90,labelpad=5)
    cb.ax.tick_params(size=5)

    plt.show()

    fig.savefig("XXX.png", format="png", dpi=330)

if __name__ == "__main__":
    main()

f:id:RuntaScience:20200517201755p:plain

グラフ-japan

日本周辺も作成しました

def main():
    fig = plt.figure(figsize=(10,10))
    plt.rcParams["font.size"] = 18

    ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
    ax.set_title("Temperature Dec. 2019")
    ax.coastlines(resolution="10m") 
    ax.add_feature(cfeature.BORDERS, linestyle=":")
    
   #緯度・経度線のグリッドの設定(参考*2)
    gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                  linewidth=2, color='gray', alpha=0.5, linestyle=':')
    gl.xlabels_top = False
    gl.ylabels_right = False
    dlon, dlat = 5, 5
    xticks = np.arange(-180, 180.1, dlon)
    yticks = np.arange(-90, 90.1, dlat)
    gl.xlocator = mticker.FixedLocator(xticks)
    gl.ylocator = mticker.FixedLocator(yticks)
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    gl.xlabel_style = {'size': 18}
    gl.ylabel_style = {'size': 18}
   
    
    #描写(参考*1)
    clevs = np.arange(-50,50,5)
    im = ax.contourf(lon, lat, temp,clevs,cmap="jet")
    ax.set_extent([128, 148, 30, 50], crs=ccrs.PlateCarree())

    #カラーバーの設定
    cb = plt.colorbar(im, orientation="vertical", pad=0.02, aspect=25, shrink=0.8)
    cb.set_label("Temperature ($^\circ$C)",rotation=90,labelpad=5)
    cb.ax.tick_params(size=5)

    plt.show()

    fig.savefig("XXX.png", format="png", dpi=330)

if __name__ == "__main__":
    main()

f:id:RuntaScience:20200517202057p:plain

Keys

1 glidlineの種類: "-", "--", "-.", ":"など
matplotlib.axes.Axes.grid — Matplotlib 3.2.1 documentation

2 colorbarの種類: "jet", "rainbow"などたくさんある。"jet_r"のように_rをつけると、色が反対になる
Choosing Colormaps in Matplotlib — Matplotlib 3.2.1 documentation

3 colorbarの設定

名前 内容
orientation vertical=縦方向, horizontal=横方向
pad グラフとカラーバーの間隔
aspect カラーバーの幅
shrink figに対するカラーバーの縦の長さ
labelpad カラーバーラベルの文字の回転
shrink カラーバーとラベルの幅    

matplotlib.pyplot.colorbar — Matplotlib 3.1.2 documentation

4 単位とLaTeX表示

runtascience.hatenablog.com

参考コード

*1 Merra-2のデータの使い方
GES DISC
*2 緯度経度のグリッド設定
Cartopy map gridlines and tick labels — cartopy 0.13.0 documentation

それでは 🌏