Python鉛直断面GMT作図

プログラムのダウンロードはこちら

PyDraw.tar.gz(202)
サンプルのCReSSの出力(必要最低限なデータ)
samp.nc.gz(59)

使用する前に

netCDF4のpythonモジュールのインストール
GMT4.4.0を利用しています。古いです。申し訳ありません。
JboundもしくはJbound2のインストール (前坂さん作成、日本の行政区分境界の緯度経度情報提供プログラム)
(Jboundを配布していないので、必要な場合相談してください。前坂さんに聞いてみます。ダメな場合、Jboundの部分をコメントアウトしてください)

使用方法

python3 main.py samp.nc RH 0 5

main.pyを好きなように変えて、表示したい図や変数を設定してください。
main.pyは第一引数に入力ファイル名, 第2引数に表示したい変数名、第3引数に時刻のスライド数
第4引数に層番号　を記載します。

サンプルの図

解説

main.pyは以下の4つのモジュールを使って、上記の水平断面と鉛直断面を書いています。

  main.py (これは参考として書いています、みなさんの使いたいように変更してください)

使用するモジュールをインポートし、コマンドラインを引数として取得する部分です。使い勝手を変更してください。

import sys
import os
import subprocess
from read_NC import read_NC,get_unitname,get_date_utc,get_date_jst
from make_grd_GMT import make_grd_GMT
from vertical_cross_section import vertical_cross_section
from vertical_cross_section import change_input_ver_cross
from draw_GMT import draw_GMT, set_deg_per1km,draw_GMT_vertical
def main():
   if len(sys.argv) < 3:
       print ("ARGV error: infile var [t] [k]")
       exit(1)
   infile     = sys.argv[1]
   if not os.path.exists(infile):
       print ("infile does not exist \n",infile)
       exit(1)
   varname    = sys.argv[2]
   if len(sys.argv) >= 4:
       t = int(sys.argv[3])
   else :
       t = 0
   if len(sys.argv) >= 5:
       k = int(sys.argv[4])
   else :
       k = 0

NetCDFを読み込む部分です。後で、read_NCのモジュールの使い方を示します

   # read NETCDF data to get 2D dataset
   data = read_NC(infile,varname,t,k)
   lon   = read_NC(infile,'LON')
   lat   = read_NC(infile,'LAT')
   udat = read_NC(infile,"U",t,k)
   vdat = read_NC(infile,"V",t,k)
   lev            = read_NC(infile,"ALT")
   height       = lev[k]
   unit_name = get_unitname(infile,varname)
   date          = get_date_utc(infile,t)
#   date          = get_date_jst(infile,t)

読み込んだ2次元の配列データからGMTで絵を描くためのgrdファイルを作成します。
詳しくはmake_grd_GMTモジュールのところで述べます。

   # make GRD file for 2-dimensional distribution
   outgrd = varname + "_HORI.grd"
   make_grd_GMT(data,outgrd,lon,lat)
   make_grd_GMT(udat,"U_HORI.grd",lon,lat)
   make_grd_GMT(vdat,"V_HORI.grd",lon,lat)

GMTで地図を書くときに、1kmあたりの経度を求めます。

   # set 1 km -> x degrees in longitude around analysis domain
   lon_per_1km = set_deg_per1km(lat)

鉛直断面の場所を設定するところです。使い勝手を考えて、始点の緯度経度、
そこからの距離（km）、真東を0度とした反時計回りの角度を入力しています。

   # input setting for vertical cross section
   slon_ver       =   138.56
   slat_ver        =    35.3
   distance_ver =   40.0 # km
   rotate_ver     =    29.0 # deg

始点、距離、回転角度の3つの情報から、始点と終点の緯度経度座標へ変換します。

   cross_section_Start_End = change_input_ver_cross(slon_ver,slat_ver,\
                             distance_ver,rotate_ver,lon_per_1km)

鉛直断面を作成するところです。出力は、GRDファイルとオプションでテキストファイルとしても出力できます

   #  vertical cross section (data creation)
   vertical_cross_section(infile,varname,cross_section_Start_End,t,\
                          need_vertical_interp="YES",\
                          lon_per_1km=lon_per_1km,textout="NO",\
                          u="U",v="V",w="W",\
                          addvar="non",uv_along_line="UVA")

