※このページは、かつて isotope.iis.u-tokyo.ac.jp/~kanon に公開していた内容を転写したものです。
ここではPython3のgtool3
パッケージ( https://github.com/hyungjun/gtool3 )を使用してgtoolファイルを読み出す方法を解説しています。このほかに、より処理速度の速いxgtool3
( https://github.com/k1bb/xgtool3/tree/main ) もありますので、用途に応じて使い分けてもよいかもしれません。 xgtool3については、以下が参考になります。
- シンプルな解説:https://qiita.com/wm-ytakano/items/ddeecf2b43a4d1b200f2#gtool3
- quick plot パッケージ(プロトタイプ):https://github.com/kanonundgigue/plot_xgt3
1.導入 / Introduction
こちらのioモジュールを利用させていただきます。
pip install git+https://github.com/hyungjun/gtool3
でインストールできます。
anaconda(miniconda)でパッケージを管理している場合も、インストールしたい環境で上記のpipコマンドを実行して下さい。
We will be using the io module from here.
You can install it with pip install git+https://github.com/hyungjun/gtool3
.
If you are using anaconda (miniconda) to manage packages, you can also run the above pip command in the environment you want to install.
なお、PyPIに置かれているバージョンは古いため、pip install gtool3
でインストールすると軸情報ファイルを読み込むことができません。
Note that the version in PyPI is old, so the axis information file cannot be read if you install it with pip install gtool3
.
早速、モジュールを読み込みましょう。
Let's load the module.
from gtool3 import gtopen as open_gtool
データの読み込みは、
Data can be loaded using:
gt = open_gtool("q")
表示してみると
Try to display:
gt.vars
OrderedDict([('Q', Q, (12, 40, 64, 128) : >f4)])
順序付き辞書 OrderedDictとなっていますね。
It's an ordered dictionary OrderedDict.
ということは、
-
keys()
: 各要素のキーkeyを取り出す -
values()
: 各要素の値valueを取り出す -
items()
: 各要素のキーkeyと値valueを取り出す
を利用できます。
This means that
-
keys()
: retrieve the key of each element. -
values()
: retrieve the value value of each element. -
items()
: retrieve the key and value of each element.
are available.
gt.vars.keys()
odict_keys(['Q'])
gt.vars.values()
odict_values([Q, (12, 40, 64, 128) : >f4])
gt.vars.items()
odict_items([('Q', Q, (12, 40, 64, 128) : >f4)])
たとえば、以下のように値を取り出せます。
For example, you can retrieve the value as follows
print(gt.vars["Q"][:,0,0,0])
[0.00018083 0.00050614 0.00030294 0.00047958 0.00166131 0.00357095
0.00362982 0.00361975 0.00229563 0.00135356 0.00053724 0.00028788]
2.ヘッダー情報を得る / Reading header information
gtoolコマンドだとngthead
で表示してくれますが、pythonだと?
The gtool command will show it in ngthead
, but in python?
gt.vars["Q"].header
[00] IDFM : 9010: DSET :jra25.gndg : ITEM :Q :
[03] EDIT1 : : EDIT2 : : EDIT3 : :
[06] EDIT4 : : EDIT5 : : EDIT6 : :
[09] EDIT7 : : EDIT8 : : FNUM : 1:
[12] DNUM : 1: TITL1 :water vapor spec: TITL2 :ific humidity :
[15] UNIT :kg/kg : ETTL1 : : ETTL2 : :
[18] ETTL3 : : ETTL4 : : ETTL5 : :
[21] ETTL6 : : ETTL7 : : ETTL8 : :
[24] TIME : ** NOTE ** : UTIM :HOUR : DATE : ** NOTE ** :
[27] TDUR : ** NOTE ** : AITM1 :GLON128 : ASTR1 : 1:
[30] AEND1 : 128: AITM2 :GGLA64 : ASTR2 : 1:
[33] AEND2 : 64: AITM3 :HETA40 : ASTR3 : 1:
[36] AEND3 : 40: DFMT :UR4 : MISS : -0.9990000E+03:
[39] DMIN : -0.9990000E+03: DMAX : -0.9990000E+03: DIVS : -0.9990000E+03:
[42] DIVL : -0.9990000E+03: STYP : 1: COPTN : :
[45] IOPTN : 0: ROPTN : 0.0000000E+00: DATE1 : ** NOTE ** :
[48] DATE2 : ** NOTE ** : MEMO1 :isotope3:/home/k: MEMO2 :anon/miroc5-iso/:
[51] MEMO3 :out/jra25.gndg/y: MEMO4 :1981 : MEMO5 : :
[54] MEMO6 : : MEMO7 : : MEMO8 : :
[57] MEMO9 : : MEMO10 : : CDATE : ** NOTE ** :
[60] CSIGN :MIROC : MDATE :20230323 192039 : MSIGN :kanon :
[63] SIZE : 327680:
** NOTE **
[24] TIME :[17365476 ... 17373492], (12)
[26] DATE :[19810116 120000 ... 19811216 120000], (12)
[27] TDUR :[744 ... 744], (12)
[47] DATE1 :[19810101 000000 ... 19811201 000000], (12)
[48] DATE2 :[19810201 000000 ... 19820101 000000], (12)
[59] CDATE :[20201223 181831 ... 20201223 203313], (12)
アイテム名を取得するには、
To get the item name:
gt.vars["Q"].header["ITEM"]
'Q'
3.簡単な絵を描いてみる / Drawing a simple figure
3.1.軸情報を取得 / Getting axes information
GTOOL形式のファイルは、headerに軸名が書いてあるだけで、自分で情報を持っていません。
GTOOL format files only have the axis name in the header, and do not have any information of their own.
軸名を確認してみましょう。
Let's check the names of axes.
xname = gt.vars["Q"].header["AITM1"] # x-axis
yname = gt.vars["Q"].header["AITM2"] # y-axis
zname = gt.vars["Q"].header["AITM3"] # z-axis
print("x-axis :",xname)
print("y-axis :",yname)
print("z-axis :",zname)
x-axis : GLON128
y-axis : GGLA64
z-axis : HETA40
このような軸の情報ファイルを、別途読み込んでくる必要があります。
You will need to load the information file for such an axis separately.
一般に、軸情報ファイルの格納先は、GTOOLがインストールされている環境であれば環境変数GTAXDIR
に格納されているため
In general, the axis information file is stored in the environment variable GTAXDIR
if GTOOL is installed in your environment.
def cmd(cmd):
# Pythonでシェルコマンドを動かし、結果を取得する関数
# 参考 : https://qiita.com/inatatsu_csg/items/40b11701d256a84a0510
import subprocess
process = (subprocess.Popen(cmd, stdout=subprocess.PIPE, shell=True).communicate()[0]).decode('utf-8')[:-1]
return process
gtaxdir = cmd("echo $GTAXDIR")
print(gtaxdir)
/home/kanon/gtool3-axis
とすることで、確認できます。
もしGTAXDIR
が空で軸情報ファイルが見当たらない場合は、こちらからダウンロードできます。
You can check it by using the following command.
If GTAXDIR
is empty and you cannot find the axis information file, you can download it from here.
# Stores axes information
dlon = open_gtool(gtaxdir+"/GTAXLOC."+xname).vars[xname][0,0,0,:]
dlat = open_gtool(gtaxdir+"/GTAXLOC."+yname).vars[yname][0,0,0,:]
dlev = open_gtool(gtaxdir+"/GTAXLOC."+zname).vars[zname][0,0,0,:]
# number of grids
nlon =len(dlon)
nlat = len(dlat)
nlev = len(dlev)
/home/kanon/.conda/envs/py2021/lib/python3.7/site-packages/gtool3/format.py:71: UserWarning: cannot cast IOPTN:[b' '] to <class 'int'>, replace with empty string.
k, values, self.fmt[k][0]))
/home/kanon/.conda/envs/py2021/lib/python3.7/site-packages/gtool3/format.py:71: UserWarning: cannot cast ROPTN:[b' '] to <class 'float'>, replace with empty string.
k, values, self.fmt[k][0]))
以上のWarningは仕様なので、気にしなくて大丈夫です。
The above warnings are specifications, so don't worry about them.
print("x-axis (grid number="+str(nlon)+"):\n",dlon)
print("y-axis (grid number="+str(nlat)+"):\n",dlat)
print("z-axis (grid number="+str(nlev)+"):\n",dlev)
x-axis (grid number=129):
[ 0. 2.8125 5.625 8.4375 11.25 14.0625 16.875 19.6875
22.5 25.3125 28.125 30.9375 33.75 36.5625 39.375 42.1875
45. 47.8125 50.625 53.4375 56.25 59.0625 61.875 64.6875
67.5 70.3125 73.125 75.9375 78.75 81.5625 84.375 87.1875
90. 92.8125 95.625 98.4375 101.25 104.0625 106.875 109.6875
112.5 115.3125 118.125 120.9375 123.75 126.5625 129.375 132.1875
135. 137.8125 140.625 143.4375 146.25 149.0625 151.875 154.6875
157.5 160.3125 163.125 165.9375 168.75 171.5625 174.375 177.1875
180. 182.8125 185.625 188.4375 191.25 194.0625 196.875 199.6875
202.5 205.3125 208.125 210.9375 213.75 216.5625 219.375 222.1875
225. 227.8125 230.625 233.4375 236.25 239.0625 241.875 244.6875
247.5 250.3125 253.125 255.9375 258.75 261.5625 264.375 267.1875
270. 272.8125 275.625 278.4375 281.25 284.0625 286.875 289.6875
292.5 295.3125 298.125 300.9375 303.75 306.5625 309.375 312.1875
315. 317.8125 320.625 323.4375 326.25 329.0625 331.875 334.6875
337.5 340.3125 343.125 345.9375 348.75 351.5625 354.375 357.1875
360. ]
y-axis (grid number=64):
[ 87.8638 85.09653 82.31291 79.525604 76.7369 73.94752
71.15775 68.36776 65.57761 62.787354 59.99702 57.20663
54.4162 51.625732 48.83524 46.044727 43.254196 40.46365
37.673088 34.882523 32.091946 29.30136 26.510769 23.720175
20.929575 18.138971 15.348365 12.557756 9.767145 6.9765334
4.1859207 1.395307 -1.395307 -4.1859207 -6.9765334 -9.767145
-12.557756 -15.348365 -18.138971 -20.929575 -23.720175 -26.510769
-29.30136 -32.091946 -34.882523 -37.673088 -40.46365 -43.254196
-46.044727 -48.83524 -51.625732 -54.4162 -57.20663 -59.99702
-62.787354 -65.57761 -68.36776 -71.15775 -73.94752 -76.7369
-79.525604 -82.31291 -85.09653 -87.8638 ]
z-axis (grid number=40):
[0.9974992 0.99149853 0.98299694 0.9719956 0.958493 0.9419898
0.92248577 0.90048254 0.87597704 0.8489725 0.8199673 0.786951
0.74692285 0.69888484 0.6428333 0.57447207 0.50147766 0.43644983
0.37985232 0.33059552 0.28772512 0.25041142 0.21794017 0.18968128
0.16508035 0.14367269 0.12504333 0.10882763 0.09471573 0.08243426
0.07174488 0.06244121 0.05427773 0.04693731 0.04015774 0.03361269
0.02659318 0.01861324 0.01060254 0.00290465]
緯度方向のみ、本来128グリッドのところ、1周つながるように0と360が重複して入っていますね。
モデルの出力結果の方は、128個のデータしか持っていないので、気をつけましょう。
Only in the latitudinal direction, the original 128 grids are duplicated with 0 and 360 grids to make a full circle.
The output result of the model has only 128 data, so be careful.
3.2.描画してみる / Drawing
今回は、[JAMSTEC山下さんのTips (GitHub)]https://yyousuke.github.io/matplotlib/)
にならい、cartopyとmatplotlibを使ってみます。
This time, let's use cartopy and matplotlib following JAMSTEC Dr.Yamashita's Tips (GitHub).
# Reading modules
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
# Read a variable
gt = open_gtool("T2")
var = gt.vars["T2"][:,0,:,:] - 273.15
そのまま描画してしまうと、3.1で指摘した通り、x方向のデータの数が軸情報ファイルと描画したい変数で異なるので、次の通りエラーになります。
If you draw it as it is, as pointed out in 3.1, the number of data in the x direction differs between the axis information file and the variable you want to draw, so you will get an error as follows.
# Drawing
fig = plt.figure() # Setting a plotting region(matplotlib)
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree()) # Using Cartopy
ax.coastlines() # Drawing coastlines
ax.contourf(dlon,dlat,var[0,:,:],transform=ccrs.PlateCarree()) # Drawing the variable
# Expected errors
# ---------------------------------------------------------------------------
#
# TypeError Traceback (most recent call last)
# ~~~
# TypeError: Length of x (129) must match number of columns in z (128)
dlon
のいちばん端だけ落としてみましょう。
Let's drop just the very end of dlon
.
# Drawing
fig = plt.figure() # Setting a plotting region(matplotlib)
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree()) # Using Cartopy
ax.coastlines() # Drawing coastlines
ax.contourf(dlon[:-1],dlat,var[0,:,:],transform=ccrs.PlateCarree()) # Drawing the variable
描けた!しかし、経度0度がきれてしまいました。
今度は、変数の経度360度のところに経度0度の情報を代入してみましょう。
Succeeded! However, longitude 0 degrees have been run out. Now, let's substitute the information of longitude 0 degrees into the variable longitude 360 degrees.
import numpy as np
ntime = len(var[:,0,0]) # number of time steps
var_cyc = np.zeros((ntime,nlat,nlon))
var_cyc[:,:,:-1] = var[:,:,:]
var_cyc[:,:,-1] = var[:,:,0]
# Drawing
fig = plt.figure() # Setting a plotting region(matplotlib)
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree()) # Using Cartopy
ax.coastlines() # Drawing coastlines
cs = ax.contourf(dlon,dlat,var_cyc[0,:,:],transform=ccrs.PlateCarree()) # Drawing the variable
cbar = plt.colorbar(cs, orientation="horizontal") # Drawing a color bar
ax.set_title("2m Temperature Jan. [℃]") # Adding a title
plt.savefig("t2jan.png") # Save the figure
4.一連の処理をモジュール化する / Modularize the processes
ここでは、gtool形式のファイルを読みこんでから、描画する前までの処理をモジュール化します。
Here, we modularize the process from reading a gtool format file to before drawing.
from gtool3 import gtopen as open_gtool
import numpy as np
import re
import warnings
warnings.simplefilter("ignore")
def cmd(cmd):
# Pythonでシェルコマンドを動かし、結果を取得する関数
# 参考 : https://qiita.com/inatatsu_csg/items/40b11701d256a84a0510
import subprocess
process = (subprocess.Popen(cmd, stdout=subprocess.PIPE, shell=True).communicate()[0]).decode('utf-8')[:-1]
return process
class ReadGtool:
def __init__(self):
self.var = "" # set input file path
self.dim = "xyzt" # set dimension: set xyt to remove z-axis
self.xcyclic = True # if connect the longitudes = 0 to 360
self.xnum = None # x-levels to select, if None, all levels are included
self.ynum = None # y-levels to select, if None, all levels are included
self.znum = None # z-levels to select, if None, all levels are included
self.tnum = None # t-levels to select, if None, all levels are included
self.axes_map = {"t": 0, "z": 1, "y": 2, "x": 3}
self.undef = -999.0 # default undef value in MIROC
def get_length(self, obj):
try:
length = len(obj)
except TypeError:
length = 0 # 任意のデフォルト値を設定するか、適切な処理を行います。
return length
def list2num(self, obj):
if isinstance(obj, list):
for item in obj:
if not isinstance(item, (int, float)):
print(f"Non-numeric item found: {type(item).__name__}")
return
return obj[0]
else:
if not isinstance(obj, (int, float)):
print(f"Non-numeric item found: {type(obj).__name__}")
return
return obj
def get_varname(self):
varn = list(open_gtool(self.var).vars.keys())[0]
return varn
def get_gtool_xyzaxis(self):
filename = self.var
varn = self.get_varname()
# get a path for axes files
gtaxdir = cmd("echo $GTAXDIR ")
# Read axes
xname = open_gtool(filename).vars[varn].header["AITM1"]
yname = open_gtool(filename).vars[varn].header["AITM2"]
zname = open_gtool(filename).vars[varn].header["AITM3"]
## Stores axes information
dlon = open_gtool(gtaxdir+"/GTAXLOC."+xname).vars[xname][0,0,0,:]
dlat = open_gtool(gtaxdir+"/GTAXLOC."+yname).vars[yname][0,0,0,:]
if zname == "SFC1":
dlev = [0]
else:
dlev = open_gtool(gtaxdir+"/GTAXLOC."+zname).vars[zname][0,0,0,:]
return dlon, dlat, dlev
def get_gtool_header(self):
filename = self.var
varn = self.get_varname()
header = open_gtool(filename).vars[varn].header
return header
def get_gtool_varlongname(self):
header = self.get_gtool_header()
varlongname = header["TITL1"]+header["TITL2"]
return varlongname
def get_gtool_varunit(self):
header = self.get_gtool_header()
varunit = header["UNIT"]
return varunit
def get_slice_and_shape(self, dim, num, original_shape):
axis = self.axes_map[dim]
if num is not None:
shape = len(num)
slice_ = num
else:
shape = original_shape[axis]
slice_ = slice(None)
if dim == "x" and self.xcyclic == True:
shape += 1
return shape, slice_
def gtopen(self):
filename = self.var
varn = self.get_varname()
# Read original data shape
dimensions_string = re.split("[()]", str(open_gtool(filename).vars.values()))[2].split(",")
original_shape = tuple(int(item.strip()) for item in dimensions_string if item.strip())
# Define the axes
num_map = {"x": self.xnum, "y": self.ynum, "z": self.znum, "t": self.tnum}
# Initialize the final shape and the slices to use
final_shape = []
slices = []
# Modify the final shape according to the specified dimensions in self.dim
for dim in self.axes_map.keys():
if dim in self.dim:
shape, slice_ = self.get_slice_and_shape(dim, getattr(self, f"{dim}num"), original_shape)
final_shape.append(shape)
slices.append(slice_)
else:
num = getattr(self, f"{dim}num")
if self.get_length(num) == 1 and self.list2num(num) >= 0:
slices.append(slice(self.list2num(num), self.list2num(num) + 1))
else:
if num is not None:
print("Warning: " + dim + "-axis was set to 0!", "Check your inputs:", dim + "levs =", num, ";", "dim =", self.dim)
slices.append(slice(0, 1))
# Create an empty variable with the final shape
var = np.zeros(final_shape)
final_slices = [slice(None) for _ in var.shape]
var0 = open_gtool(filename).vars[varn][tuple(slices)].copy()
if len(original_shape) != len(final_shape):
for dim in self.axes_map.keys():
if dim not in self.dim:
var0 = np.take(var0, 0, axis=self.axes_map[dim])
# Read original data
if ("x" in self.dim) and (self.xcyclic == True):
var[tuple(final_slices[:-1]+[slice(None, -1)])] = var0
var[...,-1] = var[...,0]
else:
var[tuple(final_slices)] = var0
# Replace undefined values with NaNs
var[var == self.undef] = np.nan
return var
4次元(3次元空間+時間)変数を読み込む場合は、
To read a 4D (3D + t) variable,
gt=ReadGtool()
gt.var="q" # set input file path
gt.dim="xyzt" # set dimension: set xyt to remove z-axis
gt.xcyclic = True # if connect the longitudes = 0 to 360
q = gt.gtopen()
dlon, dlat, dlev = gt.get_gtool_xyzaxis() # get axes
変数の内容を確認するには、
To confirm what the variable is,
varname = gt.get_gtool_varlongname()
varname
'water vapor specific humidity'
変数の単位を確認するには、
To confirm the unit,
varunit = gt.get_gtool_varunit()
varunit
'kg/kg'
それ以外のheader情報は、
For further information,
gt.get_gtool_header()
[00] IDFM : 9010: DSET :jra25.gndg : ITEM :Q :
[03] EDIT1 : : EDIT2 : : EDIT3 : :
[06] EDIT4 : : EDIT5 : : EDIT6 : :
[09] EDIT7 : : EDIT8 : : FNUM : 1:
[12] DNUM : 1: TITL1 :water vapor spec: TITL2 :ific humidity :
[15] UNIT :kg/kg : ETTL1 : : ETTL2 : :
[18] ETTL3 : : ETTL4 : : ETTL5 : :
[21] ETTL6 : : ETTL7 : : ETTL8 : :
[24] TIME : ** NOTE ** : UTIM :HOUR : DATE : ** NOTE ** :
[27] TDUR : ** NOTE ** : AITM1 :GLON128 : ASTR1 : 1:
[30] AEND1 : 128: AITM2 :GGLA64 : ASTR2 : 1:
[33] AEND2 : 64: AITM3 :HETA40 : ASTR3 : 1:
[36] AEND3 : 40: DFMT :UR4 : MISS : -0.9990000E+03:
[39] DMIN : -0.9990000E+03: DMAX : -0.9990000E+03: DIVS : -0.9990000E+03:
[42] DIVL : -0.9990000E+03: STYP : 1: COPTN : :
[45] IOPTN : 0: ROPTN : 0.0000000E+00: DATE1 : ** NOTE ** :
[48] DATE2 : ** NOTE ** : MEMO1 :isotope3:/home/k: MEMO2 :anon/miroc5-iso/:
[51] MEMO3 :out/jra25.gndg/y: MEMO4 :1981 : MEMO5 : :
[54] MEMO6 : : MEMO7 : : MEMO8 : :
[57] MEMO9 : : MEMO10 : : CDATE : ** NOTE ** :
[60] CSIGN :MIROC : MDATE :20230323 192039 : MSIGN :kanon :
[63] SIZE : 327680:
** NOTE **
[24] TIME :[17365476 ... 17373492], (12)
[26] DATE :[19810116 120000 ... 19811216 120000], (12)
[27] TDUR :[744 ... 744], (12)
[47] DATE1 :[19810101 000000 ... 19811201 000000], (12)
[48] DATE2 :[19810201 000000 ... 19820101 000000], (12)
[59] CDATE :[20201223 181831 ... 20201223 203313], (12)
## Annual mean
q_surface_ann = np.mean(q[:,0,:,:],0)
3.2と同じ方法で描画してみましょう。
Let's draw by the method in Section 3.2.
# Reading modules
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
# Drawing
fig = plt.figure() # Setting a plotting region(matplotlib)
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree()) # Using Cartopy
ax.coastlines() # Drawing coastlines
cs = ax.contourf(dlon, dlat, q_surface_ann,transform=ccrs.PlateCarree()) # Drawing the variable
cbar = plt.colorbar(cs, orientation="horizontal") # Drawing a color bar
ax.set_title(varname+" ["+varunit+"]"+"\n(Surface, annual mean)")
plt.show()
読み込む座標を限定したい場合は
gt=ReadGtool()
gt.var="q" # set input file path
gt.dim="xyzt" # set dimension: set xyt to remove z-axis
gt.xcyclic = True # if connect the longitudes = 0 to 360
gt.znum=[3,4,5] # Set grid id to read
q = gt.gtopen()
np.shape(q)
(12, 3, 64, 129)
gt=ReadGtool()
gt.var="q" # set input file path
gt.dim="xyzt" # set dimension: set xyt to remove z-axis
gt.xcyclic = True # if connect the longitudes = 0 to 360
gt.tnum=[3,4,5] # Set grid id to read
q = gt.gtopen()
np.shape(q)
(3, 40, 64, 129)
3次元(水平2次元+時間)変数を読み込む場合は、
To read a 3D (horizontal 2D + t) variable,
gt = ReadGtool()
gt.var = "T2"
gt.dim = "xyt"
gt.xcyclic = True
temp = gt.gtopen()
np.shape(temp)
(12, 64, 129)
3次元変数になっていることが確認できました。
We have confirmed that it is a 3-dimensional variable.
- 2023.4.12 Revised Section 4. (The module was updated by using ChatGPT4.)
- 2022.9.28 Revised Section 4. (Added option for extrapolation of data in latitudinal direction.)
- 2022.4.24 Added Section 4.
- 2021.3.17 Translation by DeepL
- 2021.2.14 Updated (Thanks to Mr. Takano's update of the module, the axis information files in gtool format can now be read, directly.