Xarray と julia (NCDatasets) の対応

## はじめに
　juliaのnetcdfデータを扱うパッケージに[NCDatasets](https://alexander-barth.github.io/NCDatasets.jl/stable/)がある。doiも発行されている（[https://doi.org/10.21105/joss.06504](https://joss.theoj.org/papers/10.21105/joss.06504)）。
 　pythonのnetCDF4に相当すると思われるが、比較する記事を見つけられなかった。自分も含め、netCDF4よりもxarrayを使う人もいると思うので、ここでは**ごく基本的なxarrayの関数とNCDatasetsの関数をなんとなく対応させること**を目指してまとめた。（用語や文法の説明は適当です。気が向けば追記します。）


## NCDatasetsのインストール
juliaは導入済みを想定
まずターミナルでjuliaのREPLを起動
```bash
$ julia
```
のち
```julia:julia(REPL)
julia> using Pkg
julia> Pkg.add("NCDatasets")
```
"NCDatasets"を置き換えれば任意のパッケージをインストール可能

## xarrayとjuliaの対応

### インポート
```python3:python3
import xarray as xr
```
```julia:julia
using NCDatasets
```

### netcdfデータの読み込み
```python3:python3
ds = xr.open_dataset('file.nc')
```
```julia:julia
ds = Dataset("file.nc")
```

### 変数名と次元の確認
- 配列順はpythonはCスタイル、juliaはFスタイル
- ここでは4次元気象データを想定
- pythonで時間、気圧、緯度、経度の順に作成したnetcdfだが、NCDatasetsで読み込むと経度、緯度、気圧、時間の順となっていた（再現性は不明）
```python3:python3
print([key for key, val in ds.items()]
print(ds['temperature'].dims)
```
['temperature']
('time', 'level', 'latitude', 'longitude')
```julia:julia
println(keys(ds))
```
["longitude", "latitude", "level", "time", "temperature"]

### 変数の取得
```python3:python3
temp = ds['temperature']  # どちらでも可
temp = ds.temperature  # どちらでも可
```
```julia:julia
temp = ds["temperature"][:]
temp = ds["temperature"]  # lazy loadになるらしい
```

### 変数の配列サイズの取得
- ちなみに上述の通り配列順は逆になっている
```python3:python3
print(ds.shape)
```
(504, 27, 321, 1440)
```julia:julia
println(size(ds))
```
(1440, 321, 27, 504)

### 軸情報の取得
- ここで定義した変数は下で使います
```python3: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:julia
lons = ds["longitude"][:]  # Vectorを返す
lats = ds["latitude"][:]
levs = ds["level"][:]
times = ds["time"][:]  # DateTimeを保持するVector
```

### 緯度・経度・気圧の指定（幅なし）
- 経度130º、緯度45º、500hPaを指定
- これはxarrayに軍配
- [DimensionalData.jl](https://juliapackages.com/p/dimensionaldata)などのxarray-likeなパッケージはあるが、netcdfの読み込みはできない？（2024/05/29）
- （ここでfindallを使うと選択後のArrayの次元はreduceされない）
```python3:python3
temp = temp.sel(level=500, latitude=45, longitude=130)
```
```julia: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:python3
import pandas as pd  # 多分、datetimeでもnp.datetime64でもOK
target = pd.Timestamp('2018-02-02 00')
temp = temp.sel(time=target)
```
```julia:julia
using Dates
target = DateTime(2018, 2, 2, 0)
time_idx = findfirst(times .== targed)  # 上で定義したtimes
temp = temp[:, :, :, time_idx]
```

### 緯度・経度・気圧の指定（幅あり）
- データにより緯度軸は北始まり・南始まりの違いがあるが、xarrayはこの順番に敏感
- juliaの大小比較はデータ順を考える必要ない模様、良い点
```python3:python3
temp = temp.sel(level=slice(300, 700),
                latitude=slice(50, 35),  # データは北→南で格納
                longitude=slice(135, 140))
```
```julia: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:python3
temp = temp.sel(time=slice('2018-02-01 00', '2018-02-06 23'))
```
```julia: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:python3
temp = temp.squeeze()  # 一度に複数の該当次元をreduce(削除)
```
```julia:julia
tmep = dropdims(temp, dims=1)  # 一度に一つの次元のみ
```

### 平均
- Statisticsを使う
- 処理後に次元はreduceされないのでdropdimsを使う
```python3:python3
temp_mean = temp.mean('time')  # shape -> (27, 321, 1440)
```
```julia: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](https://stackoverflow.com/questions/68556256/2d-interpolation-in-julia)）
- netcdfデータの保存（[本家サイト](https://alexander-barth.github.io/NCDatasets.jl/stable/#Create-a-netCDF-file)）


## 感想
- 普段使っている関数は概ね変換できたと思う
- pythonのnetCDF4やFortranに似てる部分あり、乗り換えコストそれなりに低そう
- xarrayは計算も作図もできるが、NCDatasetsはI/Oに特化している
- 君はDataArray？ndarray？や、君はTimestamp？datetime64？などの混乱を取り除けそうな予感


## 参考Web文献

- その他の言語の比較について
  - [pythonコードからjuliaへの移植-classとstructに着目して-](https://qiita.com/RhT/items/3e56dbc69dd6747d03ca)
  - [Numpyとの比較で学ぶJuliaのベクトル・行列の記法](https://qiita.com/take_me/items/6933e61714c5359338e0)
  - [pandasで出来る事はどこまでJuliaでできるのか？](https://zenn.dev/otwn/articles/14b90c4db5b974)
  - [Fortranから始めるJulia](https://cometscome.github.io/JuliaFromFortran/build/)
<!--
- 速度について
  - [流体計算の実行速度比較: Fortran, C++, Rust, Python, Julia](https://qiita.com/shigunodo/items/d693dc03323f9a205bb9#傾向)
  - [Fortran, Python, Juliaの計算速度比較:モンテカルロ法を例に](https://kamemori.com/research/fortran/speed_montecarlo_lcgs_ja.html)
  - [JuliaならDataFrameでforループしても十分速い](https://qiita.com/Ken-Kuroki/items/4102f599f49540f80e97)
  - [Julia言語で計算が遅くなった場合の解決法](https://gist.github.com/genkuroki/1ac59bb3e03eac12945d7040d4f98246)
-->