はじめに
juliaのnetcdfデータを扱うパッケージにNCDatasetsがある。doiも発行されている（https://doi.org/10.21105/joss.06504）。
pythonのnetCDF4に相当すると思われるが、比較する記事を見つけられなかった。自分も含め、netCDF4よりもxarrayを使う人もいると思うので、ここではごく基本的なxarrayの関数とNCDatasetsの関数をなんとなく対応させることを目指してまとめた。（用語や文法の説明は適当です。気が向けば追記します。）
NCDatasetsのインストール
juliaは導入済みを想定
まずターミナルでjuliaのREPLを起動
$ julia
のち
julia(REPL)
julia> using Pkg
julia> Pkg.add("NCDatasets")
"NCDatasets"を置き換えれば任意のパッケージをインストール可能
xarrayとjuliaの対応
インポート
python3
import xarray as xr
julia
using NCDatasets
netcdfデータの読み込み
python3
ds = xr.open_dataset('file.nc')
julia
ds = Dataset("file.nc")
変数名と次元の確認
- 配列順はpythonはCスタイル、juliaはFスタイル
- ここでは4次元気象データを想定
- pythonで時間、気圧、緯度、経度の順に作成したnetcdfだが、NCDatasetsで読み込むと経度、緯度、気圧、時間の順となっていた（再現性は不明）
python3
print([key for key, val in ds.items()]
print(ds['temperature'].dims)
['temperature']
('time', 'level', 'latitude', 'longitude')
julia
println(keys(ds))
["longitude", "latitude", "level", "time", "temperature"]
変数の取得
python3
temp = ds['temperature'] # どちらでも可
temp = ds.temperature # どちらでも可
julia
temp = ds["temperature"][:]
temp = ds["temperature"] # lazy loadになるらしい
変数の配列サイズの取得
- ちなみに上述の通り配列順は逆になっている
python3
print(ds.shape)
(504, 27, 321, 1440)
julia
println(size(ds))
(1440, 321, 27, 504)
軸情報の取得
- ここで定義した変数は下で使います
python3
lons = ds.longitude # xr.DataArrayを返す
#lons = ds.longitude.values # np.ndarrayを返す
lats = ds.latitude
levs = ds.level
times = ds.time # np.datetime64を保持するxr.DataArray
julia
lons = ds["longitude"][:] # Vectorを返す
lats = ds["latitude"][:]
levs = ds["level"][:]
times = ds["time"][:] # DateTimeを保持するVector
緯度・経度・気圧の指定（幅なし）
- 経度130º、緯度45º、500hPaを指定
- これはxarrayに軍配
- DimensionalData.jlなどのxarray-likeなパッケージはあるが、netcdfの読み込みはできない？（2024/05/29）
- （ここでfindallを使うと選択後のArrayの次元はreduceされない）
python3
temp = temp.sel(level=500, latitude=45, longitude=130)
julia
lon_idx = findfirst(lons .== 130) # 上で定義したlons
lat_idx = findfirst(lats .== 45) # 上で定義したlats
lev_idx = findfirst(levs .== 500) # 上で定義したlevs
temp = temp[lon_idx, lat_idx, lev_idx, :]
時刻の指定（幅なし）
python3
import pandas as pd # 多分、datetimeでもnp.datetime64でもOK
target = pd.Timestamp('2018-02-02 00')
temp = temp.sel(time=target)
julia
using Dates
target = DateTime(2018, 2, 2, 0)
time_idx = findfirst(times .== targed) # 上で定義したtimes
temp = temp[:, :, :, time_idx]
緯度・経度・気圧の指定（幅あり）
- データにより緯度軸は北始まり・南始まりの違いがあるが、xarrayはこの順番に敏感
- juliaの大小比較はデータ順を考える必要ない模様、良い点
python3
temp = temp.sel(level=slice(300, 700),
latitude=slice(50, 35), # データは北→南で格納
longitude=slice(135, 140))
julia
lon_idxs = findall(135. .<= lons .<= 140.)
lat_idxs = findall(35. .<= lats .<= 50.)
lev_idxs = findall(300 .<= levs .<= 700)
temp = temp[lon_idxs, lat_idxs, lev_idxs, :]
時間の指定（幅あり）
python3
temp = temp.sel(time=slice('2018-02-01 00', '2018-02-06 23'))
julia
start_time = DateTime(2018, 2, 1, 0)
end_time = DateTime(2018, 2, 6, 23)
time_idxs = findall(start_time .<= times .<= end_time)
temp = temp[:, :, :, time_idxs]
不要次元の削除
python3
temp = temp.squeeze() # 一度に複数の該当次元をreduce(削除)
julia
tmep = dropdims(temp, dims=1) # 一度に一つの次元のみ
平均
- Statisticsを使う
- 処理後に次元はreduceされないのでdropdimsを使う
python3
temp_mean = temp.mean('time') # shape -> (27, 321, 1440)
julia
using Statistics
temp_mean = mean(temp, dims=4) # size -> (1440, 321, 27, 1)
temp_mean = dropdims(temp, dims=4) # size -> (1440, 321, 27)
その他（他力本願）
- 内挿（stack overflow）
- netcdfデータの保存（本家サイト）
感想
- 普段使っている関数は概ね変換できたと思う
- pythonのnetCDF4やFortranに似てる部分あり、乗り換えコストそれなりに低そう
- xarrayは計算も作図もできるが、NCDatasetsはI/Oに特化している
- 君はDataArray？ndarray？や、君はTimestamp？datetime64？などの混乱を取り除けそうな予感