鉛直断面のGRDファイルを使って、GMTで作画するところです。詳しくはdraw_GMT_verticalのモジュールのところで説明します。

   # draw GMT for vertical cross section
   shade     = varname + "_VER.grd"
   shade_col = varname + ".cpt"
   VFig_TITLE= "vertical cross seciton of " + varname + "   T= " + date
   draw_GMT_vertical(shade=shade,shade_col=shade_col,\
                     cont="non", cont_info="-C1 -L0/10 -W4/red",
                     u="UVA_VER.grd",w="W_VER.grd",bottom=0.0,\
                     top=10.0, start=0.0,end=distance_ver,\
                     fsizX=14,fsizZ=7, PS="hoge3.ps", PNG="hoge3.png",\
                     bbase="a10:Distance\(km\):/a1:Height\(km\):neSW",\
                     title=VFig_TITLE,\
                     unit=unit_name, vec_INTV="2/0.5",\
                     vec_info="0.03c/0.11c/0.1c -W3 -G0",\
                     vec_per1cm=48,wamp=5.0,\
                     PNG_yohakuX=5.0, PNG_yohakuY=3.5,PNG_DPI=100.0)

make_grd_GMTモジュールで作成した2次元断面のGRDファイルを使ってGMTで作画するところです。詳しくはdraw_GMTモジュールのところで説明します。

   # draw GMT for horizontal distribution
   HFig_TITLE= varname + " at " + str(height) + " m" + "    T=  " + date
   shade     = varname + "_HORI.grd"
   shade_col = varname + ".cpt"
   draw_GMT(shade=shade, 
            u="U_HORI.grd",
            v="V_HORI.grd",
            shade_col=shade_col, 
            cont_info="-C1 -L0/10 -W4/blue",
            slon=138.42,
            elon=139.8,
            slat=34.7,
            lon_per_1km=lon_per_1km,
            fsizX=14,
            fsizY=13,
            PS="hoge2.ps",
            PNG="hoge2.png",
            bbase="a0.6/0.2neSW",
            JboundArea="tokai kanto tohoku koushinetsu",
            vec_INTV=0.04,
            vec_info="0.03c/0.11c/0.1c -W3 -G0",
            vec_per1cm=48,
            title=HFig_TITLE,
            unit=unit_name,
            obj_draw_data=cross_section_Start_End,
            obj_draw_command="psxy -W8/red -N",
            PNG_yohakuX=5.0, 
            PNG_yohakuY=2.0,
            PNG_DPI=100.0)

2つのPNGを縦方向にくっつけて1つのPNGに作成するところです。結果は上記に示した図となります。

   #  Merging multiple FIgures as One ##########
   command = "convert -append" + " hoge2.png hoge3.png hoge.png"
   subprocess.run(command,shell=True)
####################################

main関数を呼び出すところです。

if __name__ == "__main__":
   main()

read_NCモジュール

NetCDFの読み込みについて、私が使っているバリエーションの範囲で
自動的に包括的にデータを読み込めるように工夫をしたものです。次元変数名を
LON,LAT,LEV,TIME を仮定しています。必要であれば、適宜変更してください。

必要とする標準モジュールです。

import netCDF4
import inspect
import datetime

ファイルの中から、時刻tのタイムステップのUTC時刻を文字列として返すサブルーチン

def get_date_utc(infile,t):
   nc = netCDF4.Dataset(infile,'r')
   unixtime = nc.variables["TIME"][t]
   nc.close()
   dd = datetime.datetime.fromtimestamp(int(unixtime),datetime.timezone.utc)
   hoge = '%04d-%02d-%02d_%02d:%02d:%02d' % (dd.year,dd.month,dd.day,dd.hour,\
          dd.minute,dd.second)
   return hoge + " UTC"

ファイルの中から、時刻tのタイムステップのJST時刻を文字列として返すサブルーチン

def get_date_jst(infile,t):
   nc = netCDF4.Dataset(infile,'r')
   unixtime = nc.variables["TIME"][t]
   nc.close()
   dd = datetime.datetime.fromtimestamp(int(unixtime))
   hoge = '%04d-%02d-%02d_%02d:%02d:%02d' % (dd.year,dd.month,dd.day,dd.hour,\
          dd.minute,dd.second)
   return hoge + " JST"

ファイルの中から、ある変数に含まれるmissing_valueを読み込み、返すサブルーチン

def get_missing_value(infile,varname):
   nc = netCDF4.Dataset(infile,'r')
   if not varname in nc.variables:
       print ("not found ", varname)
       exit(1)
   missing = nc.variables[varname].missing_value
   nc.close()
   return missing

