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() 选择整数索引。
• 延迟加载 - 在 open_dataset 中使用 chunks={}，通过 Dask 处理大型文件。
• 保留元数据 - 填充 .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')

分组与重采样

时间序列与空间聚合

# 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 维数学运算与高阶数据分析的桥梁。它处理标签和元数据的能力，让科学代码能够自我说明，并显著提高可靠性。