GrADSでctlファイル無しでsdfopenできるnetcdfファイルの作り方

はじめに

• GrADSは気象・海洋分野で利用される計算・作図ソフト．歴史は長く，利用者は世界的にとても多い．Python (Cartopy)でも作図は可能だが，作図処理に時間がかかるのがネック．GrADSは非常に軽快だが「ctlファイル」というデータのメタ情報（緯度軽度や時刻などの情報）専用のファイルをいちいち作成しなくてはならない．．
• netcdfは自己記述型のバイナリデータでありメタ情報を内包している．GrADSはある規則に従ったnetcdfをctlファイルなしで読み込むことができる．しかし，これまでFortranでnetcdfファイルを作成するのは困難であったが，xarrayというpythonパッケージを利用することでとてもハードルが下がった．
• xarrayを用いてctlファイルなしで利用できるnetcdfファイルの作成法をまとめているサイトは多くなく（自分の知る限りこちらのみ，英語でもヒットしない），また作成には一部注意が必要なので，以下にまとめます．

サンプルのPythonスクリプト

ここでは4次元（時間，気圧面，緯度，経度）の4-byte big_endianのデータを読み込んで変換する．

#!/usr/bin/env python3
import numpy as np
import pandas as pd
import xarray as xr

# coordsに利用数する変数
lons = np.arange(0., 360., 1.25)
lats = np.arange(-90., 90.01, 1.25)
levs = np.array(
    [1000, 975, 950, 925, 900,
     875, 850, 825, 800, 775,
     750, 700, 650, 600, 550,
     500, 450, 400, 350, 300,
     250, 225, 200, 175, 150,
     125, 100, 70, 50, 30, 20,
     10, 7, 5, 3, 2, 1])
times = pd.date_range('2015-01-01 00', '2015-01-31 18', freq='6H')

# 変数の量はデータの次元に対応
it, ip, iy, ix = map(len, [times, levs, lats, lons])

# データの読み込み（気象データJRA-55の6-hourlyジオポテンシャル高度1ヶ月分）
with open('HGT201501.bin') as a:
    data = np.fromfile(a, dtype='>f4').reshape(it, ip, iy, ix)

# DataArrayの作成
da = xr.DataArray(data=data,
                  dims=['time', 'lev', 'lat', 'lon'],
                  coords={'time': times,
                          'lev':('lev', levs, {'units':'millibars'}),
                          'lat':('lat', lats, {'units':'degrees_north'}),
                          'lon':('lon', lons, {'units':'degrees_east'})})

# Datasetの作成とnetcdfの書き出し（多分NetCDF4）
ds = xr.Dataset(data_vars={'hgt': da})  # dする変数名はここで指定
ds.to_netcdf('HGT201501.nc')

GrADSで読み込んで利用してみる．

• sdfopenというコマンドを用いて読み込み
• q fileでデータのメタ情報を確認
• q dimで現在の図の次元の情報を確認
• set t 2でひとつタイムステップを進め，6時間間隔か確認
• set lev 500でlevsの気圧面リストが反映されているか確認
• d hgtで変数名とデータのオーダー（500 hPaだと5500 m前後）を確認

ga-> sdfopen HGT201501.nc
Scanning self-describing file:  HGT201501.nc
SDF file HGT201501.nc is open as file 1
LON set to 0 360
LAT set to -90 90
LEV set to 1000 1000
Time values set: 2015:1:1:0 2015:1:1:0
E set to 1 1
ga-> q file
File 1 :
  Descriptor: HGT201501.nc
  Binary: HGT201501.nc
  Type = Gridded
  Xsize = 288  Ysize = 145  Zsize = 37  Tsize = 124  Esize = 1
  Number of Variables = 1
     hgt  37  t,z,y,x  hgt
ga-> q dim
Default file number is: 1
X is varying   Lon = 0 to 360   X = 1 to 289
Y is varying   Lat = -90 to 90   Y = 1 to 145
Z is fixed     Lev = 1000  Z = 1
T is fixed     Time = 00Z01JAN2015  T = 1
E is fixed     Ens = 1  E = 1
ga-> set t 2
Time values set: 2015:1:1:6 2015:1:1:6
ga-> set lev 500
LEV set to 500 500
ga-> d hgt
Contouring: 4900 to 5900 interval 100

いい感じ！
これで，ctlファイルから解放された世界へ，，

Pythonで読み込む場合

netcdfはxarrayで簡単に読み込める.

import xarray as xr
ds = xr.open_dataset('HGT201501.nc')

要点の簡単な解説

シンプルでエラーが出なければ良いというスタンスであり，深いところまでは理解していないのでご注意．

時間変数について

times = pd.date_range('2015-01-01 00', '2015-01-31 18', freq='6H')

xarrayはpandas.date_rangeをうまいこと利用しているよう．date_rangeは汎用なツールであり時刻の数列をつくのに非常に便利．python標準ライブラリのdatetimeとの互換性もあるし，特にフォーマットを問わないところが好み（例えば，'2015-01-01 00'と書いても'2015/01/01 00'と書いても同様に機能する）．上記引数は左から，初期時刻，最終時刻，時間間隔（6時間）を意味している．時間間隔（freq）はmonthlyデータならfreq='MS'，yearlyデータならfreq='YS'とすると良いと思う（freqの変数リスト）．
ctlファイルのように開始時刻，時間間隔，タイムステップ数を与えたい場合は，

times = pd.date_range('2015-01-01 00', freq='6H', periods=120)

としても良い．

xarrayについて

