PythonでGTOOLファイルを処理する / Direct Reading & Drawing gtool files by Python

※このページは、かつて 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

[こちら](https://github.com/hyungjun/gtool3/)のioモジュールを利用させていただきます。
`pip install git+https://github.com/hyungjun/gtool3`でインストールできます。
anaconda(miniconda)でパッケージを管理している場合も、インストールしたい環境で上記のpipコマンドを実行して下さい。

We will be using the io module from [here](https://github.com/hyungjun/gtool3/).
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](https://pypi.org/project/gtool3/)に置かれているバージョンは古いため、`pip install gtool3`でインストールすると軸情報ファイルを読み込むことができません。

Note that the version in [PyPI](https://pypi.org/project/gtool3/) is old, so the axis information file cannot be read if you install it with `pip install gtool3`.

早速、モジュールを読み込みましょう。

Let's load the module.


```python
from gtool3 import gtopen as open_gtool
```

データの読み込みは、

Data can be loaded using:


```python
gt = open_gtool("q")
```

表示してみると

Try to display:


```python
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.


```python
gt.vars.keys()
```




    odict_keys(['Q'])




```python
gt.vars.values()
```




    odict_values([Q, (12, 40, 64, 128) : >f4])




```python
gt.vars.items()
```




    odict_items([('Q', Q, (12, 40, 64, 128) : >f4)])



たとえば、以下のように値を取り出せます。

For example, you can retrieve the value as follows


```python
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?


```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:


```python
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.  


```python
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.


```python
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`が空で軸情報ファイルが見当たらない場合は、[こちら](https://github.com/wm-ytakano/gtool3-axis)からダウンロードできます。

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](https://github.com/wm-ytakano/gtool3-axis).


```python
# 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.


```python
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)](https://github.com/yyousuke/matplotlib).


```python
# Reading modules
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
```


```python
# 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.


```python
# 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)
```

![output_63_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/83663/b0c39b77-4f14-461e-abf3-7cba1fdff2e4.png)

`dlon`のいちばん端だけ落としてみましょう。

Let's drop just the very end of `dlon`.


```python
# 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
```

![output_66_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/83663/8b4eaf17-8fcb-ed9f-6520-a962a9b46ea3.png)

描けた！しかし、経度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.


```python
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]
```


```python
# 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
```


![output_70_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/83663/eeef6b20-4c50-6fec-91a3-f7cf1ab54612.png)

## 4.一連の処理をモジュール化する / Modularize the processes

ここでは、gtool形式のファイルを読みこんでから、描画する前までの処理をモジュール化します。

Here, we modularize the process from reading a gtool format file to before drawing.


```python
from gtool3 import gtopen as open_gtool
import numpy as np
import re
import warnings
warnings.simplefilter("ignore")
```


```python
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
```


```python
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,


```python
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,


```python
varname = gt.get_gtool_varlongname()
varname
```




    'water vapor specific humidity'



変数の単位を確認するには、

To confirm the unit,


```python
varunit = gt.get_gtool_varunit()
varunit
```




    'kg/kg'



それ以外のheader情報は、

For further information,


```python
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)




```python
## Annual mean
q_surface_ann = np.mean(q[:,0,:,:],0) 
```

[3.2](#3.2.描画してみる-/-Drawing)と同じ方法で描画してみましょう。

Let's draw by the method in [Section 3.2](#3.2.描画してみる-/-Drawing).


```python
# Reading modules
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
```


```python
# 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()
```
    
![output_93_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/83663/cc9ea00a-a1c3-4edd-e01d-16a20d19b636.png)


読み込む座標を限定したい場合は


```python
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)




```python
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,


```python
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.