如何基于GRD文件与网格坐标文本创建Xarray Dataset？

1. 解析坐标文件 #

首先处理CANGRD_points_LL.txt，该文件每行对应一个网格点的像素索引（y行、x列）和实际经纬度。我们需要将零散的点转换为与GRD文件维度一致的二维经纬度网格：

import pandas as pd
import numpy as np

# 读取坐标文件，跳过首行表头，用空格分隔列
df = pd.read_csv("CANGRD_points_LL.txt", sep="\s+", skiprows=1, header=None)
df.columns = ["y_idx", "x_idx", "lat", "lon"]

# 生成二维纬度网格：按y行索引分组，提取纬度并按x列索引排序
lat_grid = df.pivot(index="y_idx", columns="x_idx", values="lat").values
# 生成二维经度网格：按y行索引分组，提取经度并按x列索引排序
lon_grid = df.pivot(index="y_idx", columns="x_idx", values="lon").values

pivot操作会自动按索引排序，确保生成的二维数组完全匹配GRD文件的行（y）、列（x）维度。

2. 修正Xarray Dataset的坐标 #

加载GRD文件后，直接替换原有的像素索引坐标为实际经纬度网格：

import glob
import rioxarray as rio
import xarray as xr
import pandas as pd

# 读取所有GRD文件（简化匹配规则，覆盖所有年月文件）
list_files = sorted(glob.glob("t?????.grd"))
# 生成时间序列（1900年1月起，每月一个数据点）
times = pd.date_range("1900/01/01", freq="M", periods=len(list_files))

# 批量加载并合并数据
datarrays = []
for rst, t in zip(list_files, times):
    # 加载单文件，移除多余的band维度
    da = rio.open_rasterio(rst, masked=True).squeeze()
    # 添加时间坐标并扩展维度
    da = da.assign_coords(time=t).expand_dims("time")
    datarrays.append(da)

# 合并所有数据并重命名维度
ds = xr.concat(datarrays, dim="time").rename({"band": "tas", "y": "y_idx", "x": "x_idx"})

# 替换为实际经纬度坐标
ds = ds.assign_coords(
    lat=(["y_idx", "x_idx"], lat_grid),
    lon=(["y_idx", "x_idx"], lon_grid)
)

# 替换维度名为经纬度（可选，让数据集更直观）
ds = ds.swap_dims({"y_idx": "lat", "x_idx": "lon"}).transpose("time", "lat", "lon")

关键注意事项 #

• 维度匹配：如果加载后发现经纬度显示异常，检查GRD文件的y/x维度是否与坐标文件的y_idx/x_idx对应，可交换pivot的index和columns参数调整。
• 缺失值处理：masked=True会保留GRD文件中的原始缺失值掩码，保证后续分析的准确性。
• 性能优化：若文件数量过大，可尝试xr.open_mfdataset结合自定义引擎，但批量加载的方式更直观易调试。

验证结果 #

通过绘图确认坐标是否正确：

ds["tas"].isel(time=0).plot()  # 绘制第一个月的温度异常，检查坐标轴为实际经纬度

内容的提问来源于stack exchange，提问作者User18981898198119