4次元データ(T,Z,Y,X)の時間と高度を指定し、2次元断面(スカラー)を描画するためのpython (ver1)
version1: anaconda-2020.11 (CENTOS8環境)
お手軽版（colorscaleの調整ができないが、とにかく描画ができる）
main_ver1.py(51)
import sys
import os
import datetime
import inspect
import re
import argparse
import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import matplotlib as mpl
import matplotlib.pyplot as plt
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cartopy.feature as cfeature
def main():
parser = argparse.ArgumentParser(description='draw netcdf')
parser.add_argument('arg1',help='input file setting')
parser.add_argument('arg2',help='variable setting')
parser.add_argument('-t','--timenum',help='time slide from initial time')
parser.add_argument('-T','--TIME',help='YYYY-M2-D2_H2:M2:S2(UTC)')
parser.add_argument('-z','--levelnum',help='level number to be drawn')
parser.add_argument('-Z','--LEV',help='level (m) to be drawn')
parser.add_argument('-O','--OUTFIG',help='output file name')
args = parser.parse_args()
#print(args.arg1)
#print(args.arg2)
#print(args.timenum)
#print(args.TIME)
#print(args.levelnum)
#print(args.LEV)
#print(args.OUTFIG)
infile = args.arg1
if not os.path.exists(infile):
print ("infile does not exist \n",infile)
exit(1)
#open NETCDF as XARRAY
data = xr.open_dataset(infile)
# If you need, you can check all CLASS in XARRAY
#for aaa in inspect.getmembers(data):
# print(aaa)
# get variable name
varname = args.arg2
# check if the variable is included in NETCDF
if not varname in data.variables:
print("variable ", varname, " is not found")
exit(1)
# draw region setting
slon = 128.2
elon = 131.8
slat = 30.0
elat = 34.0
#############################################################
# check the dimension name for TIME
#############################################################
dimnameT = "NOTDEFINED"
for ddd in data.dims:
if ddd == "TIME":
dimnameT = "TIME"
elif ddd == "time":
dimnameT = "time"
elif ddd == "Time":
dimnameT = "Time"
if dimnameT == "NOTDEFINED":
print("TIME DIMENSION NAME is not found")
exit(1)
if args.timenum is None and args.TIME is None:
print("chose your time from following list \n")
for i, a in enumerate(data[dimnameT].values):
print(i,a)
print("enter TNUM:")
t = int(input())
TTT = data[dimnameT].values[t]
timesel={dimnameT : t}
V3D = data[varname].isel(timesel)
elif args.timenum is None:
print("YYYY-MM-DD_HH:MN:SC format \n")
ptn = r"\d\d\d\d-\d\d-\d\dT\d\d:\d\d:\d\d"
if not re.match(ptn,args.TIME):
print("TIME INPUT FORMAT WRONG")
print("Y4-M2-D2TH2:M2:S2")
print("ex 2020-07-06T03:00:00")
exit(1)
TTT = args.TIME
timesel={dimnameT : args.TIME}
V3D = data[varname].sel(timesel)
else:
t = int(args.timenum)
timesel={dimnameT : t}
TTT = data[dimnameT].values[t]
V3D = data[varname].isel(timesel)
#############################################################
# check the dimension name for Z
#############################################################
dimnameZ = "NOTDEFINED"
for ddd in data.dims:
if ddd == "LEV":
dimnameZ = "LEV"
elif ddd == "ALT":
dimnameZ = "ALT"
if dimnameZ == "NOTDEFINED":
print("Z DIMENSION NAME is not found")
exit(1)
if args.levelnum is None and args.LEV is None:
print("chose your lev from following list \n")
for k, a in enumerate(data[dimnameZ].values):
print(k,a)
print("enter ZNUM:")
k = int(input())
ZZZ = data[dimnameZ].values[k]
levsel={dimnameZ : k}
V2D = V3D.isel(levsel)
elif args.levelnum is None:
ZZZ = args.LEV
levsel={dimnameZ : ZZZ}
V2D = V3D.sel(levsel)
else:
k = int(args.levelnum)
levsel={dimnameZ : k}
ZZZ = data[dimnameZ].values[k]
V2D = V3D.isel(levsel)
# get max-min in 2dim array maxV = max(list(map(lambda x: max(x), V2D.values))) minV = min(list(map(lambda x: min(x), V2D.values))) diff = (maxV-minV)/5 diff = np.round(diff,1) maxV = maxV*1.01 minV = minV - abs(minV)*0.01 maxV = np.round(maxV,1) minV = np.round(minV,1) print(maxV,minV,diff) clevels =np.arange(minV,maxV,diff)
############################################################# # DRAW by CARTOPY ############################################################# # configure for whole figure domain xFsiz = (elon -slon) / 0.0105 yFsiz = (elat -slat) / 0.009 Fxsiz = 8 Fysiz = Fxsiz * (yFsiz/xFsiz) fig = plt.figure(figsize=(Fxsiz,Fysiz))
# configure for a single figure in the whole figure domain
ax = fig.add_subplot(111,projection=
ccrs.PlateCarree())
# draw two-dimensional distribution
V2D.plot.contourf(ax=ax,transform=ccrs.PlateCarree(),
levels=clevels,cmap='rainbow',add_colorbar=True,
extend='neither')
# set drawing lon.lat
ax.set_extent([slon,elon,slat,elat],crs=ccrs.PlateCarree())
# setting for xlabel and ylabel dlon = int(((elon-slon)/4.0)*10.0)/10.0 dlat = int(((elat-slat)/4.0)*10.0)/10.0 xticks0 = np.arange(slon,elon+0.01,dlon) xticks = np.round(xticks0,1) yticks0 = np.arange(slat,elat+0.01,dlat) yticks = np.round(yticks0,1) ax.set_xticks(xticks,crs=ccrs.PlateCarree()) ax.set_yticks(yticks,crs=ccrs.PlateCarree()) lon_formatter = LongitudeFormatter(zero_direction_label=True) lat_formatter = LatitudeFormatter() ax.xaxis.set_major_formatter(lon_formatter) ax.yaxis.set_major_formatter(lat_formatter)
# Boundary for prefecture
countries_10m = cfeature.NaturalEarthFeature('cultural',
'admin_1_states_provinces_lines',
'10m', edgecolor='black',facecolor='none')
# draw coastline
ax.coastlines()
# draw boundary of prefecture
ax.add_feature(countries_10m, linestyle=':', linewidth=0.2, alpha=0.8)
if args.OUTFIG is not None:
plt.savefig(args.OUTFIG)
# draw figure
plt.show()
if __name__ == "__main__":
main()
使い方
python main_ver1.py -h
で必要な引数などが表示される。
python main_ver1.py sample.nc U python main_ver1.py sample.nc U -t 0 -z 2 -O hoge.ps
概要
xarray でnetcdfを開き、時刻リストや高度リストが表示されるのでユーザーが選択し、
ある変数の2次元平面のデータの切り出しをxarrayで行う。描画もxarrayのインスタンスであるV2Dに対して、
xarrayのplot.contourfを使用。大変お手軽。add_coloarba=Trueにすると、自動的にカラーバーを表示。
さらに図の説明なども全部自動で書いてくれる。colorscaleさえ調整できれば完璧なのに。。。
作業用としてはこれでいいかも。図の縦横比ができるだけ１となるようにするか、縦長の図をつくるように
心がけることが大事。
V2D.plot.contourf(ax=ax,transform=ccrs.PlateCarree(),
levels=clevels,cmap='rainbow',add_colorbar=False,
extend='neither')
