How to read a stars object with over 32768 bands?

I have a very large dataset consisting of one attribute, simulated daily from 1970 to 2100, defined on a rather fine geographic grid. It has been given to me as a netCDF file, which I would like to read and analyze in an R script. The data is too big to fully fit in memory, so I wrote a script that does the analysis with stars proxy objects and the purrr package. It has worked for similar smaller datasets. However, this dataset seems too big - there are 45956 bands, one for each time step., and it seems like the read_stars() command has an upper limit to how many bands an object can have. This is what my code looks like after loading the proper librairies, where data_path points to a single .nc file:

data_full <- read_stars(data_path, proxy = TRUE)

It returns the following:

Warning message:
In CPL_read_gdal(as.character(x), as.charater(options), as.characters(driver), :
  GDAL Message 1 : Limiting number of bands to 32768 instead of 45956

Then the data is cropped and stops around 2050. I would like to have the full data in the data_full variable. Is is posible to increase the bands limits? Or are there alternative ways of doing this?

Solution

Try setting GDAL_MAX_BAND_COUNT to 65536

Python:

gdal.SetConfigOption('GDAL_MAX_BAND_COUNT',65536)

bash:

export GDAL_MAX_BAND_COUNT=65536