xarray

Xarray - N 次元ラベル付き配列

Xarray は、多次元データに pandas のような操作性を提供します。Pangeo エコシステムの中核であり、NetCDF、GRIB、Zarr 形式を扱ううえで欠かせません。

使用する場面

• 多次元の科学データ（時間、緯度、経度、レベル、アンサンブル）を扱う場合。
• 気候、天気、海洋学のデータセット（NetCDF ファイル）を分析する場合。
• メモリに収まらない大規模なデータセットを扱う場合（Dask との連携を利用）。
• 軸インデックスではなく次元名に基づいて、複雑なブロードキャストやアラインメントを行う場合。
• メタデータ（単位、説明）をデータオブジェクト内に直接保存する場合。
• リモートセンシングや地理空間画像を分析する場合。

リファレンスドキュメント

公式ドキュメント：https://docs.xarray.dev/
チュートリアル：https://tutorial.xarray.dev/
検索パターン：xr.DataArray、xr.Dataset、ds.sel、ds.groupby、ds.resample、xr.open_dataset

基本原則

DataArray と Dataset

構造	説明	たとえ
DataArray	ラベル付きの単一の N 次元配列です。	pandas.Series の N 次元版です。
Dataset	複数の DataArray を格納する辞書のようなコンテナーです。	pandas.DataFrame の N 次元版です。

主な概念

• 次元：軸の名前です（例：x、y、time）。
• 座標：次元に関連付けられた値です（例：実際のタイムスタンプや緯度の値）。
• 属性：任意のメタデータです（例：units='Kelvin'、standard_name='air_temperature'）。

クイックリファレンス

インストール

pip install xarray netCDF4 dask zarr

標準的なインポート

import xarray as xr
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

基本パターン - 作成

import xarray as xr
import numpy as np

# Create a DataArray
data = np.random.rand(4, 3)
times = pd.date_range("2023-01-01", periods=4)
lons = [-120, -110, -100]

da = xr.DataArray(
    data, 
    coords={"time": times, "lon": lons}, 
    dims=("time", "lon"),
    name="temp",
    attrs={"units": "degC"}
)

# Convert to Dataset
ds = da.to_dataset()
print(ds)

重要なルール

✅ 推奨事項

• 名前付き次元を使用する - 整数の軸ではなく、常に dim=('time', 'lat', 'lon') を使用します。
• ラベルで選択する - 座標値には .sel() を、インデックスの整数には .isel() を使用します。
• 遅延読み込み - 大きなファイルを Dask で扱うには、open_dataset で chunks={} を使用します。
• メタデータを保持する - データの自己記述性を確保するため、.attrs に値を設定します。
• アラインメント - ブロードキャストは Xarray に任せます。座標値に基づいてデータが自動的にアラインメントされます。
• アクセサーの活用 - 日時操作には .dt を、文字列操作には .str を使用します。

❌ 非推奨事項

• 整数によるインデックス指定 - data[0, :, 5] は（読みにくく、壊れやすいため）避けます。.isel(time=0, lon=5) を使用します。
• エンコーディングを無視する - NetCDF に保存する際は、圧縮やスケーリングについて ds.encoding を確認します。
• 手動ループ - 時間ステップごとにループせず、.groupby() または .resample() を使用します。
• Dask を忘れない - RAM より大きなデータセットでは、Dask がインストールされ、チャンクが定義されていることを確認します。

アンチパターン（絶対に避ける）

# ❌ BAD: Positional indexing (What is axis 1? Lat or Lon?)
mean_val = ds.variable.mean(axis=1)

# ✅ GOOD: Named dimension reduction (Clear and robust)
mean_val = ds.variable.mean(dim='lat')

# ❌ BAD: Manual time slicing with list comprehensions
# subset = [ds.sel(time=t) for t in my_times if t > '2020-01-01']

# ✅ GOOD: Built-in slicing
subset = ds.sel(time=slice('2020-01-01', '2021-12-31'))

# ❌ BAD: Losing metadata during numpy conversion
raw_data = ds.temp.values # Now it's just a numpy array, units are gone!

# ✅ GOOD: Keep in Xarray as long as possible
processed = ds.temp * 10 # Units and coords are preserved

選択とインデックス指定

sel と isel

# Select by coordinate values
subset = ds.sel(lat=45.0, lon=slice(-100, -80))

# Select by index (integer)
first_step = ds.isel(time=0)

# Nearest neighbor lookup
point = ds.sel(lat=42.1, lon=-71.2, method="nearest")

# Multi-dimensional selection
high_temp_days = ds.where(ds.temp > 30, drop=True)

計算と数学

ブロードキャストとアラインメント

# Xarray aligns automatically by coordinate names
da1 = xr.DataArray([1, 2], coords=[[1, 2]], dims=['x'])
da2 = xr.DataArray([1, 2, 3], coords=[[1, 2, 3]], dims=['y'])

# result is a 2x3 matrix
result = da1 + da2 

# Mathematical operations preserve coordinates
log_temp = np.log(ds.temp)
anomalies = ds.temp - ds.temp.mean(dim='time')

GroupBy とリサンプリング

時系列と空間集約

# Monthly means
monthly = ds.resample(time="1MS").mean()

# Climatology (group by month regardless of year)
climatology = ds.groupby("time.month").mean()