ファイルの中から、ある変数の単位を読み込み、文字列として返すサブルーチン

def get_unitname(infile,varname):
   nc = netCDF4.Dataset(infile,'r')
   if not varname in nc.variables:
       print ("not found ", varname)
       exit(1)
   unitname = nc.variables[varname].units
   nc.close()
   return unitname

ファイルの中から、ある変数が持つ次元変数の組み合わせを調べるサブルーチン。
TIME,LEV,LAT.LON の4つの次元のありなしを4つの要素をもつ1次元配列として返す。
TIME,LAT,LONの3次元の場合、[1,0,1,1]を返す。LEV,LAT,LONの3次元の場合は[0,1,1,1]を返す。

def check_NC_type(infile,varname):
   nc = netCDF4.Dataset(infile,'r')   
   if not varname in nc.variables:
       print ("not found ", varname)
       exit(1)       
   dim = nc.variables[varname].dimensions
   nc.close()
   dim_check = [0,0,0,0]
   if "TIME" in dim:
       dim_check[0] = 1
   else:
       dim_check[0] = 0
   if "LEV" in dim or "ALT" in dim:
       dim_check[1] = 1
   else:
       dim_check[1] = 0
   if "LAT" in dim:
       dim_check[2] = 1
   else:
       dim_check[2] = 0
   if "LON" in dim:
       dim_check[3] = 1
   else:
       dim_check[3] = 0
   return dim_check

2次元平面で時間変化する変数を読み込むサブルーチン。
ただし、tを指定しない場合、すべての時間を読み込む

def read_NC_T2D(infile,varname,t=-1):
   nc = netCDF4.Dataset(infile,'r')
   if not varname in nc.variables:
       print ("not found ", varname)
       exit(1)
   if t!= -1:
       DATA = nc.variables[varname][t][:][:]
   else:
       DATA = nc.variables[varname][:][:][:]
   nc.close()
   return DATA

3次元空間が時間変化する変数を読み込むサブルーチン。
ただし、tを指定しない場合、すべての時間を読み込む

def read_NC_T3D(infile,varname,t,k=-1):
   nc = netCDF4.Dataset(infile,'r')
   if not varname in nc.variables:
       print ("not found ", varname)
       exit(1)
   if k!=-1:
       DATA = nc.variables[varname][t][k][:][:]
   else:
       DATA = nc.variables[varname][t][:][:][:]
   nc.close()
   return DATA

2次元空間が時間変化しない変数を読み込むサブルーチン。

def read_NC_2D(infile,varname):
   nc = netCDF4.Dataset(infile,'r')
   if not varname in nc.variables:
       print ("not found ", varname)
       exit(1)
   DATA2D = nc.variables[varname][:][:]
   nc.close()
   return DATA2D

1次元変数を読み込むサブルーチン。

def read_NC_1D(infile,varname):
   nc = netCDF4.Dataset(infile,'r')
   if not varname in nc.variables:
       print ("not found ", varname)
       exit(1)
   DATA1D = nc.variables[varname][:]
   nc.close()
   return DATA1D

NetCDFを読み込むmainサブルーチン

def read_NC(infile, varname,t=-1,k=-1):
   !その変数がどんな次元変数を持つか調べる
   type = check_NC_type(infile,varname)
   ndim = sum(type)
   if type   == [1,0,1,1] :
       print (varname,":2D evolution -> 2D snapshot")
       DATA = read_NC_T2D(infile,varname,t)
   elif type == [1,1,1,1] :
       if(k!=-1):
           print (varname,":3D evolution -> 2D snapshot")
           DATA = read_NC_T3D(infile,varname,t,k)
       else:
           print (varname,":3D evolution -> 3D snapshot")
           DATA = read_NC_T3D(infile,varname,t)
   elif type == [0,1,1,1] :
       if(k!=-1):
           print (varname,":3D static -> 2D snapshot")
           DATA = read_NC_T2D(infile,varname,k)
       else:
           print (varname,":3D static -> 3D snapshot")
           DATA = read_NC_T2D(infile,varname)
   elif type == [0,0,1,1] :
       print (varname,":2D static")
       DATA = read_NC_2D(infile,varname)
   elif ndim == 1:
       print (varname,":1dimensional")
       DATA = read_NC_1D(infile,varname)
   return DATA

make_grd_GMT モジュール

vertical_cross_section モジュール

draw_GMT モジュール