如何读取Berkeley Earth NetCDF文件中年为单位的时间坐标？

问题分析 #

Berkeley Earth的时间坐标采用小数年份格式：整数部分为公元年，小数部分是该年已过时间的占比（例如1850.041667对应1850年1月中旬，1850.125对应1850年2月中旬，间隔为1/12年，属于月度数据）。pd.to_datetime无法直接解析这种格式，且之前用origin='julian'的思路错误——儒略日是从公元前4713年开始的累计天数，和这里的小数年份完全不是同一体系，必然触发范围报错。

解决方案 #

方法1：小数年份转datetime（精准适配格式）

结合numpy和pandas拆分年份与小数占比，用平均年长度（365.25天）适配闰年误差，转换为标准datetime类型：

import xarray as xr
import pandas as pd
import numpy as np

flname = "Land_and_Ocean_LatLon1.nc"
ds = xr.open_dataset(flname)

def decimal_year_to_datetime(decimal_year):
    # 拆分整数年份和小数占比
    year = np.floor(decimal_year).astype(int)
    # 计算该年已过天数（365.25平衡平年/闰年误差）
    day_of_year = (decimal_year - year) * 365.25
    # 转换为datetime对象
    return pd.to_datetime(year, format="%Y") + pd.to_timedelta(day_of_year, unit="D")

# 替换原时间坐标
ds["time"] = decimal_year_to_datetime(ds["time"].values)

# 验证间隔：输出应为约30天的月度间隔
print(np.diff(ds["time"]))

方法2：直接生成月度中点时间（更简洁）

已知数据是月度频率，可直接生成对应年份的月度中旬时间，无需计算小数：

# 生成1850年1月至2022年12月的月度中旬时间序列
ds["time"] = xr.cftime_range(start="1850-01-15", end="2022-12-15", freq="MS") + pd.Timedelta(days=14)

关键提示 #

• 不要用pd.to_datetime(..., unit='D', origin='julian')：儒略日数值远大于1850-2022，完全不匹配当前时间格式。
• 365.25天的计算方式：对于月度精度的数据，足以平衡平年和闰年的天数误差。

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