Create .nc file from multiple tif files : metadata time

Question

I'm trying to create a single .nc file with 3 dimensions : x,y and time. In order to do so I have multiple tif files named by their date : 2013-04-24_Turb.tif for the 24-04-2013. I would like to do this using xarray but I'm a little bit lost. Here is how my tif files are structured :

<xarray.DataArray (band: 1, y: 2093, x: 1104)>
[2310672 values with dtype=float32]
Coordinates:
  * band         (band) int32 1
  * x            (x) float64 7.331e+05 7.332e+05 ... 7.662e+05 7.662e+05
  * y            (y) float64 -3.704e+05 -3.704e+05 ... -4.331e+05 -4.332e+05
    spatial_ref  int32 0
Attributes:
    AREA_OR_POINT:  Area
    scale_factor:   1.0
    add_offset:     0.0

I have used this code :

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

filenames = glob.glob('*.tif')

transformed_filenames = [filename.replace("-", "") for filename in filenames]

def time_index_from_filenames(filenames):
    '''helper function to create a pandas DatetimeIndex
       Filename example: 20150520_0164.tif'''
    return pd.DatetimeIndex([pd.Timestamp(f[:8]) for f in filenames])


time = xr.Variable('time', time_index_from_filenames(transformed_filenames))
chunks = {'x': 1104, 'y': 2093, 'band': 1}
da = xr.concat([rxr.open_rasterio(f, chunks=chunks) for f in filenames], dim=time)
da.to_netcdf('OC3_AR.nc')

But I don't know what to do next

This is the result I got with this code :

<xarray.Dataset>
Dimensions:                        (band: 1, x: 1104, y: 2093, time: 290)
Coordinates:
  * band                           (band) int32 1
  * x                              (x) float64 7.331e+05 7.332e+05 ... 7.662e+05
  * y                              (y) float64 -3.704e+05 ... -4.332e+05
  * time                           (time) datetime64[ns] 2013-04-24 ... 2024-...
Data variables:
    spatial_ref                    int32 ...
    __xarray_dataarray_variable__  (time, band, y, x) float32 ...
PS C:\Users\aurel\OneDrive\Documents\Memoire\Datasets\NDTI_mosaic> 

But when I open the .nc file in QGIS my date are not in the well indicated (see screenshot)

Does someone know how to handle it ?

Answer 1

It looks like your script does what you need it to do. It is just that QGIS doesn't interpret datetime objects in an intuitive way. So what is it that is wrong with it?

You can try using a different software to open your netCDF, for example Panoply should be able to decode the time dimension properly.

I have some quality of life suggestions (assume ds is your dataset):

1. Drop the superfluous dimension with length 1, in this case your band using ds.drop('band') or ds.squeeze(drop = True)
2. Give your data variable a name. Since you use xr.concat this could be done by assigning a name to the dataarray using ds.name = 'NDTI'