簡単に一般的な説明を試みる．xarrayにはDatasetとDataArrayというクラスがあり，この2つを抑えればとりあえず大丈夫と思っている．簡単にいうと，Datasetは「弁当箱」で，DataArrayはそこに入れる「おにぎり」と思っていただきたい（？？）．ここではジオポテンシャル高度という一つのデータしか扱っていないが，本来netcdfは複数のデータを格納することができ，それぞれのデータに対応したDataArrayを作る必要がある．つまり，おにぎりは複数あり，それぞれ1種類の具が詰まっている，ということ．最後に，複数のDataArrayに変数名をつけて一つにまとめたものがDatasetである．弁当箱の中でおにぎりの具を区別できるよう名前シールを付ける感じに近い．
（注，Datasetの作成手順は複数あるが，以下では個人的にしっくりきた手順を解説します．）
xarrayについては「中の人」の日本語によるとても詳細な説明があります．
Xarrayを用いたデータ解析 - Qiita

xr.DataArray

da = xr.DataArray(data=data,
                  dims=['time', 'lev', 'lat', 'lon'],
                  coords={'time': times,
                          'lev':('lev', levs, {'units':'millibars'}),
                          'lat':('lat', lats, {'units':'degrees_north'}),
                          'lon':('lon', lons, {'units':'degrees_east'})})

• dataには格納するデータを指定，とりあえずndarrayを入れておけば良いと思っている．DataFrameもイケる気がしているが，Dataset.from_dataframeというメソッドもあるらしい（詳細は知らない）．

• dimsにはデータの変数軸の名前を指定する．ndarrayでは変数軸（axis）は0,1,2,...という整数で管理するが，xarrayではDataFrameのように各軸に名前をつけることができより明示的なコーディングを可能にする．

• coordsには変数軸に対応する変数を指定する．ただし，'units':部はattrsと呼ばれ，とても重要であり，緯度はdegrees_north，経度はdegrees_eastと書かないとGrADSで読み込んでくれない．unitsに関しては以下のサイトを参考にした．
  *make coords for xr.dataarray - Python Tips & Memos by Physical Oceanography & Climate Laboratory in Hokkaido University

• 例には示してないが，attrsというキーワードで変数の単位やその他の情報を辞書で指定できる，例えば，

                  attrs={'units':'geopotential meter'}

としてデータの単位などを記述することができる．（以前エラーが出たトラウマで自分は書かないことが多い．．）

変数の次元が少ない場合

よろしくdimsやcoordsを削れば良いのだが，場合によりnp.newaxis（Noneのエイリアス）を用いて4次元（T, Z, Y, X）にすると良い．より明示的になる（し，unitsの文言を忘れ過去のテンプレを探し彷徨う必要がなくなる）ので．
例えば，2020年1月1日0時UTCの気圧面500 hPaの2次元データとして保存する場合，

da = xr.DataArray(data=data[np.newaxis, np.newaxis, :, :],
                  dims=['time', 'lev', 'lat', 'lon'],
                  coords={'time': pd.date_range('2020-01-01 00', freq='6H', periods=1),
                          'lev':('lev', [500], {'units':'millibars'}),
                          'lat':('lat', lats, {'units':'degrees_north'}),
                          'lon':('lon', lons, {'units':'degrees_east'})})

np.newaxisについては，例えば以下，
np.newaxis Method - DelftStack

xr.Dataset

ds = xr.Dataset(data_vars={'hgt': da})  # dする変数名はここで指定
ds.to_netcdf('HGT201501.nc')

data_varsに変数名と変数（データ名とデータ）を辞書で紐付けてものを与える．

複数の変数（データ）のある場合，例えば，変数軸を共有した3つの異なるデータ（data1, data2, data3）があり，それらを一つのnetcdfへまとめたい場合，

# データ配列，データ名，単位をリスト化
multi_data = [data1, data2, data3]
names = ['wind', 'temperature', 'humidity']
units = ['m/s', 'K', '%']

# 空リストに複数のDataArrayを一旦まとめる
das = []
for data, unit in zip(multi_data, units):
    da = xr.DataArray(data=data,
                      dims=['time', 'lev', 'lat', 'lon'],
                      coords={'time': times,
                              'lev':('lev', levs, {'units':'millibars'}
                              'lat':('lat', lats, {'units':'degrees_north'}),
                              'lon':('lon', lons, {'units':'degrees_east'})},
                      attrs={'units':unit})
    das.append(da)

# DataArrayのリストに変数名を紐付けてDatasetの作成
data_dict = {n: d for n, d in zip(names, das)}  # 辞書の内包表記
ds = xr.Dataset(data_vars=data_dict)

みたいな感じにまとめると良い．個人的に，上記のようにDatasetへごちゃごちゃ記述しないのが好み（多分，Datasetのキーワードを利用した方が合理的な場合もあると思うが）．
辞書の内包表記は以下を参考にした．
Python で 2つのリストを辞書に変換する - DelftStack

よく分かっていない点

• netcdfにはCF Conventionsなるメタデータの記述規則があり，元データをよく確認し，attrsにその記述をした方が丁寧と思う．オフィシャルにデータ共有しないのであればそこまでこだわらなくて良いと思っている．
• Pythonではfloatはデフォルトで64bitなので，データの有効数字によってはハイスペックすぎるかもしれないが，ここではケアしていない．上でも引用しているがこちらの記述は参考になるかも．
• xarrayではctlファイルのtemplate（年月日のパターンマッチで複数のファイルを同一として扱う）的な取り扱いも可能らしい（使う機会がなく調べたことがない）．
• 異なる時刻や異なるデータ格子のDataArrayを同時に格納する際，Datasetに特別な（包括的な）coordsを指定した方が良いと思うが，試したことはない．（あるいはDatasetを先に作成した方がよい場合もありそう）
• 欠損値については，暇ができたら付け足す予定