# Calculate anomalies relative to climatology
anomalies = ds.groupby("time.month") - climatology

# Rolling window (Moving average)
rolling_mean = ds.rolling(time=7, center=True).mean()

ファイル I/O（NetCDF、Zarr）

読み込みと書き込み

# Open a single file
ds = xr.open_dataset("weather_data.nc")

# Open multiple files (MFDataset)
ds_all = xr.open_mfdataset("data/*.nc", combine="by_coords", chunks={'time': 100})

# Write to NetCDF
ds.to_netcdf("output.nc")

# Write to Zarr (Cloud optimized)
ds.to_zarr("data.zarr")

プロット

Matplotlib を高水準でラップ

# 1D plot
ds.temp.sel(lat=0, lon=0, method='nearest').plot()

# 2D map
ds.temp.isel(time=0).plot(cmap='RdBu_r', robust=True)

# Faceting (Subplots)
ds.temp.isel(time=slice(0, 4)).plot(col="time", col_wrap=2)

pandas および NumPy との連携

# To Pandas
df = ds.to_dataframe()

# From Pandas
new_ds = xr.Dataset.from_dataframe(df)

# To NumPy (Lose coordinates)
arr = ds.temp.values

# Interoperability
# Xarray objects work in many SciPy/NumPy functions
from scipy.signal import detrend
detrended = xr.apply_ufunc(detrend, ds.temp, input_core_dims=[['time']], output_core_dims=[['time']])

応用：大規模データ向け Dask

メモリに収まらないデータの計算

# Opening with chunks creates Dask arrays
ds = xr.open_dataset("huge_file.nc", chunks={'time': 500, 'lat': 100, 'lon': 100})

# Computation is now lazy
result = ds.temp.mean(dim='time') # Returns immediately

# Trigger computation
final_val = result.compute()

実践的なワークフロー

1. 全球気温偏差ワークフロー

def calculate_temp_anomaly(filepath):
    """Calculate monthly anomalies from NetCDF data."""
    ds = xr.open_dataset(filepath)
    
    # 1. Compute climatology (mean for each month of the year)
    climatology = ds.temp.groupby("time.month").mean("time")
    
    # 2. Subtract climatology from original data
    anomalies = ds.temp.groupby("time.month") - climatology
    
    # 3. Global mean anomaly
    # Weighted by cos(lat) because grid cells get smaller at poles
    weights = np.cos(np.deg2rad(ds.lat))
    weights.name = "weights"
    anom_weighted = anomalies.weighted(weights)
    
    return anom_weighted.mean(("lat", "lon"))

# ts_anomaly = calculate_temp_anomaly("global_temps.nc")

2. 複数モデルのアンサンブル分析

def analyze_ensemble(file_list):
    """Combine multiple model runs into a single dataset with a 'model' dimension."""
    datasets = [xr.open_dataset(f) for f in file_list]
    model_names = ["Model_A", "Model_B", "Model_C"]
    
    # Concatenate along a new dimension
    combined = xr.concat(datasets, dim=pd.Index(model_names, name="model"))
    
    # Calculate ensemble mean and spread
    ens_mean = combined.mean(dim="model")
    ens_std = combined.std(dim="model")
    
    return ens_mean, ens_std

3. 衛星画像処理（NDVI）

def calculate_ndvi(ds):
    """Calculate NDVI from Red and NIR bands in an Xarray Dataset."""
    # NDVI = (NIR - Red) / (NIR + Red)
    red = ds.sel(band='red')
    nir = ds.sel(band='nir')
    
    ndvi = (nir - red) / (nir + red)
    ndvi.attrs['long_name'] = "Normalized Difference Vegetation Index"
    
    return ndvi

パフォーマンスの最適化

チャンク化の戦略

# ❌ Problem: Small chunks lead to massive overhead
# ds = ds.chunk({'time': 1, 'lat': 1, 'lon': 1})

# ✅ Solution: Aim for 10MB - 100MB per chunk
ds = ds.chunk({'time': -1, 'lat': 100, 'lon': 100})

apply_ufunc によるベクトル化

# Wrap a custom numpy function to work on Xarray objects efficiently
def my_complex_stat(x):
    return np.median(x) * np.std(x)

result = xr.apply_ufunc(
    my_complex_stat, 
    ds.temp,
    input_core_dims=[['time']], # The dimension to map the function over
    vectorize=True,
    dask="parallelized"
)

よくある落とし穴と解決策

座標の不一致

# ❌ Problem: DataArrays don't align due to floating point jitter in lat/lon
# ✅ Solution: Use .interp_like() or .reindex_like()
ds2_aligned = ds2.interp_like(ds1)

values によるメモリリーク

# ❌ Problem: Calling .values on a huge Dask-backed array crashes the machine
# ✅ Solution: Use .compute() or subset first
subset_val = ds.temp.isel(time=0).values # This is safer

スライスの問題（開始／終了）

# ❌ Problem: slice(10, 0) returns empty because order is wrong
# ✅ Solution: Check if your index is ascending or descending
# ds.sortby('lat').sel(lat=slice(-90, 90))

Xarray は、生の N 次元演算と高水準のデータ分析をつなぐ橋渡し役です。ラベルやメタデータを扱えるため、科学技術コードが自己記述的になり、信頼性も大幅に向上します。