Problem loading timestamps correct

[simjoh]

Hi,
I've a net cdf with time format as days since 1899-12-31. But I cannot get xarray to load the dates correct.
the output about the variable and the time variable and the global attributes
ncdump.exe -h thefile.nc
double StdTime(StdTime) ;
StdTime:bounds = "StdTime_bounds" ;
StdTime:long_name = "time" ;
StdTime:standard_name = "time" ;
StdTime:units = "days since 1899-12-30 00:00:00 UTC" ;
// global attributes:
:Conventions = "CF-1.6" ;
:gdal_version_string = "GDAL 3.7.0e (3070046), libnetcdf 4.9.2, HDF5 library version: 1.14.1-2" ;

I've tried multiple different options an combinations. My impression is that the following should be enough?
xr.open_dataset(file,decode_times=False,)['StdTime'].values

But the output is
array([0., 0., 0., 0., 0., 0])

If I open the dataset in gdal I get the times correct.

Any pointers appreciated I'm fairly new to xarray. I'm using xarray=2024.2.0

[andersy005]

@simjoh, can you share here the output of

import xarray as xr 
xr.show_versions()

?

I've tried multiple different options an combinations. My impression is that the following should be enough?
xr.open_dataset(file,decode_times=False,)['StdTime'].values

since it sounds like you want decoded time values, i'm curious as to why you have decode_times=False?

[simjoh]

Thanks for your effort.

ds = xr.open_dataset(file,decode_cf=False)
print(ds['StdTime'].attrs)

{'bounds': 'StdTime_bounds', 'long_name': 'time', 'standard_name': 'time', 'units': 'days since 1899-12-30 00:00:00 UTC'}

print(ds['StdTime'].encoding)

{'dtype': dtype('float64'), 'zlib': True, 'szip': False, 'zstd': False, 'bzip2': False, 'blosc': False, 'shuffle': True, 'complevel': 1, 'fletcher32': False, 'contiguous': False, 'chunksizes': (365,), 'preferred_chunks': {'StdTime': 365}, 'source': '\\UNC_PATH\Temperature_2023.nc', 'original_shape': (365,)}

[spencerkclark]

But I found no way xarray to get it to work with this dataset.

Could you provide more detail? Is there an error message? What are the expected times? Could you show the full output of ds.info()?

With the information provided so far, at least on my end, xarray decodes the times as all 1899-12-30, which seems reasonable, but there may be something else tripping things up:

>>> da = xr.DataArray(
...     [0., 0., 0., 0., 0., 0],
...     dims=["StdTime"],
...     attrs={"units": "days since 1899-12-30 00:00:00 UTC", "bounds": "StdTime_bounds", "long_name": "time", "standard_name": "time"},
...     name="StdTime",
... )
>>> xr.decode_cf(da.to_dataset())
<xarray.Dataset> Size: 48B
Dimensions:  (StdTime: 6)
Coordinates:
  * StdTime  (StdTime) datetime64[ns] 48B 1899-12-30 1899-12-30 ... 1899-12-30
Data variables:
    *empty*

[simjoh]

