Python | 处理海洋2C 数据 | 非标准时间格式

写在前面

最近，师弟在用Python读取某海洋2C数据时，突然冒出一个报错：“时间单位的参考日期无效，当前日期 00:00:00.0”。这让我回想起，似乎在很久很久以前，我处理SMAP和Argo数据时也遇到过类似的问题。为了不让未来的师弟师妹们再次被这种问题困扰，我决定在这里记录一下这个坑。

希望未来再看到这条记录时，能少叹气、多微笑，毕竟已经知道怎么解决了嘛！

问题复现

在正常使用xarray读取nc文件时，会出现以下报错：

import xarray as xr
path= r'I:/H2C_*.nc'
ds  = xr.open_dataset(path)

ValueError: invalid reference date for time units: current day 00:00:00.0

ValueError: unable to decode time units 'secs since current day 00:00:00.0' with 'the default calendar'. 

Try opening your dataset with decode_times=False or installing cftime if it is not installed.`

遇到问题就解决问题嘛，上面的报错信息给的也挺详细，包括解决方法，说让你在打开dataset的时候可以加上’decode_times=False’.试试加上确实可以成功打开了

ds  = xr.open_dataset(path,decode_times=False)

此外，我发现可以直接实现netCDF4进行打开，就不会有上面的错误提示。

import netCDF4 as nc
ds  = nc.Dataset(path)

好了，这是遇到的第一个问题，文件打开时出现的错误。

下面是第二个问题：时间转换异常。

下面我希望将数据的时间转化为datetime的格式，方面后续的处理。运行下面的代码进行转换

path= r'I:/H2C_*.nc'
ds  = xr.open_dataset(path,decode_times=False)
time = ds['time'].values
time_c = pd.to_datetime(time)

虽然代码没有报错，但是他的转换结果是有问题的，因为我已知读取的这个文件时间应该是在2021年。所以这是什么原因导致的转换的时间出现异常呢？

问题解释

在一般的netcdf文件中，对于该文件的时间，会给定一个参考的时间标准，然后从标准进行计算，得到相应的时间序列。

在测绘里面有点类似于85国家高程基准的概念。

怎么找到这个时间从参考日期呢，可以直接在time的变量的属性信息里面找到（一般来说都会写在里面，没写的除非特殊说明）

ds.time
Out[139]: 
<xarray.DataArray 'time' (time: 2355)> Size: 19kB
array([6.818707e+08, 6.818707e+08, 6.818707e+08, ..., 6.818738e+08,
       6.818738e+08, 6.818738e+08])
Coordinates:
  * time     (time) float64 19kB 6.819e+08 6.819e+08 ... 6.819e+08 6.819e+08
    lat      (time) float64 19kB ...
    lon      (time) float64 19kB ...
Attributes:
    long_name:           time (sec since 2000-01-01)
    standard_name:       time
    calendar:            gregorian
    tai_utc_difference:  -37.0
    leap_second:         0000-00-00 00:00:00
    units:               seconds since 2000-01-01 00:00:00.0
    comment:             [tai_utc_difference] is the difference between TAI a...

在Attributes-units可以发现，其参考时间为’seconds since 2000-01-01 00:00:00.0’.下面就可以将这个作为基准时间，然后将原本的时间进行转换，主要使用netCDF4.num2date（）函数

path= r'I:/H2C_*.nc'
ds  = xr.open_dataset(path,decode_times=False)
time = ds['time']
time_array = nc.num2date(time,units='sec since 2000-01-01',calendar=time.calendar)
time_array

转换后时间就正常了：

time_array
Out[150]: 
array([cftime.DatetimeGregorian(2021, 8, 10, 0, 31, 55, 472829, has_year_zero=False),
       cftime.DatetimeGregorian(2021, 8, 10, 0, 31, 56, 727732, has_year_zero=False),
       cftime.DatetimeGregorian(2021, 8, 10, 0, 31, 57, 477053, has_year_zero=False),
       ...,
       cftime.DatetimeGregorian(2021, 8, 10, 1, 23, 57, 42110, has_year_zero=False),
       cftime.DatetimeGregorian(2021, 8, 10, 1, 23, 58, 42587, has_year_zero=False),
       cftime.DatetimeGregorian(2021, 8, 10, 1, 23, 59, 43059, has_year_zero=False)],
      dtype=object)

然后再用pandas.to_datetime()进行相应转换

datetime_array = np.array([datetime.datetime(t.year, t.month, t.day, t.hour, t.minute, t.second, t.microsecond) for t in time_array])
time_convert = pd.to_datetime(datetime_array)

time_convert
Out[154]: 
              ['2021-08-10 00:31:55.472829', '2021-08-10 00:31:56.727732',
               '2021-08-10 00:31:57.477053', '2021-08-10 00:31:59.762333',
               '2021-08-10 00:32:06.855991', '2021-08-10 00:32:07.055824',
               '2021-08-10 00:32:10.808600', '2021-08-10 00:32:11.058392',
               '2021-08-10 00:32:12.207441', '2021-08-10 00:32:15.635650',
               ...
               '2021-08-10 01:23:32.429303', '2021-08-10 01:23:33.429862',
               '2021-08-10 01:23:45.886553', '2021-08-10 01:23:46.086657',
               '2021-08-10 01:23:47.937610', '2021-08-10 01:23:55.791510',
               '2021-08-10 01:23:56.041632', '2021-08-10 01:23:57.042110',
               '2021-08-10 01:23:58.042587', '2021-08-10 01:23:59.043059'],
              dtype='datetime64[ns]', length=2355, freq=None)

可以将上述步骤封装为一个函数，方便对于多个文件进行批量处理：

import xarray as xr
import netCDF4 as nc
import numpy as np
import pandas as pd
import datetime
def convert_nc_time_to_datetime(path):
    """
    将NetCDF文件中的时间变量转换为pandas的Datetime对象。

    参数:
        path (str): NetCDF文件的路径。

    返回:
        pd.DatetimeIndex: 转换后的时间数组。
    """
    # 打开NetCDF文件并读取时间变量
    ds = xr.open_dataset(path, decode_times=False)
    time = ds['time']
    
    # 将时间变量转换为datetime对象数组
    time_array = nc.num2date(time, units=time.units, calendar=time.calendar)
    
    # 将datetime对象数组转换为np.datetime64对象
    datetime_array = np.array([datetime.datetime(t.year, t.month, t.day, t.hour, t.minute, t.second, t.microsecond) for t in time_array])
    
    # 将datetime对象数组转换为pandas的DatetimeIndex
    time_convert = pd.to_datetime(datetime_array)
    
    return time_convert
# 使用示例
path = r'I:/H2C_*.nc'
time_convert = convert_nc_time_to_datetime(path)
print(time_convert)

其他场景应用

了解了nc.num2date这个函数后，可以干什么呢？对于两组数据的units不一致时，如果你要将两组数据进行误差比较，就可以用到这个函数，比如我现在有一组SMAP盐度卫星数据以及一组Argo盐度数据，希望将两组数据进行时间匹配以及误差比较。就可以利用这个函数，将Argo的时间标准('days since 1950-01-01 00:00:00')转换为SMAP的时间标准('seconds since 2000-1-1 00:00:00')：

dt=nc.num2date(time,units='days since 1950-01-01 00:00:00')     # Argo 时间标准
argo_dt=nc.date2num(dt,units='seconds since 2000-1-1 00:00:00') #与smap中的时间对应

关于相关函数的使用可以参考引用的链接

https://www.programcreek.com/python/example/89488/netCDF4.num2date

https://unidata.github.io/netcdf4-python/

https://pandas.pydata.org/docs/reference/api/pandas.to_datetime.html