could xarray read multiple bil files

[shuai-zhou]

I was using data from PRISM; the data format it provides is .bil. I can read individual .bil file using xr.open_rasterio, but could xarray read multiple .bil files using method such as xr.open_mfdataset

[keewis]

you should be able to use xr.open_mfdataset(paths_or_glob, engine="rasterio"), which requires rioxarray to be installed. Also note that xr.open_rasterio is deprecated in favor of xr.open_dataset(path, engine="rasterio").

[shuai-zhou]

Thanks for the help. I used the following code to read the .bil files:

ds = xr.open_mfdataset('./DataExample/bil/PRISM_ppt_stable_4kmD2_20200101_20201231_bil/PRISM_ppt_stable_4kmD2*.bil', combine='nested', compat='override', engine='rasterio')

But the xarray.Dataset does not have a time coordinate; it has a band coordinate which equals 1. I was wondering what is the difference between band and time in this context, and can I retain the time coordinate in the xarray.Datasetwhen reading those .bil files? Thanks, I really appreciate your help.

[keewis]

can you post the repr of two individual files?

If the individual files don't contain the datetime as a dimension you might have to extract it from the filepath in a preprocess function:

def add_time_dimension(ds):
    filepath = ds.encoding["source"]
    # extract the time value from the filename and convert it to datetime.datetime, probably using `dateutil`
    time = ...
    return ds.assign_coords(time=time).expand_dims("time")

xr.open_mfdataset(
    pattern,
    engine="rasterio",
    combine="nested",
    concat_dim="time",
    compat="override",
    preprocess=add_time_dimension,
)

[shuai-zhou]

I stored the files on my local drive, the final names are something like this:

PRISM_ppt_stable_4kmD2_20200101_bil.bil
PRISM_ppt_stable_4kmD2_20200102_bil.bil
...
PRISM_ppt_stable_4kmD2_20201230_bil.bil
PRISM_ppt_stable_4kmD2_20201231_bil.bil

As you mentioned, the .bil files do not have a time dimension, I followed your def and extracted the date for them:

import os
from datetime import datetime

def add_time_dimension(ds):
    filepath = ds.encoding['source']
    # Extrat the time value from the filename and
    # convert it to datetime.datetime using `datetime`
    filename = os.path.basename(filepath)
    time = filename.rsplit('_')[-2].split('.')[0]
    date = datetime.strptime(time, '%Y%m%d').strftime('%Y-%m-%d')
    return ds.assign_coords(time=time).expand_dims('time')

Then, read those files into a dataset (I copied the June, July, August in a folder and only read the summer data):

ds = xr.open_mfdataset('./DataExample/bil/'
                       'PRISM_summer/'
                       'PRISM_ppt_stable_4kmD2*.bil',
                       engine='rasterio',
                       combine='nested',
                       concat_dim='time',
                       preprocess=add_time_dimension)

This returned a dataset with a time dimension:

But when I was trying to aggregate the dataset using avg = ds.groupby('time.season').mean(...), I got the error: AttributeError: 'IndexVariable' object has no attribute 'season'. There is a similar issue reported here.

[keewis]

the issue is this:

    time = filename.rsplit('_')[-2].split('.')[0]
    date = datetime.strptime(time, '%Y%m%d').strftime('%Y-%m-%d')
    return ds.assign_coords(time=time).expand_dims('time')

which should be:

    time = filename.rsplit('_')[-2].split('.')[0]
    date = datetime.strptime(time, '%Y%m%d')
    return ds.assign_coords(time=date).expand_dims('time')

in other words: to be able to use "time.season" as a groupby key the time dimension has to have a datetime dtype

[shuai-zhou]

Right, I forgot that. Thank you so much for your help.