@spencerkclark Sorry for late response. The output of ds.info() is
'''
xarray.Dataset {
dimensions:
StdTime = 365 ;
x = 676 ;
y = 1560 ;
nbnds = 2 ;

variables:
float64 StdTime(StdTime) ;
StdTime:bounds = StdTime_bounds ;
StdTime:long_name = time ;
StdTime:standard_name = time ;
StdTime:units = days since 1899-12-30 00:00:00 UTC ;
float64 x(x) ;
x:long_name = x-coordinate in cartesian system ;
x:standard_name = projection_x_coordinate ;
x:grid_mapping = transverse_mercator ;
x:units = Meter ;
|S1 transverse_mercator() ;
transverse_mercator:grid_mapping_name = transverse_mercator ;
transverse_mercator:longitude_of_central_meridian = 15.0 ;
transverse_mercator:false_easting = 500000.0 ;
transverse_mercator:false_northing = 0.0 ;
transverse_mercator:latitude_of_projection_origin = 0.0 ;
transverse_mercator:scale_factor_at_central_meridian = 0.9996 ;
transverse_mercator:long_name = CRS definition ;
transverse_mercator:longitude_of_prime_meridian = 0.0 ;
transverse_mercator:semi_major_axis = 6378137.0 ;
transverse_mercator:inverse_flattening = 298.257222101 ;
transverse_mercator:spatial_ref = PROJCS["SWEREF99 TM",GEOGCS["SWEREF99",DATUM["SWEREF99",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6619"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4619"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","3006"]] ;
transverse_mercator:crs_wkt = PROJCS["SWEREF99 TM",GEOGCS["SWEREF99",DATUM["SWEREF99",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6619"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4619"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","3006"]] ;
float64 y(y) ;
y:long_name = y-coordinate in cartesian system ;
y:standard_name = projection_y_coordinate ;
y:grid_mapping = transverse_mercator ;
y:units = Meter ;
float32 MaxTemp(StdTime, y, x) ;
MaxTemp:Description = Daily max temperature 2 m above surface ;
MaxTemp:center = European Center for Medium-Range Weather Forecasts (RSMC) ;
MaxTemp:center_code = 98 ;
MaxTemp:code = 167 ;
MaxTemp:grib_name = 2T ;
MaxTemp:level = SFC ;
MaxTemp:level_code = 1 ;
MaxTemp:level_description = SFC (Ground or water surface) ;
MaxTemp:subcenter_code = 0 ;
MaxTemp:table_code = 128 ;
MaxTemp:grid_mapping = transverse_mercator ;
MaxTemp:units = Celsius ;
float32 MeanTemp(StdTime, y, x) ;
MeanTemp:Description = Daily mean temperature 2 m above surface ;
MeanTemp:center = European Center for Medium-Range Weather Forecasts (RSMC) ;
MeanTemp:center_code = 98 ;
MeanTemp:code = 167 ;
MeanTemp:grib_name = 2T ;
MeanTemp:level = SFC ;
MeanTemp:level_code = 1 ;
MeanTemp:level_description = SFC (Ground or water surface) ;
MeanTemp:subcenter_code = 0 ;
MeanTemp:table_code = 128 ;
MeanTemp:grid_mapping = transverse_mercator ;
MeanTemp:units = Celsius ;
float64 StdTime_bounds(StdTime, nbnds) ;

// global attributes:
:Conventions = CF-1.6 ;
:gdal_version_string = GDAL 3.7.0e (3070046), libnetcdf 4.9.2, HDF5 library version: 1.14.1-2 ;
'''

Regarding the output of the dataset I would assume that xarray gets the same data as gdal.
With gdal I use the following.
'''
gds =gdal.Open(file)
subdatasets = gds.GetSubDatasets()
gdal.Open(subdatasets[2][0]).ReadAsArray()
'''
which produces the output
'''
array([[45291. , 45291.95833333],
[45290. , 45290.999999 ],.
.
.
[44927. , 44927.999999 ]])
'''
Which are two timestamps for each time. I understand that it could be hard for xarray to parse those into timestamps since it is all intervals. But if I use decode=False I could convert it myself from the first column of the numpy array. But all I get is an array of zeros.

[spencerkclark]

Thanks for the further information. Sorry I'm not familiar with GDAL—based on the shape I take it this corresponds with the output of the raw data stored in the StdTime_bounds variable? Does xarray decode that variable properly as a 2D array of times?

[simjoh]

@spencerkclark you got it! Indeed the timeintervals are stored in that variable insted of StdTime as I assumed. Sorry for the troubles. I think I will manage to convert timestamps.

[spencerkclark]

No worries, yeah, it's weird that in this file the values of StdTime are inconsistent with those of StdTime_bounds, but xarray will obviously take the file at its word. Glad to help sort it out.