PythonでGTOOLファイルを処理する / Direct Reading & Drawing gtool files by Python¶
- 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.
目次 / Contents¶
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となっていますね。
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()
gt.vars.values()
gt.vars.items()
たとえば、以下のように値を取り出せます。
For example, you can retrieve the value as follows
print(gt.vars["Q"][:,0,0,0])
2.ヘッダー情報を得る / Reading header information¶
gtoolコマンドだとngtheadで表示してくれますが、pythonだと？
The gtool command will show it in ngthead, but in python?
gt.vars["Q"].header
アイテム名を取得するには、
To get the item name:
gt.vars["Q"].header["ITEM"]
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)
このような軸の情報ファイルを、別途読み込んでくる必要があります。
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)
とすることで、確認できます。
もし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)
以上の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)
緯度方向のみ、本来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
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
変数の単位を確認するには、
To confirm the unit,
varunit = gt.get_gtool_varunit()
varunit
それ以外のheader情報は、
For further information,
gt.get_gtool_header()
## 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)
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次元(水平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)
3次元変数になっていることが確認できました。
We have confirmed that it is a 3-dimensional variable.